diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index 14f15369..7b6750a7 100644 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -3,20 +3,20 @@ ! Testing input. These are used to perform non-Monte Carlo tests. testflag F ! Set to T to create a single crater with user-defined impactor properties -testimp 15000.0 ! Diameter of test impactor (m) +testimp 100000 ! 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 +testxoffset 0 ! x-axis offset of crater center from grid center (m) - Default 0.0 +testyoffset 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 T ! MC run constrained by non-MC 'real' craters given in a list -realcraterlist craterlist.in ! list of 'real' craters for Quasi-MC runs +quasimc F ! MC run constrained by non-MC 'real' craters given in a list +realcraterlist qmctest2.in ! list of 'real' craters for Quasi-MC runs ! IDL driver in uts -interval 628.0 +interval 10.0 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)) @@ -33,10 +33,10 @@ runtype single ! Run type: options are normal / ! statistical: surface is reset between intervals ! CTEM required inputs -seed 76535 ! Random number generator seed +seed 765 ! Random number generator seed gridsize 500 ! Size of grid in pixels numlayers 10 ! Number of perched layers -pix 12.32e3 ! Pixel size (m) +pix 12.32e3 ! Pixel size (m) mat rock ! Material (rock or ice) ! Bedrock scaling parameters mu_b 0.55e0 ! Experimentally derived parameter for bedrock crater scaling law @@ -60,7 +60,7 @@ doseismic F ! Perform seismic shaking calcul ! 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 +maxcrat 0.15 ! Fraction of gridsize that maximum crater can be - Default 1.0 (0.15 for full MC; 3.15e-2 for Quasi-MC) 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 @@ -71,7 +71,7 @@ psi 2.000 fe 4.00 ejecta_truncation 24.0 doregotrack T -dorays T +dorays F superdomain F dorealistic F domixing T diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 5d031596..118cdd5d 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -94,6 +94,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt real(SP),dimension(60) :: agetot type(regolisttype),pointer :: current => null() real(DP) :: age_resolution + integer(I4B) :: age_counter nmixingtimes = 0 @@ -143,8 +144,9 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! Reset age clock = 0.0_DP finterval = 1.0_DP / real(ntotcrat,kind=DP) - age = user%interval + age = user%interval * user%numintervals age_resolution = age / real(MAXAGEBINS) + domain%age_counter = 1 ! Reset coverage map domain%tallycoverage = 0 @@ -160,6 +162,9 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt icrater = icrater + 1 crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=SP) pbarpos = nint(real(icrater) / real(ntotcrat) * PBARRES) + if (crater%timestamp > (domain%age_counter*age_resolution)) then + domain%age_counter = domain%age_counter + 1 + end if !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 diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 3dda6811..e105c279 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -159,6 +159,7 @@ module module_globals integer(I4B) :: vnum ! size of velocity distribution array 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(I2B) :: age_counter ! Current age bin for the production of new melt integer(I4B) :: local ! integer(I4B) :: nqmc ! current Quasi-MC crater for regolith distribution logical :: currentqmc ! is the current crater a Quasi-MC crater? diff --git a/src/regolith/regolith_interior.f90 b/src/regolith/regolith_interior.f90 index ef4046ad..173b7f18 100644 --- a/src/regolith/regolith_interior.f90 +++ b/src/regolith/regolith_interior.f90 @@ -74,6 +74,7 @@ subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltshee newlayer%distvol(domain%nqmc) = newlayer%meltvolume else newlayer%distvol(1+domain%rcnum) = newlayer%meltvolume + newlayer%age(domain%age_counter) = newlayer%meltvolume end if call util_push_array(surf(xpi,ypi)%regolayer,newlayer) end do diff --git a/src/regolith/regolith_melt_glass.f90 b/src/regolith/regolith_melt_glass.f90 index df47d1c4..10b69fd6 100644 --- a/src/regolith/regolith_melt_glass.f90 +++ b/src/regolith/regolith_melt_glass.f90 @@ -178,14 +178,17 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad newlayer%distvol(:) = 0.0_SP if(domain%currentqmc) then newlayer%distvol(domain%nqmc) = newlayer%meltvolume + else + newlayer%distvol(1+domain%rcnum) = newlayer%meltvolume + newlayer%age(domain%age_counter) = newlayer%meltvolume end if - n_age = max(ceiling(age / age_resolution), 1) - if (lrad >= RAD_GP * crater%rad) then - newlayer%age(n_age) = melt / (user%pix * user%pix) - else - newlayer%age(n_age) = 0.0_SP - end if + ! n_age = max(ceiling(age / age_resolution), 1) + ! if (lrad >= RAD_GP * crater%rad) then + ! newlayer%age(n_age) = melt / (user%pix * user%pix) + ! else + ! newlayer%age(n_age) = 0.0_SP + ! end if return end subroutine regolith_melt_glass diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index 1fe3584c..77e2f74c 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -269,7 +269,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, tots = totseb newlayer%thickness = ebh newlayer%comp = min(totmare/tots, 1.0_DP) - newlayer%age(:) = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP) + !newlayer%age(:) = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP) else rbody = sqrt(xints(2)**2 + yints(2)**2) @@ -313,7 +313,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%thickness = ebh newlayer%comp = min(totmare/tots, 1.0_DP) newlayer%age(:) = newlayer%age(:) + age_collector(:) - newlayer%age(:) = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP) + !newlayer%age(:) = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP) ! if (newlayer%meltfrac > 1.0_DP) then ! write(*,*) "Melt fraction >1! (Traverse)", xpi,ypi,crater%timestamp,crater%fcrat,crater%xlpx,crater%ylpx,& ! newlayer%meltvolume, newlayer%totvolume, newlayer%ejm, newlayer%ejmf, totvol diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index 9460260c..3057cde7 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -66,11 +66,11 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector tots = tots + vsgly totmare = totmare + vsgly * current(M)%comp recyratio = vsgly / (user%pix**2) /current(M)%thickness - age_collector(:) = age_collector(:) + current(M)%age(:) * recyratio ratio = vsgly / current(M)%totvolume if (ratio > 1) then ratio = 1.0_DP end if + age_collector(:) = age_collector(:) + (current(M)%age(:) * ratio) headmeltvol = current(M)%meltvolume * ratio meltinejecta = meltinejecta + headmeltvol distvol(:) = distvol(:) + (current(M)%distvol(:) * ratio) @@ -85,11 +85,11 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector tothead = tothead + vhead * vratio totmarehead = totmarehead + vhead * vratio * current(N)%comp recyratio = vhead * vratio / (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio ratio = (vhead * vratio) / current(N)%totvolume if (ratio > 1) then ratio = 1.0_DP end if + age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio) headmeltvol = current(N)%meltvolume * ratio meltinejecta = meltinejecta + headmeltvol distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) @@ -103,11 +103,11 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector totmarehead = totmarehead + (vsgly-tothead) * current(N)%comp tothead = vsgly recyratio = (vsgly - tothead) / (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio ratio = (vsgly-tothead) / current(N)%totvolume if (ratio > 1) then ratio = 1.0_DP end if + age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio) headmeltvol = current(N)%meltvolume * ratio meltinejecta = meltinejecta + headmeltvol distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90 index a834e8c4..cf72a475 100644 --- a/src/regolith/regolith_streamtube_lineseg.f90 +++ b/src/regolith/regolith_streamtube_lineseg.f90 @@ -83,12 +83,12 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) ebh_recyl = (1.0_DP - (newlayer%meltvolume/newlayer%totvolume)) * newlayer%thickness - vsh / user%pix**2 recyratio = max(ebh_recyl,0.0_DP) / current(N-1)%thickness - age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio ratio = max(vsgly-vsh,0.0_DP) / current(N-1)%totvolume if (ratio > 1) then ratio = 1.0_DP end if - vol = vol + sum(current(N-1)%age(:)) * recyratio + age_collector(:) = age_collector(:) + (current(N-1)%age(:) * ratio) + !vol = vol + sum(current(N-1)%age(:)) * recyratio linmelt = current(N-1)%meltvolume * ratio meltinejecta = meltinejecta + linmelt distvol(:) = distvol(:) + (current(N-1)%distvol(:) * ratio) @@ -97,12 +97,12 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vseg = regolith_streamtube_volume_func(eradi,xmints,rip1,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N-1)%thickness - age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio ratio = max(vseg-vsh,0.0_DP) / current(N-1)%totvolume if (ratio > 1) then ratio = 1.0_DP end if - vol = vol + sum(current(N-1)%age(:)) * recyratio + age_collector(:) = age_collector(:) + (current(N-1)%age(:) * ratio) + !vol = vol + sum(current(N-1)%age(:)) * recyratio linmelt = current(N-1)%meltvolume * ratio meltinejecta = meltinejecta + linmelt distvol(:) = distvol(:) + (current(N-1)%distvol(:) * ratio) @@ -131,12 +131,12 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vseg = regolith_streamtube_volume_func(eradi,max(xmints,rstart),rend,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rstart,rend) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio + age_collector(:) = age_collector(:) + current(N)%age(:) * ratio ratio = max((vseg-vsh),0.0_DP) / current(N)%totvolume if (ratio > 1) then ratio = 1.0_DP end if - vol = vol + sum(current(N)%age(:)) * recyratio + !vol = vol + sum(current(N)%age(:)) * recyratio linmelt = current(N)%meltvolume * ratio meltinejecta = meltinejecta + linmelt distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) @@ -148,12 +148,12 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vseg = regolith_streamtube_volume_func(eradi,max(xmints,rend),rstart,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rend,rstart) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio ratio = max((vseg-vsh),0.0_DP) / current(N)%totvolume if (ratio > 1) then ratio = 1.0_DP end if - vol = vol + sum(current(N)%age(:)) * recyratio + age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio) + !vol = vol + sum(current(N)%age(:)) * recyratio linmelt = current(N)%meltvolume * ratio meltinejecta = meltinejecta + linmelt distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) @@ -177,9 +177,12 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vseg = regolith_streamtube_volume_func(eradi,max(xmints,ri),rip1,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio ratio = max((vseg-vsh),0.0_DP) / current(N)%totvolume - vol = vol + sum(current(N)%age(:)) * recyratio + if (ratio > 1) then + ratio = 1.0_DP + end if + age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio) + !vol = vol + sum(current(N)%age(:)) * recyratio linmelt = current(N)%meltvolume * ratio meltinejecta = meltinejecta + linmelt distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) @@ -194,6 +197,10 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vmare = vmare + (vsgly - totseb) * current(N)%comp totseb = vsgly ratio = vsgly / current(N)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if + age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio) linmelt = current(N)%meltvolume * ratio meltinejecta = meltinejecta + linmelt distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index ca271fac..7b015b7d 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -122,8 +122,8 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer meltinejecta = surfi%regolayer(M)%meltvolume * ratio distvol(:) = distvol(:) + (surfi%regolayer(M)%distvol(:) * ratio) totvol = vseg - age_collector(:) = age_collector(:) + surfi%regolayer(M)%age(:) * recyratio - vol = vol + sum(age_collector(:)) + age_collector(:) = age_collector(:) + (surfi%regolayer(M)%age(:) * ratio) + !vol = vol + sum(age_collector(:)) ! write(*,*) '1',eradi, xmints, xsfints, & ! vseg/user%pix**2, (vseg-vsh)/user%pix**2, recyratio end if @@ -195,8 +195,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer if (rrighti > xmints) then vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rrighti,rrightf) recyratio = max((vsgly2 -vsh),0.0_DP)/ (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio - vol = vol + sum(current(N)%age(:)) * recyratio + !vol = vol + sum(current(N)%age(:)) * recyratio ratio2 = max((vsgly2-vsh),0.0_DP) / current(N)%totvolume if (ratio2 > 1) then ratio2 = 1.0_DP @@ -209,8 +208,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer if (rlefti > xmints) then vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rlefti,rleftf) recyratio = max((vsgly1-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio - vol = vol + sum(current(N)%age(:)) * recyratio + !vol = vol + sum(current(N)%age(:)) * recyratio ratio = max((vsgly1-vsh),0.0_DP) / current(N)%totvolume if (ratio > 1) then ratio = 1.0_DP @@ -224,8 +222,10 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer !distvol(:) = distvol(:) + (current(N)%meltdist(:)*(max((vsgly1-vsh),0.0_DP)+max((vsgly2-vsh2),0.0_DP))) if (ratio + ratio2 < 1.0_DP) then distvol(:) = distvol(:) + (current(N)%distvol(:)*(ratio)) + (current(N)%distvol(:)*(ratio2)) + age_collector(:) = age_collector(:) + (current(N)%age(:)*(ratio)) + (current(N)%age*(ratio2)) else distvol(:) = distvol(:) + current(N)%distvol(:) + age_collector(:) = age_collector(:) + current(N)%age(:) end if totvol = totvol + vsgly1 + vsgly2 !N = N - 1 diff --git a/src/regolith/regolith_traverse_streamtube.f90 b/src/regolith/regolith_traverse_streamtube.f90 index cdceb053..f777aabc 100644 --- a/src/regolith/regolith_traverse_streamtube.f90 +++ b/src/regolith/regolith_traverse_streamtube.f90 @@ -90,11 +90,11 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne vseg = regolith_streamtube_volume_func(eradi,max(xmints,ri),rip1,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) recyratio = max(vseg-vsh,0.0_DP) / (user%pix**2) / surfi%regolayer(N)%thickness - age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio ratio = max(vseg-vsh,0.0_DP) / surfi%regolayer(N)%totvolume if (ratio > 1) then ratio = 1.0_DP end if + age_collector(:) = age_collector(:) + (surfi%regolayer(N)%age(:) * ratio) meltinejecta = meltinejecta + (surfi%regolayer(N)%meltvolume * ratio) distvol(:) = distvol(:) + (surfi%regolayer(N)%distvol(:) * ratio) totvol = totvol + vseg diff --git a/src/util/util_traverse_pop_array.f90 b/src/util/util_traverse_pop_array.f90 index e6e650f4..243376c2 100644 --- a/src/util/util_traverse_pop_array.f90 +++ b/src/util/util_traverse_pop_array.f90 @@ -86,30 +86,8 @@ subroutine util_traverse_pop_array(user,regolayer,traverse_depth,poppedarray) ! copy regolayer from 1 to maxi to temp variable, then deallocate regolayer, then movealloc templayer onto regolayer <--may need temp array allocate(oldregodata,source=regolayer(1:maxi)) deallocate(regolayer) - call move_alloc(oldregodata,regolayer) ! right intents? + call move_alloc(oldregodata,regolayer) - ! if (z <= depth) then - ! dz = depth - z - - !*****the following lines may still be needed, especially if they deal with thickness:***** - - ! oldregodata = regolayer(maxi) - ! oldregodata%thickness = z - ! oldregodata%age(:) = z / regolayer(maxi)%thickness * regolayer(maxi)%age(:) - ! recyratio = dz / regolayer(maxi)%thickness - ! regolayer(maxi)%age(:) = recyratio * regolayer(maxi)%age(:) - ! regolayer(maxi)%thickness = dz - - !******************** - - - !call util_push_array(poppedarray,oldregodata) <--not needed; just editing in place - ! else - ! z = z - regolayer%thickness - ! call util_pop_array(regolayer,oldregodata) - ! call util_push_array(poppedarray,oldregodata) - ! depth = regolayer%thickness - ! end if return end subroutine util_traverse_pop_array