Moved the regotracking part of crater_form_interior to a new subroutine, regolith_interior, called after ejecta_emplace

src/Makefile.am

@@ -126,6 +126,7 @@ regolith/regolith_melt_zone_superdomain.f90\
 regolith/regolith_streamtube_volume_func.f90\
 regolith/regolith_shock_damage_zone.f90\
 regolith/regolith_shock_damage.f90\
+regolith/regolith_interior.f90\
 porosity/porosity_form_interior.f90\
 realistic/realistic_perlin_noise.f90\
 main/CTEM.f90

src/crater/crater_emplace.f90

@@ -51,7 +51,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
+subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval)
    use module_globals
    use module_util
    use module_porosity   
@@ -64,6 +64,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
    type(cratertype),intent(inout) :: crater
    type(domaintype),intent(inout) :: domain
    real(DP),intent(out) :: deltaMtot
+   integer(I4B),intent(out) :: incval
 
    ! Internal variables
    real(DP) :: lradsq,newelev, x_relative, y_relative 
@@ -81,6 +82,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
    fradsq = crater%frad**2
    deltaMtot = 0.0_DP !ejbmass
    incsq = inc**2
+   incval = inc
 
    ! This loop may not be parallelizable because of the linked list operation inside crater_form_interior
    do j=-inc,inc  ! Do the loop in pixel space

src/crater/crater_form_interior.f90

@@ -75,17 +75,11 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
    elchange  = newdem - surfi%dem
    deltaMi = elchange
    surfi%dem = newdem
+   surfi%abselc = abs(elchange)
 
    !change ejecta coverage
    surfi%ejcov = max(surfi%ejcov + elchange,0.0_DP)
 
-   if (user%doregotrack) then
-      call util_traverse_pop_array(surfi%regolayer,abs(elchange),poppedarray)
-      !call util_destroy_list(poppedlist)
-      deallocate(poppedarray)
-   end if
-
-
    return
 end subroutine crater_form_interior
 

src/crater/crater_populate.f90

@@ -78,7 +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
+   integer(I4B)            :: nmixingtimes, incval
 
    ! ejecta blanket array
    type(ejbtype),dimension(EJBTABSIZE) :: ejb       ! Ejecta blanket lookup table
@@ -249,7 +249,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
          call crater_averages(user,surf,crater)
 
          ! Place crater onto the surface
-         call crater_emplace(user,surf,crater,domain,ejbmass)
+         call crater_emplace(user,surf,crater,domain,ejbmass,incval)
          if (abs(ejbmass) < 2*tiny(1.0_DP)) cycle
 
          call ejecta_distance_estimate(user,crater,domain,crater%ejdis) ! Fast but imprecise estimate of the total ejecta distance
@@ -272,6 +272,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
             ejtble = 0
          end if
 
+         if (user%doregotrack) call regolith_interior(user,surf,crater,incval)
+
          if (user%dorealistic) call crater_realistic_topography(user,surf,crater,domain,ejecta_dem) 
          deallocate(ejecta_dem)
 

src/crater/module_crater.f90

@@ -97,14 +97,15 @@ end subroutine crater_averages
    end interface
 
    interface
-      subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
+      subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(inout) :: surf
       type(cratertype),intent(inout) :: crater
       type(domaintype),intent(inout) :: domain
       real(DP),intent(out) :: deltaMtot
+      integer(I4B),intent(out) :: incval
       end subroutine crater_emplace
    end interface
 

src/globals/module_globals.f90

@@ -92,6 +92,7 @@ module module_globals
    real(DP) :: ejcov                ! Ejecta coverage
    real(DP) :: dem                  ! Digital elevation model
    real(DP) :: mantle               ! Height of mantle (should be smaller than dem)
+   real(DP) :: abselc               ! abs(elchange)
    ! type(regolisttype), pointer :: regolayer => null() ! Pointer to the top of the regolith layer stack
    ! type(regolisttype), pointer :: porolayer => null() ! Pointer to the top of the porosity layer stack
    !type(regolisttype), pointer :: regolayer => null() ! Pointer to the top of the regolith layer stack

src/init/init_surf.f90

@@ -31,6 +31,7 @@ subroutine init_surf(user,surf,domain)
 
    surf%ejcov  = 0.0_DP
    surf%dem    = 0.0_DP
+   surf%abselc = 0.0_DP
    do layer = 1,user%numlayers
       surf%diam(layer)   = 0.0_DP
       surf%xl(layer)     = 0.0_SP

src/regolith/module_regolith.f90

@@ -320,5 +320,15 @@ function regolith_shock_damage(erad,deltar,xmints,xsfints,xleft,xright) result(v
      real(DP)                    :: vsh
      end function regolith_shock_damage
    end interface
+
+   interface
+      subroutine regolith_interior(user,surf,crater,incval)
+      use module_globals
+      type(usertype),intent(in) :: user
+      type(surftype),dimension(:,:),intent(inout) :: surf
+      type(cratertype),intent(in) :: crater
+      integer(I4B),intent(in)     :: incval
+      end subroutine regolith_interior
+   end interface
 
 end module

src/regolith/regolith_interior.f90

@@ -0,0 +1,68 @@
+!**********************************************************************************************************************************
+!
+!  Unit Name   : regolith_interior
+!  Unit Type   : subroutine
+!  Project     : CTEM
+!  Language    : Fortran 2003
+!
+!  Description : Pops off and destroys layers in the crater inteorior and fills with melt sheet
+!  
+!
+!  Input
+!    Arguments : user,surf,crater,inc
+!
+!  Output
+!    Arguments : surf
+!           
+! 
+!  Notes       :
+!
+!**********************************************************************************************************************************
+subroutine regolith_interior(user,surf,crater,incval)
+    use module_globals 
+    use module_util
+    use module_regolith, EXCEPT_THIS_ONE => regolith_interior
+    implicit none
+
+    !Arguments
+    type(usertype),intent(in) :: user
+    type(surftype),dimension(:,:),intent(inout) :: surf
+    type(cratertype),intent(in) :: crater
+    integer(I4B),intent(in)     :: incval
+
+    !internal variables
+    integer(I4B) xpi,ypi,i,j,inc,incsq,iradsq
+    real(DP) :: lradsq, x_relative, y_relative, xp, yp 
+    type(regodatatype),dimension(:),allocatable :: poppedarray
+
+    !Executable code
+
+    inc = incval
+
+    do j=-inc,inc
+        do i=-inc,inc
+            iradsq = i**2 + j**2
+
+            xpi = crater%xlpx + i
+            ypi = crater%ylpx + j
+
+            xp = xpi * user%pix
+            yp = ypi * user%pix
+
+            call util_periodic(xpi,ypi,user%gridsize)
+            x_relative = (crater%xl - xp)
+            y_relative = (crater%yl - yp)
+            lradsq = x_relative**2 + y_relative**2
+
+            if (lradsq > crater%frad**2) cycle
+            call util_traverse_pop_array(surf(xpi,ypi)%regolayer,surf(xpi,ypi)%abselc,poppedarray)
+            deallocate(poppedarray)
+
+            !Adding the melt sheet goes here! :)
+        end do
+    end do
+
+    return
+end subroutine regolith_interior
+
+