made it so it runs subpixel after every crater

src/crater/crater_populate.f90

@@ -78,6 +78,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
    integer(I4B)            :: oldpbarpos
    real(DP),dimension(:,:),allocatable   :: ejecta_dem
    real(DP)                :: hmax, hmin
+   integer(I4B)            :: nmixingtimes
 
    ! ejecta blanket array
    type(ejbtype),dimension(EJBTABSIZE) :: ejb       ! Ejecta blanket lookup table
@@ -93,6 +94,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
    type(regolisttype),pointer                        :: current => null()
    real(DP)              :: age_resolution
 
+   nmixingtimes = 0
+
    if (user%testflag) then
       write(*,*) "Generating a test crater"
       write(*,*) "Dimp = ",user%testimp
@@ -318,7 +321,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
 
       ! Do periodic subpixel processes on the whole grid
       if (.not.user%testflag) then
-         if ((domain%subpixelcoverage / real(user%gridsize**2,kind=DP) > SUBPIXELCOVERAGE).or.(icrater == ntotcrat)) then
+         !if ((domain%subpixelcoverage / real(user%gridsize**2,kind=DP) > SUBPIXELCOVERAGE).or.(icrater == ntotcrat)) then
+         if (makecrater) then
             domain%subpixelcoverage = 0
             write(message,*) "Subpixel"
             call io_updatePbar(message)
@@ -333,6 +337,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
                call regolith_depth_model(user,domain,finterval,nflux,p)
                call regolith_subcrater_mix(user,surf,domain,nflux,finterval,p)
                age = age - finterval * user%interval
+               nmixingtimes = nmixingtimes + 1
             end if 
 
             icrater_last_subpixel = icrater
@@ -382,5 +387,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
 
    end if
 
+   write(*,*) 'Total times mixing was called:',nmixingtimes
+
    return
 end subroutine crater_populate

src/globals/module_globals.f90

@@ -52,7 +52,7 @@ module module_globals
 integer(I4B),parameter :: MAXLAYER=20          ! Maximum number of layers (you need roughly 1-2 layers per order of magnitude of 
                                                ! resolution
 real(DP),parameter :: TALLYCOVERAGE = 0.01_DP  ! The total area coverage to reach before a tally step is executed
-real(DP),parameter :: SUBPIXELCOVERAGE = 0.025_DP ! The total area coverage to reach before a subpixel evaluate step is executed: 0.05_DP
+real(DP),parameter :: SUBPIXELCOVERAGE = 0.00025_DP ! The total area coverage to reach before a subpixel evaluate step is executed: 0.05_DP
 real(DP),parameter :: COOKIESIZE = 3.0_DP      ! Relative size of old crater to new crater that cookie cutting is applied
                                                ! Only craters smaller than COOKIESIZE times the new crater are cookie cut
 integer(I2B),parameter :: MAXAGEBINS=60       ! Maximum number of bins in age distribution reset by impact melting

src/regolith/regolith_depth_model.f90

@@ -37,6 +37,7 @@ subroutine regolith_depth_model(user,domain,finterval,nflux,p)
    real(DP)     :: h, dr, f, fmin, fmax
    real(DP)     :: ntotsubcrat
    logical :: underflow
+   real(DP) :: psumfunc
 
    if (ieee_support_underflow_control(psum)) call ieee_set_underflow_mode(gradual=.true.)
    ! Smallest crater size in sub-pixel crater regime (regolith scaling column in "nflux")
@@ -96,7 +97,8 @@ subroutine regolith_depth_model(user,domain,finterval,nflux,p)
          f    = 0.5_DP * dr * ( fmin + fmax )
          psum = psum + f   
       end do
-      if (abs(psum) > epsilon(psum)) then
+      psumfunc = exp(-1.0_DP * PI * psum * t)
+      if (psumfunc > epsilon(psum)) then
          p(2,i) = 1.0_DP - exp(-1.0_DP * PI * psum * t)
       else
          p(2,i) = 1.0_DP