Overhauled melt tracking and made it much simpler

MintonGroup · Feb 3, 2023 · 3a3b282 · 3a3b282
1 parent 2291272
commit 3a3b282
Show file tree

Hide file tree

Showing 7 changed files with 54 additions and 102 deletions.
diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90
@@ -387,7 +387,5 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
 
    end if
 
-   write(*,*) 'Total times mixing was called:',nmixingtimes
-
    return
 end subroutine crater_populate
diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90
@@ -35,7 +35,7 @@ module module_regolith
 save
 
    interface
-      subroutine regolith_melt_zone(user,crater,dimp,vimp,rmelt,depthb, volm)
+      subroutine regolith_melt_zone(user,crater,dimp,vimp,rmelt,depthb,volm)
       use module_globals
       type(usertype),intent(in) :: user
       type(cratertype),intent(inout) :: crater
@@ -77,7 +77,7 @@ end subroutine regolith_transport
 
    interface 
       subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,&
-                 rm,vsq,age,age_resolution, volm)
+                 rm,vsq,age,age_resolution,volm)
       use module_globals 
       implicit none
       type(usertype),intent(in) :: user
@@ -95,13 +95,14 @@ end subroutine regolith_streamtube
 
    interface 
       subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,newlayer,vmare,totseb,&
-                 age_collector,xmints,xsfints,rsh,depthb,mixedregodata)
+                 age_collector,xmints,xsfints,rsh,depthb,meltinejecta)
       use module_globals 
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),intent(inout) :: surfi
       real(DP),intent(in)            :: deltar,ri,rip1,eradi,erado
-      type(regodatatype),intent(inout) :: newlayer, mixedregodata
+      type(regodatatype),intent(inout) :: newlayer
+      real(DP),intent(out)            :: meltinejecta
       real(DP),intent(out) :: vmare,totseb
       real(SP),dimension(:),intent(inout) :: age_collector
       real(DP),intent(in)             :: xmints
@@ -111,13 +112,14 @@ end subroutine regolith_traverse_streamtube
 
    interface 
       subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer,vmare,totseb,& 
-                 age_collector,xmints,xsfints,vol,mixedregodata) 
+                 age_collector,xmints,xsfints,vol,meltinejecta) 
       use module_globals 
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),intent(in) :: surfi
       real(DP),intent(in)       :: deltar,ri,rip1,eradi
-      type(regodatatype),intent(inout)    :: newlayer, mixedregodata
+      type(regodatatype),intent(inout)    :: newlayer
+      real(DP),intent(out)                :: meltinejecta
       real(DP),intent(out)                :: vmare,totseb
       real(SP),dimension(:),intent(inout) :: age_collector
       real(DP),intent(in)                 :: xmints
@@ -128,13 +130,14 @@ end subroutine regolith_subpixel_streamtube
 
    interface 
       subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad,eradi,deltar,newlayer,vmare,&
-      totseb,age_collector,xmints,xsfints,depthb,mixedregodata)
+      totseb,age_collector,xmints,xsfints,depthb,meltinejecta)
       use module_globals 
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),intent(in) :: surfi
       real(DP),intent(in) :: thetast,ri,rip1,zmin,zmax,erad,eradi,deltar
-      type(regodatatype),intent(inout) :: newlayer,mixedregodata
+      type(regodatatype),intent(inout) :: newlayer
+      real(DP),intent(inout)           :: meltinejecta
       real(DP),intent(inout) :: vmare,totseb
       real(SP),dimension(:),intent(inout) :: age_collector
       real(DP),intent(in)             :: xmints
@@ -143,12 +146,12 @@ end subroutine regolith_streamtube_lineseg
    end interface 
 
    interface 
-      subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,mixedregodata)
+      subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,meltinejecta)
       use module_globals 
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),intent(in) :: surfi
-      type(regodatatype),intent(inout) :: mixedregodata
+      real(DP),intent(inout)    :: meltinejecta
       real(DP),intent(in) :: deltar
       real(DP),intent(inout) :: totmare,tots
       real(SP),dimension(:),intent(inout) :: age_collector

diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90
@@ -89,7 +89,8 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
    integer(I4B) :: i,j,k,toti,totj,toty,cnt,xstpi,ystpi
    real(DP)     :: vtot,vseg,ri,rip1,xc,yc,thetast
    real(DP)     :: vst,vbody,rbody,vmare,totmare,totseb,tots
-   type(regodatatype) :: newlayer, mixedregodata
+   real(DP)     :: meltinejecta
+   type(regodatatype) :: newlayer
 
    ! Constrain the tangital tube's volume with CTEM result
    real(DP)     :: k1,k2,k3,k4,c1,c2
@@ -117,10 +118,6 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
 
    ! Executalbe code
 
