From a8df2d43f032a9d5928cef720018a3813ffb6ff1 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Thu, 16 Feb 2023 16:13:42 -0500 Subject: [PATCH] Moved the regotracking part of crater_form_interior to a new subroutine, regolith_interior, called after ejecta_emplace --- src/Makefile.am | 1 + src/crater/crater_emplace.f90 | 4 +- src/crater/crater_form_interior.f90 | 8 +--- src/crater/crater_populate.f90 | 6 ++- src/crater/module_crater.f90 | 3 +- src/globals/module_globals.f90 | 1 + src/init/init_surf.f90 | 1 + src/regolith/module_regolith.f90 | 10 +++++ src/regolith/regolith_interior.f90 | 68 +++++++++++++++++++++++++++++ 9 files changed, 91 insertions(+), 11 deletions(-) create mode 100644 src/regolith/regolith_interior.f90 diff --git a/src/Makefile.am b/src/Makefile.am index 52a0bbbc..2ea4799b 100755 --- a/src/Makefile.am +++ b/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 diff --git a/src/crater/crater_emplace.f90 b/src/crater/crater_emplace.f90 index 885143cd..bdd55c43 100644 --- a/src/crater/crater_emplace.f90 +++ b/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 diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90 index aacafde6..0b2b5334 100644 --- a/src/crater/crater_form_interior.f90 +++ b/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 diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 6ae9faf1..141daa66 100644 --- a/src/crater/crater_populate.f90 +++ b/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) diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index b946c25b..f6fa0bb3 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -97,7 +97,7 @@ 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 @@ -105,6 +105,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 end subroutine crater_emplace end interface diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 3260127a..cc856a18 100644 --- a/src/globals/module_globals.f90 +++ b/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 diff --git a/src/init/init_surf.f90 b/src/init/init_surf.f90 index 955ef97e..bde018f6 100644 --- a/src/init/init_surf.f90 +++ b/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 diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90 index f7e01a87..fbe18f30 100644 --- a/src/regolith/module_regolith.f90 +++ b/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 diff --git a/src/regolith/regolith_interior.f90 b/src/regolith/regolith_interior.f90 new file mode 100644 index 00000000..a1b01ba0 --- /dev/null +++ b/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 + +