Melt distributions are now tracked through the streamtube

MintonGroup · Feb 17, 2023 · 985f49b · 985f49b
1 parent a8df2d4
commit 985f49b
Show file tree

Hide file tree

Showing 6 changed files with 45 additions and 16 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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