compiles with melt tracking

src/regolith/regolith_streamtube.f90

@@ -117,6 +117,10 @@ 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
@@ -254,9 +258,6 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
       vseg = regolith_streamtube_volume_func(eradi,0.0_DP,eradi,deltar)
       newlayer%thickness = vseg/(user%pix**2)
       call util_periodic(xstpi,ystpi,user%gridsize)
-      mixedregodata%meltfrac = surf(xstpi,ystpi)%regolayer%meltfrac
-      mixedregodata%meltvolume = surf(xstpi,ystpi)%regolayer%meltvolume
-      mixedregodata%totvolume = surf(xstpi,ystpi)%regolayer%totvolume
       call regolith_subpixel_streamtube(user,surf(xstpi,ystpi),deltar,ri,rip1,eradi,newlayer,vmare,totseb,&
            age_collector,xmints,xsfints,vol,mixedregodata)
 
@@ -281,11 +282,8 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
          vseg = regolith_streamtube_volume_func(eradi,rbody,eradi,deltar)
          newlayer%thickness = vseg/(user%pix**2)
          call util_periodic(xstpi,ystpi,user%gridsize)
-         mixedregodata%meltfrac = surf(xstpi,ystpi)%regolayer%meltfrac
-         mixedregodata%meltvolume = surf(xstpi,ystpi)%regolayer%meltvolume
-         mixedregodata%totvolume = surf(xstpi,ystpi)%regolayer%totvolume
          call regolith_traverse_streamtube(user,surf(xstpi,ystpi),deltar,rbody,eradi,eradi,erado,newlayer,vmare,&
-            totseb,age_collector,xmints,xsfints,rsh,depthb)
+            totseb,age_collector,xmints,xsfints,rsh,depthb,mixedregodata)
          totmare = totmare + vmare
          tots = tots + totseb
       end if           
@@ -301,9 +299,6 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
             vseg = regolith_streamtube_volume_func(eradi,ri,rip1,deltar)
             newlayer%thickness = vseg/(user%pix**2)
             call util_periodic(xstpi,ystpi,user%gridsize)
-            mixedregodata%meltfrac = surf(xstpi,ystpi)%regolayer%meltfrac
-            mixedregodata%meltvolume = surf(xstpi,ystpi)%regolayer%meltvolume
-            mixedregodata%totvolume = surf(xstpi,ystpi)%regolayer%totvolume
             call regolith_traverse_streamtube(user,surf(xstpi,ystpi),deltar,ri,rip1,eradi,erado,newlayer,vmare,&
                totseb,age_collector,xmints,xsfints,rsh,depthb,mixedregodata)
             totmare = totmare + vmare

src/regolith/regolith_streamtube_head.f90

@@ -25,7 +25,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
    ! arguemnts
    type(usertype),intent(in) :: user
    type(surftype),intent(in) :: surfi
-   type(regodatatype),intent(inout) :: regodatatype
+   type(regodatatype),intent(inout) :: mixedregodata
    real(DP),intent(in) :: deltar
    real(DP),intent(inout) :: totmare,tots
    real(SP),dimension(:),intent(inout) :: age_collector
@@ -69,50 +69,51 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
       mixedregodata%totvolume = mixedregodata%totvolume + vsgly
       mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol
       mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-      if (mixedregoda%meltfrac > 1.0_DP) then
+      if (mixedregodata%meltfrac > 1.0_DP) then
          write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
       end if
    else ! head is not intersected with layers. 
 
-   do
-      ! if (.not. associated(current%next)) exit
-
-      if (zend < zmax) then 
-         vhead = regolith_circle_sector_func(deltar,zstart,zend)
-         tothead = tothead + vhead * vratio 
-         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 (mixedregoda%meltfrac > 1.0_DP) then
-            write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
+      do
+         ! if (.not. associated(current%next)) exit
+
+         if (zend < zmax) then 
+            vhead = regolith_circle_sector_func(deltar,zstart,zend)
+            tothead = tothead + vhead * vratio 
+            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
+            !current => current%next
+            N = N - 1
+            z = z + current(N)%thickness
+            zstart = zend
+            zend = z
+         else 
+            totmarehead = totmarehead + (vsgly-tothead) * current(N)%comp
+            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
+            exit
          end if
