age tracking for local melt seems to work

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

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

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?

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

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

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

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)

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)

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

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