From 985f49b78a7e481149c9e33eed8a002d1b7bdf98 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Fri, 17 Feb 2023 14:49:21 -0500 Subject: [PATCH] Melt distributions are now tracked through the streamtube --- src/regolith/module_regolith.f90 | 12 ++++++++---- src/regolith/regolith_streamtube.f90 | 17 ++++++++++++----- src/regolith/regolith_streamtube_head.f90 | 6 +++++- src/regolith/regolith_streamtube_lineseg.f90 | 9 ++++++++- src/regolith/regolith_subpixel_streamtube.f90 | 11 ++++++++--- src/regolith/regolith_traverse_streamtube.f90 | 6 ++++-- 6 files changed, 45 insertions(+), 16 deletions(-) diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90 index fbe18f30..2f0e140d 100644 --- a/src/regolith/module_regolith.f90 +++ b/src/regolith/module_regolith.f90 @@ -95,7 +95,7 @@ end subroutine regolith_streamtube interface subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,depthb,meltinejecta,totvol) + age_collector,xmints,xsfints,depthb,meltinejecta,totvol,distvol) use module_globals implicit none type(usertype),intent(in) :: user @@ -107,12 +107,13 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints real(DP),intent(in) :: xsfints, depthb + real(SP),dimension(:),intent(inout) :: distvol end subroutine regolith_traverse_streamtube end interface interface subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,vol,meltinejecta,totvol) + age_collector,xmints,xsfints,vol,meltinejecta,totvol,distvol) use module_globals implicit none type(usertype),intent(in) :: user @@ -125,12 +126,13 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer real(DP),intent(in) :: xmints real(DP),intent(in) :: xsfints real(DP),intent(inout) :: vol + real(SP),dimension(:),intent(inout) :: distvol end subroutine regolith_subpixel_streamtube end interface interface subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad,eradi,deltar,newlayer,vmare,& - totseb,age_collector,xmints,xsfints,depthb,meltinejecta,totvol) + totseb,age_collector,xmints,xsfints,depthb,meltinejecta,totvol,distvol) use module_globals implicit none type(usertype),intent(in) :: user @@ -142,11 +144,12 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints real(DP),intent(in) :: xsfints, depthb + real(SP),dimension(:),intent(inout) :: distvol end subroutine regolith_streamtube_lineseg end interface interface - subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,meltinejecta,totvol) + subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,meltinejecta,totvol,distvol) use module_globals implicit none type(usertype),intent(in) :: user @@ -155,6 +158,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector real(DP),intent(in) :: deltar real(DP),intent(inout) :: totmare,tots real(SP),dimension(:),intent(inout) :: age_collector + real(SP),dimension(:),intent(inout) :: distvol end subroutine regolith_streamtube_head end interface diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index 9f5aa1b5..48ecb89c 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -91,6 +91,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, real(DP) :: vst,vbody,rbody,vmare,totmare,totseb,tots real(DP) :: meltinejecta, totvol type(regodatatype) :: newlayer + real(SP),dimension(:),allocatable :: distvol ! Constrain the tangital tube's volume with CTEM result real(DP) :: k1,k2,k3,k4,c1,c2 @@ -238,6 +239,8 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, depthb = crater%imp / 2.0_DP meltinejecta = 0.0_DP totvol = 0.0_DP + allocate(distvol(domain%rcnum)) + distvol(:) = 0.0_SP ! if (eradc<=user%testimp) then ! write(*,*) lrad/crater%frad, user%testimp, crater%frad, rm, deltar, eradc, eradi, erado, ebh, newlayer%meltfrac @@ -258,7 +261,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%thickness = vseg/(user%pix**2) call util_periodic(xstpi,ystpi,user%gridsize) call regolith_subpixel_streamtube(user,surf(xstpi,ystpi),deltar,ri,rip1,eradi,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,vol,meltinejecta,totvol) + age_collector,xmints,xsfints,vol,meltinejecta,totvol,distvol) newlayer%age(:) = newlayer%age(:) + age_collector(:) @@ -269,6 +272,8 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%age(:) = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP) newlayer%meltvolume = newlayer%meltvolume + meltinejecta newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume + distvol(:) = distvol(:) + (newlayer%ejm*newlayer%meltdist(:)) + newlayer%meltdist(:) = distvol(:) / newlayer%totvolume if (newlayer%meltfrac > 1.0_DP) then write(*,*) "Melt fraction >1! (Subpixel)", xpi,ypi,crater%timestamp,crater%fcrat,crater%xlpx,crater%ylpx,& newlayer%meltvolume, newlayer%totvolume, newlayer%ejm, newlayer%ejmf, totvol @@ -285,7 +290,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%thickness = vseg/(user%pix**2) call util_periodic(xstpi,ystpi,user%gridsize) call regolith_traverse_streamtube(user,surf(xstpi,ystpi),deltar,rbody,eradi,eradi,erado,newlayer,vmare,& - totseb,age_collector,xmints,xsfints,depthb,meltinejecta,totvol) + totseb,age_collector,xmints,xsfints,depthb,meltinejecta,totvol,distvol) totmare = totmare + vmare tots = tots + totseb end if @@ -302,7 +307,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%thickness = vseg/(user%pix**2) call util_periodic(xstpi,ystpi,user%gridsize) call regolith_traverse_streamtube(user,surf(xstpi,ystpi),deltar,ri,rip1,eradi,erado,newlayer,vmare,& - totseb,age_collector,xmints,xsfints,depthb,meltinejecta,totvol) + totseb,age_collector,xmints,xsfints,depthb,meltinejecta,totvol,distvol) totmare = totmare + vmare tots = tots + totseb end if @@ -311,7 +316,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, xstpi = crater%xlpx + nint(eradc*xl/lrad/user%pix) ystpi = crater%ylpx + nint(eradc*yl/lrad/user%pix) call util_periodic(xstpi,ystpi,user%gridsize) - call regolith_streamtube_head(user,surf(xstpi,ystpi),deltar,totmare,tots,age_collector,meltinejecta,totvol) + call regolith_streamtube_head(user,surf(xstpi,ystpi),deltar,totmare,tots,age_collector,meltinejecta,totvol,distvol) newlayer%thickness = ebh newlayer%comp = min(totmare/tots, 1.0_DP) @@ -320,6 +325,8 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, if (newlayer%ejmf < 1.0_DP) then newlayer%meltvolume = newlayer%meltvolume + meltinejecta newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume + distvol(:) = distvol(:) + (newlayer%ejm*newlayer%meltdist(:)) + newlayer%meltdist(:) = distvol(:) / newlayer%totvolume end if if (newlayer%meltfrac > 1.0_DP) then write(*,*) "Melt fraction >1! (Traverse)", xpi,ypi,crater%timestamp,crater%fcrat,crater%xlpx,crater%ylpx,& @@ -330,7 +337,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, call util_push_array(surf(xpi,ypi)%regolayer,newlayer) - deallocate(xints,yints) + deallocate(xints,yints,distvol) return end subroutine regolith_streamtube diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index 6d16bd00..ac95a08c 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -18,7 +18,7 @@ ! Notes : The stream tube's head is always attached to the surface. ! !********************************************************************************************************************************** -subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,meltinejecta,totvol) +subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,meltinejecta,totvol,distvol) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_streamtube_head implicit none @@ -29,6 +29,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector real(DP),intent(in) :: deltar real(DP),intent(inout) :: totmare,tots real(SP),dimension(:),intent(inout) :: age_collector + real(SP),dimension(:),intent(inout) :: distvol ! internal variables !type(regolisttype),pointer :: current @@ -67,6 +68,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio headmeltvol = current(N)%meltfrac * recyratio * vsgly meltinejecta = meltinejecta + headmeltvol + distvol(:) = distvol + (current(N)%meltdist(:)*recyratio*vsgly) totvol = totvol + vsgly else ! head is not intersected with layers. @@ -81,6 +83,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio headmeltvol = current(N)%meltfrac * recyratio * tothead meltinejecta = meltinejecta + headmeltvol + distvol(:) = distvol + (current(N)%meltdist(:)*recyratio*tothead) totvol = totvol + tothead !current => current%next N = N - 1 @@ -94,6 +97,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio headmeltvol = current(N)%meltfrac * recyratio * tothead meltinejecta = meltinejecta + headmeltvol + distvol(:) = distvol + (current(N)%meltdist(:)*recyratio*tothead) totvol = totvol + tothead exit end if diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90 index d2dbf577..77646e4f 100644 --- a/src/regolith/regolith_streamtube_lineseg.f90 +++ b/src/regolith/regolith_streamtube_lineseg.f90 @@ -19,7 +19,7 @@ ! !********************************************************************************************************************************** subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad,eradi,deltar,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,depthb,meltinejecta,totvol) + age_collector,xmints,xsfints,depthb,meltinejecta,totvol,distvol) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_streamtube_lineseg implicit none @@ -33,6 +33,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints real(DP),intent(in) :: xsfints, depthb + real(SP),dimension(:),intent(inout) :: distvol ! internal variables !type(regolisttype),pointer :: current @@ -85,6 +86,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vol = vol + sum(current(N-1)%age(:)) * recyratio linmelt = current(N-1)%meltfrac * vsgly * recyratio meltinejecta = meltinejecta + linmelt + distvol(:) = distvol(:) + (current(N-1)%meltdist(:)*vsgly*recyratio) totvol = totvol + vsgly else if (ri <= xmints .and. rip1 > xmints) then vseg = regolith_streamtube_volume_func(eradi,xmints,rip1,deltar) @@ -94,6 +96,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vol = vol + sum(current(N-1)%age(:)) * recyratio linmelt = current(N-1)%meltfrac * vseg * recyratio meltinejecta = meltinejecta + linmelt + distvol(:) = distvol(:) + (current(N-1)%meltdist(:)*vseg*recyratio) totvol = totvol + vseg end if exit @@ -123,6 +126,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vol = vol + sum(current(N)%age(:)) * recyratio linmelt = current(N)%meltfrac * vseg * recyratio meltinejecta = meltinejecta + linmelt + distvol(:) = distvol(:) + (current(N)%meltdist(:)*vseg*recyratio) totvol = totvol + vseg end if @@ -135,6 +139,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vol = vol + sum(current(N)%age(:)) * recyratio linmelt = current(N)%meltfrac * vseg * recyratio meltinejecta = meltinejecta + linmelt + distvol(:) = distvol(:) + (current(N)%meltdist(:)*vseg*recyratio) totvol = totvol + vseg end if !current => current%next @@ -159,6 +164,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vol = vol + sum(current(N)%age(:)) * recyratio linmelt = current(N)%meltfrac * vseg * recyratio meltinejecta = meltinejecta + linmelt + distvol(:) = distvol(:) + (current(N)%meltdist(:)*vseg*recyratio) totvol = totvol + vseg end if @@ -171,6 +177,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad totseb = vsgly linmelt = current(N)%meltfrac * vsgly meltinejecta = meltinejecta + linmelt + distvol(:) = distvol(:) + (current(N)%meltdist(:)*vsgly) totvol = totvol + vsgly exit end if diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index 48770072..46758758 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -51,7 +51,7 @@ ! !********************************************************************************************************************************** subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,vol,meltinejecta,totvol) + age_collector,xmints,xsfints,vol,meltinejecta,totvol,distvol) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_subpixel_streamtube implicit none @@ -67,6 +67,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer real(DP),intent(in) :: xmints real(DP),intent(in) :: xsfints real(DP),intent(inout) :: vol + real(SP),dimension(:),intent(inout) :: distvol ! Internal variables @@ -79,7 +80,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer integer(I4B) :: N ! Stream tube's distance from the edge of a melt zone - real(DP) :: zm, recyratio, xmints1, vseg + real(DP) :: zm, recyratio, xmints1, vseg, recyratio2 ! Parameters for calculating shocked segment of stream tube real(DP) :: x_up_sh, x_low_sh, vsh @@ -102,6 +103,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer totvol = 0.0_DP mvl = 0.0_DP mvr = 0.0_DP + recyratio2 = 0.0_DP ! Two cases: subpixel is inside the first layer, and its volume is simply the landing ejecta blanket. if (zend>=zmax) then @@ -112,6 +114,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,0.0_DP,eradi) recyratio = max(vseg-vsh,0.0 )/ (user%pix**2) / (surfi%regolayer(N)%thickness) meltinejecta = surfi%regolayer(N)%meltfrac * vseg * recyratio + distvol = distvol + (surfi%regolayer(N)%meltdist(:)*vseg*recyratio) totvol = vseg age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio vol = vol + sum(age_collector(:)) @@ -189,6 +192,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio mvr = vsgly2 * current(N)%meltfrac * recyratio + recyratio2 = recyratio end if if (rlefti > xmints) then @@ -201,6 +205,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer !current => current%next meltinejecta = meltinejecta + mvl + mvr + distvol = distvol + (current(N)%meltdist(:)*((vsgly1*recyratio)+(vsgly2*recyratio2))) totvol = totvol + vsgly1 + vsgly2 N = N - 1 z = z + current(N)%thickness @@ -218,7 +223,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer end do end if - call regolith_streamtube_head(user,surfi,deltar,vmare,totseb,age_collector,meltinejecta,totvol) + call regolith_streamtube_head(user,surfi,deltar,vmare,totseb,age_collector,meltinejecta,totvol,distvol) return end subroutine regolith_subpixel_streamtube diff --git a/src/regolith/regolith_traverse_streamtube.f90 b/src/regolith/regolith_traverse_streamtube.f90 index b4f6e7af..7b56a4bd 100644 --- a/src/regolith/regolith_traverse_streamtube.f90 +++ b/src/regolith/regolith_traverse_streamtube.f90 @@ -19,7 +19,7 @@ ! !********************************************************************************************************************************** subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,depthb,meltinejecta,totvol) + age_collector,xmints,xsfints,depthb,meltinejecta,totvol,distvol) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_traverse_streamtube implicit none @@ -34,6 +34,7 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints real(DP),intent(in) :: xsfints, depthb + real(SP),dimension(:),intent(inout) :: distvol ! Internal variables @@ -89,13 +90,14 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne recyratio = max(vseg-vsh,0.0_DP) / (user%pix**2) / surfi%regolayer(N)%thickness age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio meltinejecta = meltinejecta + surfi%regolayer(N)%meltfrac * vseg * recyratio + distvol(:) = distvol(:) + (surfi%regolayer(N)%meltdist(:)*vseg*recyratio) totvol = totvol + vseg end if else call regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad,eradi,deltar,& - newlayer,vmare,totseb,age_collector,xmints,xsfints,depthb,meltinejecta,totvol) + newlayer,vmare,totseb,age_collector,xmints,xsfints,depthb,meltinejecta,totvol,distvol) end if