diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index 8eaf1f8a..6c924f20 100644 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -11,12 +11,12 @@ testyoffset 0.0 ! y-axis offset of crater center fro tallyonly F ! Tally the craters without generating any craters testtally F quasimc T ! MC run constrained by non-MC 'real' craters given in a list -realcraterlist craterlist.in ! list of 'real' craters for Quasi-MC runs +realcraterlist qmctest.in ! list of 'real' craters for Quasi-MC runs ! IDL driver in uts -interval 628 +interval 0.1 numintervals 1 ! Total number of intervals (total time = interval * numintervals) <--when runtype is 'single' restart F ! Restart a previous run impfile NPFextrap.dat ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1)) diff --git a/examples/global-lunar-bombardment/temp.in b/examples/global-lunar-bombardment/temp.in new file mode 100644 index 00000000..328c21ed --- /dev/null +++ b/examples/global-lunar-bombardment/temp.in @@ -0,0 +1,76 @@ +! CTEM Input file + + +! Testing input. These are used to perform non-Monte Carlo tests. +testflag T! F ! Set to T to create a single crater with user-defined impactor properties +testimp 10 ! 15000.0 ! Diameter of test impactor (m) +testvel 18.3e3 ! Velocity of test crater (m/s) +testang 45.0 ! Impact angle of test crater (deg) - Default 90.0 +testxoffset 0.0 ! x-axis offset of crater center from grid center (m) - Default 0.0 +testyoffset 0.0 ! y-axis offset of crater center from grid center (m) - Default 0.0 +tallyonly F ! Tally the craters without generating any craters +testtally F +quasimc F! T ! MC run constrained by non-MC 'real' craters given in a list +realcraterlist qmctest.in ! list of 'real' craters for Quasi-MC runs + + + +! IDL driver in uts +interval 1 ! 0.1 +numintervals 1 ! 1 ! Total number of interval 1 !s (total time = interval 1 ! * numintervals 1 !) <--when runtype is 'single' +restart F ! Restart a previous run +impfile NPFextrap.dat ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1)) +popupconsole F ! Pop up console window every output interval 1 ! +saveshaded F ! Output shaded relief images +saverego T ! Output regolith map images +savepres F ! Output simplified console display images (presentation-compatible images) +savetruelist T ! Save the true cumulative crater distribution for each interval 1 ! (large file size) +sfdcompare LOLASethCraterCatalogv8gt20-binned.dat ! File name for the SFD comparison file used in the console display +shadedminh -85.0 ! Minimum height for shaded relief map (m) (Default - automatic) +shadedmaxh 85.0 ! Maximum height for shaded relief map (m) (Default - automatic) +runtype single ! Run type: options are normal / statistical + ! single: craters accumulate in successive interval 1 !s + ! statistical: surface is reset between interval 1 !s + +! CTEM required inputs +seed 76535 ! Random number generator seed +gridsize 2000 ! Size of grid in pixels +numlayers 10 ! Number of perched layers +pix 3.08e3 ! Pixel size (m) +mat rock ! Material (rock or ice) +! Bedrock scaling parameters +mu_b 0.55e0 ! Experimentally derived parameter for bedrock crater scaling law +kv_b 0.20e0 ! Experimentally derived parameter for bedrock crater scaling law +trho_b 2250.0e0 ! Target density (bedrock) (kg/m**3) +ybar_b 0.0e6 ! Bedrock strength (Pa) +! Regolith scaling parameters +mu_r 0.55e0 ! Experimentally derived parameter for regolith crater scaling law +kv_r 0.20e0 ! Experimentally derived parameter for regolith crater scaling law +trho_r 2250.0e0 ! Target density (regolith) (kg/m**3) +ybar_r 0.00e6 ! Regolith strength (Pa) +! Body parameters +gaccel 1.62e0 ! Gravitational acceleration at target (m/s**2) +trad 1737.35e3 ! Target radius (m) +prho 2500.0e0 ! Projectile density (kg/m**3) +sfdfile production.dat ! Impactor SFD file (col 1: Dimp (m), col 2: ! impactors > D +velfile lunar-MBA-impactor-velocities.dat ! Impactor velocity dist file + +! Seismic shaking input (required if seismic shaking is set to T, otherwise ignored) +doseismic F ! Perform seismic shaking calculations with each impact - Default F + +! Optional inputF These have internally set default values that work reasonable well. Comment them out with +deplimit 9e99 ! Depth limit for craters (m) - Default is to ignore. +maxcrat 3.15e-2 ! Fraction of gridsize that maximum crater can be - Default 1.0 +killatmaxcrater F ! Stop the run if a crater larger than the maximum is produced - Default F +basinimp 35.0e3 ! Size of impactor to switch to lunar basin scaling law - Default is to ignore +docollapse T ! Do slope collapse - Default T +dosoftening T ! Do ejecta softening - Default T +doangle T ! Vary the impact angle. Set to F to have only vertical impacts - Default T +Kd1 0.0001 +psi 2.000 +fe 4.00 +ejecta_truncation 4.0 +doregotrack T +dorays F +superdomain F +dorealistic F diff --git a/src/Makefile.am b/src/Makefile.am index 59c00bce..de331416 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -11,7 +11,7 @@ AM_FCFLAGS = $(IDEBUG) #ifort debug flags #gfortran optimized flags -#AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace +AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace #gfortran debug flags #AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form -fsanitize-address-use-after-scope -fstack-check -fsanitize=bounds-strict -fsanitize=undefined -fsanitize=signed-integer-overflow -fsanitize=object-size -fstack-protector-all diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index e3087269..ddd5e32b 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -151,6 +151,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt oldpbarpos = 0 do while (icrater < ntotcrat) makecrater = .true. + domain%currentqmc = .false. timestamp_old = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=SP) icrater = icrater + 1 crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=SP) @@ -158,6 +159,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt !if in quasiMC mode: check to see if it's time for a real crater if (user%doquasimc) then if ((user%rctime > timestamp_old) .and. (user%rctime < crater%timestamp)) then + domain%currentqmc = .true. write(message,*) "Real @ ", crater%timestamp call io_updatePbar(message) user%testflag = .true. @@ -189,6 +191,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt if (user%doquasimc) then if (crater%timestamp > user%rctime) then user%testflag = .false. + domain%nqmc = domain%rccount domain%rccount = domain%rccount + 1 if (domain%rccount > domain%rcnum) then write(message,*) "Real crater list complete." diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index ab42f236..50fc657e 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -155,8 +155,10 @@ module module_globals real(DP) :: ejbres ! Ejecta blanket lookup table resolution integer(I4B) :: pnum ! size of production function array integer(I4B) :: vnum ! size of velocity distribution array - integer(I4B) :: rcnum ! size of real crater list array for quasi-MC - integer(I4B) :: rccount ! Current quasimc crater + integer(I4B) :: rcnum ! size of real crater list array for quasi-MC + integer(I4B) :: rccount ! quasimc crater for iterating through the list of craters + integer(I4B) :: nqmc ! current Quasi-MC crater for regolith distribution + logical :: currentqmc ! is the current crater a Quasi-MC crater? real(DP) :: vescsq ! Escape velocity at target integer(I4B) :: vlo ! Index of lowest valid velocity in the velocity distribution file integer(I4B) :: vhi ! Index of highest valid velocity in the velocity distribution file @@ -302,6 +304,7 @@ module module_globals character(*),parameter :: DATFILE = 'ctem.dat' character(*),parameter :: MASSFILE = 'impactmass.dat' character(*),parameter :: RCFILE = 'craterlist.dat' +character(*),parameter :: MDFILE = 'surface_meltdist.dat' ! Global variables integer(I4B),parameter :: PBCLIM = 1 ! periodic boundary condition limit diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90 index f07e3681..c574b7cb 100644 --- a/src/io/io_write_regotrack.f90 +++ b/src/io/io_write_regotrack.f90 @@ -16,7 +16,7 @@ ! Notes : ! !********************************************************************************************************************************** -subroutine io_write_regotrack(user,surf) +subroutine io_write_regotrack(user,surf,domain) use module_globals use module_io, EXCEPT_THIS_ONE => io_write_regotrack implicit none @@ -24,6 +24,7 @@ subroutine io_write_regotrack(user,surf) ! Arguments type(usertype),intent(in) :: user type(surftype),dimension(:,:),intent(in) :: surf + type(domaintype),intent(in) :: domain ! Regotrack Internals integer(I4B) :: i,j,k @@ -32,11 +33,12 @@ subroutine io_write_regotrack(user,surf) integer(I4B), parameter :: FREGO = 11 integer(I4B), parameter :: FCOMP = 12 integer(I4B), parameter :: FAGE = 13 + integer(I4B), parameter :: FMD = 14 !type(regolisttype),pointer :: current => null() type(regodatatype),dimension(:),allocatable :: current integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num real(DP),dimension(:),allocatable :: meltfrac, thickness, comp - real(SP),dimension(:,:),allocatable :: age + real(SP),dimension(:,:),allocatable :: age, meltdist integer(kind=8) :: recsize real(DP) :: dtmp real(SP) :: stmp @@ -49,6 +51,7 @@ subroutine io_write_regotrack(user,surf) open(FREGO,file=REGOFILE,status='replace',form='unformatted') open(FCOMP,file=COMPFILE,status='replace',form='unformatted') open(FAGE,file=AGEFILE,status='replace',form='unformatted') + open(FMD,file=MDFILE,status='replace',form='unformatted') ! First pass to get stack numbers stacks_num(:,:) = 0 @@ -71,7 +74,7 @@ subroutine io_write_regotrack(user,surf) do i=1,user%gridsize !current => surf(i,j)%regolayer N = stacks_num(i,j) - allocate(meltfrac(N),thickness(N),comp(N),age(MAXAGEBINS,N)) + allocate(meltfrac(N),thickness(N),comp(N),age(MAXAGEBINS,N),meltdist(domain%rcnum,N)) allocate(current,source=surf(i,j)%regolayer) do k=1,N ! meltfrac(k) = current%regodata%meltfrac @@ -83,6 +86,9 @@ subroutine io_write_regotrack(user,surf) thickness(k) = current(k)%thickness comp(k) = current(k)%comp age(:,k) = current(k)%age(:) + write(*,*) size(current(k)%age) + write(*,*) size(current(k)%meltdist) + meltdist(:,k) = current(k)%meltdist(:) end do deallocate(current) @@ -90,13 +96,15 @@ subroutine io_write_regotrack(user,surf) write(FREGO) thickness(:) write(FCOMP) comp(:) write(FAGE) age(:,:) - deallocate(meltfrac,thickness,comp,age) + write(FMD) meltdist(:,:) + deallocate(meltfrac,thickness,comp,age,meltdist) end do end do close(FMELT) close(FREGO) close(FCOMP) close(FAGE) + close(FMD) recsize = sizeof(itmp) * user%gridsize * user%gridsize open(LUN,file=STACKNUMFILE,status='replace',form='unformatted',recl=recsize,access='direct') diff --git a/src/io/io_write_surf.f90 b/src/io/io_write_surf.f90 index 45686009..cd11acc2 100644 --- a/src/io/io_write_surf.f90 +++ b/src/io/io_write_surf.f90 @@ -16,7 +16,7 @@ ! Notes : ! !********************************************************************************************************************************** -subroutine io_write_surf(user,surf) +subroutine io_write_surf(user,surf,domain) use module_globals use module_io, EXCEPT_THIS_ONE => io_write_surf implicit none @@ -24,6 +24,7 @@ subroutine io_write_surf(user,surf) ! Arguments type(usertype),intent(in) :: user type(surftype),dimension(:,:),intent(in) :: surf + type(domaintype),intent(in) :: domain ! Internals integer(I4B) :: i @@ -67,7 +68,7 @@ subroutine io_write_surf(user,surf) end do close(LUN) - if (user%doregotrack) call io_write_regotrack(user,surf) + if (user%doregotrack) call io_write_regotrack(user,surf,domain) ! if (user%docrustal_thinning) then ! recsize = sizeof(itmp) * user%gridsize * user%gridsize diff --git a/src/io/module_io.f90 b/src/io/module_io.f90 index 67e18ec3..3ea2873a 100644 --- a/src/io/module_io.f90 +++ b/src/io/module_io.f90 @@ -69,20 +69,22 @@ end subroutine io_read_regotrack end interface interface - subroutine io_write_surf(user,surf) + subroutine io_write_surf(user,surf,domain) use module_globals implicit none type(usertype),intent(in) :: user type(surftype),dimension(:,:),intent(in) :: surf + type(domaintype),intent(in) :: domain end subroutine io_write_surf end interface interface - subroutine io_write_regotrack(user,surf) + subroutine io_write_regotrack(user,surf,domain) use module_globals implicit none type(usertype),intent(in) :: user - type(surftype),dimension(:,:),intent(in) :: surf + type(surftype),dimension(:,:),intent(in) :: surf + type(domaintype),intent(in) :: domain end subroutine io_write_regotrack end interface diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90 index bd58eff7..5cdc3d08 100644 --- a/src/main/CTEM.f90 +++ b/src/main/CTEM.f90 @@ -147,7 +147,7 @@ program CTEM call io_write_tally(truedist,truelist(:,1:ntrue),obsdist,obslist,oposlist,depthdiam,degradation_state) if (.not.user%tallyonly) then write(*,*) "Writing surface files" - call io_write_surf(user,surf) + call io_write_surf(user,surf,domain) end if if (user%testflag) then ! Draw a profile across the crater diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90 index d9941907..09fb2349 100644 --- a/src/regolith/module_regolith.f90 +++ b/src/regolith/module_regolith.f90 @@ -246,10 +246,11 @@ end subroutine regolith_depth_model end interface interface - subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt) + subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt) use module_globals type(usertype),intent(in) :: user type(cratertype),intent(in) :: crater + type(domaintype),intent(in) :: domain real(DP),intent(in) :: age real(DP),intent(in) :: age_resolution real(DP),intent(in) :: ebh diff --git a/src/regolith/regolith_melt_glass.f90 b/src/regolith/regolith_melt_glass.f90 index 94b0b5e3..07f21be9 100644 --- a/src/regolith/regolith_melt_glass.f90 +++ b/src/regolith/regolith_melt_glass.f90 @@ -55,7 +55,7 @@ ! Notes : ! !********************************************************************************************************************************** -subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt) +subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt) use module_globals use module_util use module_regolith, EXCEPT_THIS_ONE => regolith_melt_glass @@ -64,6 +64,7 @@ subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad, ! Arguments type(usertype),intent(in) :: user type(cratertype),intent(in) :: crater + type(domaintype),intent(in) :: domain real(DP),intent(in) :: age real(DP),intent(in) :: age_resolution real(DP),intent(in) :: ebh @@ -164,6 +165,7 @@ subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad, volm1 = regolith_streamtube_volume_func(eradi,0.0_DP,xmints,deltar) melt = volm1 - volv1 newlayer%meltfrac = melt/(vst-volv1) + newlayer%meltdist(domain%nqmc) = newlayer%meltfrac end if n_age = max(ceiling(age / age_resolution), 1) diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index ae9a1c15..1de21914 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -226,7 +226,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, ! Purpose 2: Once we have the size information of a stream tube, we can ! calculate the distal melt: the precursor of glass spherules within a ! stream tube. The result is contained in a linked list "newlayer". - call regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints, volm) + call regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints, volm) ! if (eradc>rm) then ! write(*,*) 'eradc > rm!' ! write(*,*) ebh, exp(ejb(k)%thick) diff --git a/src/regolith/regolith_superdomain.f90 b/src/regolith/regolith_superdomain.f90 index ba48f482..3c3d76ce 100644 --- a/src/regolith/regolith_superdomain.f90 +++ b/src/regolith/regolith_superdomain.f90 @@ -60,7 +60,7 @@ subroutine regolith_superdomain(user,crater,domain,regolayer,ejdistribution,xpi, erad = cvpg**(user%mu_b) * (lrad**(-0.5_DP * user%mu_b)) * (crater%rad)**(user%mu_b * 0.5_DP + 1.0_DP) vej = cvpg * sqrt(user%gaccel * crater%grad) * (erad / crater%grad)**(-1.0_DP / user%mu_b) !equation 18 in Richardson 2009 lrad = ( vej **2 ) / user%gaccel !assume ejection angle is 45 degree. - call regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,erad,lrad,deltar,newlayer,xmints,melt) + call regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad,lrad,deltar,newlayer,xmints,melt) call util_push_array(regolayer,newlayer) return