-   mixedregodata%totvolume = 0
-   mixedregodata%meltvolume = 0
-   mixedregodata%meltfrac = 0
-
    ! ****** Interpolate radial distance, erad, for a given pixel *******
    ! outeredge = crater%frad + domain%ejbres * (EJBTABSIZE - 0.5_DP)
    ! inneredge = crater%frad + 0.5_DP * domain%ejbres
@@ -259,7 +256,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,mixedregodata)
+           age_collector,xmints,xsfints,vol,meltinejecta)
 
       newlayer%age(:) = newlayer%age(:) + age_collector(:)
 
@@ -268,8 +265,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
       newlayer%thickness = ebh
       newlayer%comp      = min(totmare/tots, 1.0_DP)
       newlayer%age(:)    = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP)
-      newlayer%meltvolume = newlayer%meltvolume + mixedregodata%meltvolume
-      newlayer%totvolume  = newlayer%totvolume + mixedregodata%totvolume
+      newlayer%meltvolume = newlayer%meltvolume + meltinejecta
       newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume
 
    else
@@ -283,7 +279,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,rsh,depthb,mixedregodata)
+            totseb,age_collector,xmints,xsfints,rsh,depthb,meltinejecta)
          totmare = totmare + vmare
          tots = tots + totseb
       end if           
@@ -300,7 +296,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,rsh,depthb,mixedregodata)
+               totseb,age_collector,xmints,xsfints,rsh,depthb,meltinejecta)
             totmare = totmare + vmare
             tots = tots + totseb 
          end if
@@ -309,15 +305,17 @@ 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,mixedregodata)
+      call regolith_streamtube_head(user,surf(xstpi,ystpi),deltar,totmare,tots,age_collector,meltinejecta)
 
       newlayer%thickness = ebh
       newlayer%comp      = min(totmare/tots, 1.0_DP)
       newlayer%age(:)    = newlayer%age(:) + age_collector(:)
       newlayer%age(:)    = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP)
-      newlayer%meltvolume = newlayer%meltvolume + mixedregodata%meltvolume
-      newlayer%totvolume  = newlayer%totvolume + mixedregodata%totvolume
+      newlayer%meltvolume = newlayer%meltvolume + meltinejecta
       newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume
+      if (newlayer%meltfrac > 1.0_DP) then
+         write(*,*) "Melt fraction >1!", xpi,ypi,crater%timestamp,crater%fcrat,crater%xlpx,crater%ylpx
+      end if
 
   end if
 

diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90
@@ -18,14 +18,14 @@
 !  Notes       :  The stream tube's head is always attached to the surface. 
 !
 !**********************************************************************************************************************************
-subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,mixedregodata)
+subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,meltinejecta)
    use module_globals 
    use module_regolith, EXCEPT_THIS_ONE => regolith_streamtube_head
    implicit none
    ! arguemnts
    type(usertype),intent(in) :: user
    type(surftype),intent(in) :: surfi
-   type(regodatatype),intent(inout) :: mixedregodata
+   real(DP),intent(inout) :: meltinejecta
    real(DP),intent(in) :: deltar
    real(DP),intent(inout) :: totmare,tots
    real(SP),dimension(:),intent(inout) :: age_collector
@@ -66,12 +66,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
       recyratio = vsgly / (user%pix**2) /current(N)%thickness
       age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio
       headmeltvol = current(N)%meltfrac * recyratio * vsgly
-      mixedregodata%totvolume = mixedregodata%totvolume + vsgly
-      mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol
-      mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-      if (mixedregodata%meltfrac > 1.0_DP) then
-         write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
-      end if
+      meltinejecta = meltinejecta + headmeltvol
    else ! head is not intersected with layers. 
 
       do
@@ -83,13 +78,8 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
             totmarehead = totmarehead + vhead * vratio * current(N)%comp
             recyratio = vhead * vratio / (user%pix**2) / current(N)%thickness
             age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio
-            mixedregodata%totvolume = mixedregodata%totvolume + tothead
             headmeltvol = current(N)%meltfrac * recyratio * tothead
-            mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol
-            mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-            if (mixedregodata%meltfrac > 1.0_DP) then
-              write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
-            end if
+            meltinejecta = meltinejecta + headmeltvol
             !current => current%next
             N = N - 1
             z = z + current(N)%thickness
@@ -100,13 +90,8 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
             tothead = vsgly
             recyratio = (vsgly - tothead) / (user%pix**2) / current(N)%thickness
             age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio
-            mixedregodata%totvolume = mixedregodata%totvolume + tothead
             headmeltvol = current(N)%meltfrac * recyratio * tothead
-            mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol
-            mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-            if (mixedregodata%meltfrac > 1.0_DP) then
-               write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
-            end if
+            meltinejecta = meltinejecta + headmeltvol
             exit
          end if
       end do

diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90
@@ -19,15 +19,16 @@
 !
 !**********************************************************************************************************************************
 subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad,eradi,deltar,newlayer,vmare,totseb,&
-           age_collector,xmints,xsfints,depthb,mixedregodata)
+           age_collector,xmints,xsfints,depthb,meltinejecta)
    use module_globals 
    use module_regolith, EXCEPT_THIS_ONE => regolith_streamtube_lineseg
    implicit none
    ! arguemnts
    type(usertype),intent(in) :: user
    type(surftype),intent(in) :: surfi
    real(DP),intent(in) :: thetast,ri,rip1,zmin,zmax,erad,eradi,deltar
-   type(regodatatype),intent(inout) :: newlayer,mixedregodata
+   type(regodatatype),intent(inout) :: newlayer
+   real(DP),intent(inout)           :: meltinejecta
    real(DP),intent(inout) :: vmare,totseb
    real(SP),dimension(:),intent(inout) :: age_collector
    real(DP),intent(in)             :: xmints
@@ -82,26 +83,16 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                 recyratio        = max(ebh_recyl,0.0_DP) / current(N-1)%thickness
                 age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio
                 vol              = vol + sum(current(N-1)%age(:)) * recyratio
-                mixedregodata%totvolume = mixedregodata%totvolume + (vsgly-vsh)
                 linmelt = current(N-1)%meltfrac * vsgly * recyratio
-                mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
-                mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                if (mixedregodata%meltfrac > 1.0_DP) then
-                  write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
-               end if
+                meltinejecta = meltinejecta + linmelt
              else if (ri <= xmints .and. rip1 > xmints) then 
                      vseg             = regolith_streamtube_volume_func(eradi,xmints,rip1,deltar)
                      vsh              = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) 
                      recyratio        = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N-1)%thickness 
                      age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio
                      vol              = vol + sum(current(N-1)%age(:)) * recyratio
-                     mixedregodata%totvolume = mixedregodata%totvolume + (vseg-vsh)
                      linmelt = current(N-1)%meltfrac * vseg * recyratio
-                     mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
-                     mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                     if (mixedregodata%meltfrac > 1.0_DP) then
-                        write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
-                     end if
+                     meltinejecta = meltinejecta + linmelt
              end if
              exit
           else
@@ -128,13 +119,8 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                   recyratio        = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness
                   age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio
                   vol              = vol + sum(current(N)%age(:)) * recyratio
-                  mixedregodata%totvolume = mixedregodata%totvolume + (vseg-vsh)
                   linmelt = current(N)%meltfrac * vseg * recyratio
-                  mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
-                  mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                  if (mixedregodata%meltfrac > 1.0_DP) then
-                     write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
-                  end if
+                  meltinejecta = meltinejecta + linmelt
                end if
 
                ! A segment coming from the side of emerging location of a streamtube rip1 
@@ -144,13 +130,8 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                   recyratio        = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness
                   age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio
                   vol              = vol + sum(current(N)%age(:)) * recyratio
-                  mixedregodata%totvolume = mixedregodata%totvolume + (vseg-vsh)
                   linmelt = current(N)%meltfrac * vseg * recyratio
-                  mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
-                  mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                  if (mixedregodata%meltfrac > 1.0_DP) then
-                     write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
-                  end if
+                  meltinejecta = meltinejecta + linmelt
                end if
                !current => current%next
                N = N - 1
@@ -172,13 +153,8 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                   recyratio        = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness
                   age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio
                   vol              = vol + sum(current(N)%age(:)) * recyratio
-                  mixedregodata%totvolume = mixedregodata%totvolume + (vseg-vsh)
                   linmelt = current(N)%meltfrac * vseg * recyratio
-                  mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
-                  mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                  if (mixedregodata%meltfrac > 1.0_DP) then
-                     write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
-                  end if
+                  meltinejecta = meltinejecta + linmelt
                end if
 
                exit
@@ -188,13 +164,8 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                vsgly = regolith_streamtube_volume_func(eradi,ri,rip1,deltar)
                vmare = vmare + (vsgly - totseb) * current(N)%comp
                totseb = vsgly
-               mixedregodata%totvolume = mixedregodata%totvolume + vsgly
                linmelt = current(N)%meltfrac * vsgly * recyratio
-               mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
-               mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-               if (mixedregodata%meltfrac > 1.0_DP) then
-                  write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
-               end if 
+               meltinejecta = meltinejecta + linmelt
                exit
        end if
-Original file line number
+Diff line change
@@ Expand Up @@
        end if
-       write(*,*) 'Total times mixing was called:',nmixingtimes
        return
     end subroutine crater_populate