From 7ca58d75ba0ea9238d12c6f8de31dc746704b4a9 Mon Sep 17 00:00:00 2001
From: Austin Michael Blevins <blevins2@bell-a065.rcac.purdue.edu>
Date: Mon, 23 Jan 2023 11:33:38 -0500
Subject: [PATCH] made it so it runs subpixel after every crater

---
 src/crater/crater_populate.f90        | 9 ++++++++-
 src/globals/module_globals.f90        | 2 +-
 src/regolith/regolith_depth_model.f90 | 4 +++-
 3 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90
index 5e230b9c..ea40c8a0 100644
--- a/src/crater/crater_populate.f90
+++ b/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
diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90
index 50fc657e..34147a7b 100644
--- a/src/globals/module_globals.f90
+++ b/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
diff --git a/src/regolith/regolith_depth_model.f90 b/src/regolith/regolith_depth_model.f90
index 8f37eff8..6b424055 100644
--- a/src/regolith/regolith_depth_model.f90
+++ b/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