From a0be7090d1073466e0aabd746b31ab8acffcf08f Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Tue, 21 Feb 2023 11:44:46 -0500 Subject: [PATCH] added melt sheet emplacement. Compiles but not tested (want to make mixing optional for test craters first) --- src/crater/crater_emplace.f90 | 6 ++++-- src/crater/crater_populate.f90 | 9 +++++---- src/crater/module_crater.f90 | 4 ++-- src/ejecta/ejecta_emplace.f90 | 6 +++++- src/ejecta/module_ejecta.f90 | 3 ++- src/regolith/module_regolith.f90 | 6 ++++-- src/regolith/regolith_interior.f90 | 31 ++++++++++++++++++++++++------ 7 files changed, 47 insertions(+), 18 deletions(-) diff --git a/src/crater/crater_emplace.f90 b/src/crater/crater_emplace.f90 index bdd55c43..3a38043e 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,incval) +subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval,nmeltsheet) use module_globals use module_util use module_porosity @@ -64,7 +64,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval) type(cratertype),intent(inout) :: crater type(domaintype),intent(inout) :: domain real(DP),intent(out) :: deltaMtot - integer(I4B),intent(out) :: incval + integer(I4B),intent(out) :: incval,nmeltsheet ! Internal variables real(DP) :: lradsq,newelev, x_relative, y_relative @@ -83,6 +83,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval) deltaMtot = 0.0_DP !ejbmass incsq = inc**2 incval = inc + nmeltsheet = 0 ! 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 @@ -108,6 +109,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval) if (lradsq > crater%frad**2) cycle call crater_form_interior(user,surf(xpi,ypi),crater,x_relative,y_relative,newelev,deltaMi) deltaMtot = deltaMtot + deltaMi + nmeltsheet = nmeltsheet + 1 ! do porosity computation if (user%doporosity) ! It is still important to consider the physical meaning of frad and rad. diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 141daa66..ee354074 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -78,7 +78,8 @@ 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, incval + integer(I4B) :: nmixingtimes, incval, nmeltsheet + real(DP) :: vmeltsheet ! ejecta blanket array type(ejbtype),dimension(EJBTABSIZE) :: ejb ! Ejecta blanket lookup table @@ -249,7 +250,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,incval) + call crater_emplace(user,surf,crater,domain,ejbmass,incval,nmeltsheet) 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 @@ -267,12 +268,12 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt call ejecta_table_define(user,crater,domain,ejb,ejtble) !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) end if - call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,age,age_resolution,ejecta_dem) + call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,age,age_resolution,ejecta_dem,vmeltsheet) else ejtble = 0 end if - if (user%doregotrack) call regolith_interior(user,surf,crater,incval) + if (user%doregotrack) call regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltsheet) 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 f6fa0bb3..307d0b78 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,incval) + subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval,nmeltsheet) use module_globals implicit none type(usertype),intent(in) :: user @@ -105,7 +105,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot,incval) type(cratertype),intent(inout) :: crater type(domaintype),intent(inout) :: domain real(DP),intent(out) :: deltaMtot - integer(I4B),intent(out) :: incval + integer(I4B),intent(out) :: incval,nmeltsheet end subroutine crater_emplace end interface diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 2411112c..ce8436e7 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -75,7 +75,7 @@ ! The cutoff of ejecta thickness is still buggy. ! !********************************************************************************************************************************** -subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange) +subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange,vmeltsheet) use module_globals use module_util use module_io @@ -95,6 +95,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r real(DP),intent(in) :: age real(DP),intent(in) :: age_resolution real(DP),dimension(:,:),allocatable,intent(out) :: cumulative_elchange + real(DP),intent(out) :: vmeltsheet ! Internal variables real(DP) :: lrad,lradsq @@ -129,6 +130,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r if (user%doregotrack) call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm,totmelt) vmelt = 0.0_DP + vmeltsheet = 0.0_DP crater%vdepth = crater%rimheight + crater%floordepth crater%vrim = crater%rimheight @@ -290,6 +292,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r write(*,*) 'Total Melt: ', totmelt write(*,*) 'ejected / total melt:', vmelt/totmelt end if + vmeltsheet = totmelt - vmelt + ejbmass = sum(cumulative_elchange) ! Create buffer to prevent infinite hole bug diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index 2fb6bda5..dbf416b6 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/src/ejecta/module_ejecta.f90 @@ -25,7 +25,7 @@ module module_ejecta save interface - subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange) + subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange,vmeltsheet) use module_globals implicit none type(usertype),intent(in) :: user @@ -38,6 +38,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r real(DP),intent(in) :: age real(DP),intent(in) :: age_resolution real(DP),dimension(:,:),allocatable,intent(inout) :: cumulative_elchange + real(DP),intent(out) :: vmeltsheet end subroutine ejecta_emplace end interface diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90 index 2f0e140d..befe69a6 100644 --- a/src/regolith/module_regolith.f90 +++ b/src/regolith/module_regolith.f90 @@ -326,12 +326,14 @@ end function regolith_shock_damage end interface interface - subroutine regolith_interior(user,surf,crater,incval) + subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltsheet) use module_globals type(usertype),intent(in) :: user type(surftype),dimension(:,:),intent(inout) :: surf + type(domaintype),intent(in) :: domain type(cratertype),intent(in) :: crater - integer(I4B),intent(in) :: incval + integer(I4B),intent(in) :: incval, nmeltsheet + real(DP),intent(in) :: vmeltsheet end subroutine regolith_interior end interface diff --git a/src/regolith/regolith_interior.f90 b/src/regolith/regolith_interior.f90 index a1b01ba0..5159086b 100644 --- a/src/regolith/regolith_interior.f90 +++ b/src/regolith/regolith_interior.f90 @@ -9,16 +9,16 @@ ! ! ! Input -! Arguments : user,surf,crater,inc +! Arguments : user,surf,domain,crater,inc,nmeltsheet,vmeltsheet ! ! Output ! Arguments : surf ! ! -! Notes : +! Notes : nmeltsheet is the number of pixels that will contain the melt sheet; vmeltsheet is the melt sheet volume ! !********************************************************************************************************************************** -subroutine regolith_interior(user,surf,crater,incval) +subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltsheet) use module_globals use module_util use module_regolith, EXCEPT_THIS_ONE => regolith_interior @@ -27,16 +27,22 @@ subroutine regolith_interior(user,surf,crater,incval) !Arguments type(usertype),intent(in) :: user type(surftype),dimension(:,:),intent(inout) :: surf + type(domaintype),intent(in) :: domain type(cratertype),intent(in) :: crater - integer(I4B),intent(in) :: incval + integer(I4B),intent(in) :: incval, nmeltsheet + real(DP),intent(in) :: vmeltsheet !internal variables integer(I4B) xpi,ypi,i,j,inc,incsq,iradsq - real(DP) :: lradsq, x_relative, y_relative, xp, yp + real(DP) :: lradsq, x_relative, y_relative, xp, yp, hmeltsheet type(regodatatype),dimension(:),allocatable :: poppedarray + type(regodatatype) :: newlayer !Executable code + hmeltsheet = vmeltsheet / (nmeltsheet*user%gridsize*user%gridsize) + allocate(newlayer%meltdist(domain%rcnum)) + inc = incval do j=-inc,inc @@ -58,10 +64,23 @@ subroutine regolith_interior(user,surf,crater,incval) call util_traverse_pop_array(surf(xpi,ypi)%regolayer,surf(xpi,ypi)%abselc,poppedarray) deallocate(poppedarray) - !Adding the melt sheet goes here! :) + !fill top layer with melt sheet of given thickness hmeltsheet + newlayer%meltfrac = 1.0_DP + newlayer%ejm = 0.0_DP + newlayer%ejmf = 0.0_DP + newlayer%thickness = hmeltsheet + newlayer%meltvolume = vmeltsheet / nmeltsheet + newlayer%totvolume = newlayer%meltvolume + newlayer%meltdist(:) = 0.0_SP + if(domain%currentqmc) then + newlayer%meltdist(domain%nqmc) = newlayer%meltfrac + end if + call util_push_array(surf(xpi,ypi)%regolayer,newlayer) end do end do + deallocate(newlayer%meltdist) + return end subroutine regolith_interior