-         !current => current%next
-         N = N - 1
-         z = z + current(N)%thickness
-         zstart = zend
-         zend = z
-      else 
-         totmarehead = totmarehead + (vsgly-tothead) * current(N)%comp
-         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 (mixedregoda%meltfrac > 1.0_DP) then
-            write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
-         end if
-         exit
-      end if
-   end do
+      end do
+
 
-   tots = tots + tothead
-   totmare = totmare + totmarehead
+      tots = tots + tothead
+      totmare = totmare + totmarehead
 
    end if
 

src/regolith/regolith_streamtube_lineseg.f90

@@ -86,7 +86,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                 linmelt = current(N-1)%meltfrac * vsgly * recyratio
                 mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
                 mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                if (mixedregoda%meltfrac > 1.0_DP) then
+                if (mixedregodata%meltfrac > 1.0_DP) then
                   write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
                end if
              else if (ri <= xmints .and. rip1 > xmints) then 
@@ -99,7 +99,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                      linmelt = current(N-1)%meltfrac * vseg * recyratio
                      mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
                      mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                     if (mixedregoda%meltfrac > 1.0_DP) then
+                     if (mixedregodata%meltfrac > 1.0_DP) then
                         write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
                      end if
              end if
@@ -132,7 +132,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                   linmelt = current(N)%meltfrac * vseg * recyratio
                   mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
                   mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                  if (mixedregoda%meltfrac > 1.0_DP) then
+                  if (mixedregodata%meltfrac > 1.0_DP) then
                      write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
                   end if
                end if
@@ -148,7 +148,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                   linmelt = current(N)%meltfrac * vseg * recyratio
                   mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
                   mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                  if (mixedregoda%meltfrac > 1.0_DP) then
+                  if (mixedregodata%meltfrac > 1.0_DP) then
                      write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
                   end if
                end if
@@ -176,7 +176,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                   linmelt = current(N)%meltfrac * vseg * recyratio
                   mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
                   mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-                  if (mixedregoda%meltfrac > 1.0_DP) then
+                  if (mixedregodata%meltfrac > 1.0_DP) then
                      write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
                   end if
                end if
@@ -189,10 +189,10 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
                vmare = vmare + (vsgly - totseb) * current(N)%comp
                totseb = vsgly
                mixedregodata%totvolume = mixedregodata%totvolume + vsgly
-               linmelt = surfi(N)%meltfrac * vsgly * recyratio
+               linmelt = current(N)%meltfrac * vsgly * recyratio
                mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt
                mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-               if (mixedregoda%meltfrac > 1.0_DP) then
+               if (mixedregodata%meltfrac > 1.0_DP) then
                   write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)"
                end if 
                exit

src/regolith/regolith_subpixel_streamtube.f90

@@ -88,9 +88,6 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer
 
    !current => surfi%regolayer
    allocate(current,source=surfi%regolayer)
-   mixedregodata%totvolume = 0
-   mixedregodata%meltvolume = 0
-   mixedregodata%meltfrac = 0
    N = size(current)
    z = surfi%regolayer(N)%thickness
    zstart = 0.0_DP
@@ -202,7 +199,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer
          mixedregodata%meltvolume = mixedregodata%meltvolume + mvl + mvr
          mixedregodata%totvolume = mixedregodata%totvolume + vsgly
          mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume
-         if (mixedregoda%meltfrac > 1.0_DP) then
+         if (mixedregodata%meltfrac > 1.0_DP) then
             write(*,*) "ERROR! mixedregodata%meltfrac >1! (SUBPIXEL)"
          end if
          N = N - 1