melt tracking WIP sent to cluster for debugging (does not compile)

src/globals/module_globals.f90

@@ -72,6 +72,10 @@ module module_globals
    real(DP) :: porosity   ! Porosity: Maximum 1, Minium 0.
    real(DP) :: damage     ! Damage  : Maximum 1, Minium 0.
    real(DP) :: depth      ! Absolute location with respect to the initial surface. 
+   real(DP) :: meltvolume
+   real(DP) :: totvolume
+   real(DP) :: ejm !ejected melt
+   real(DP) :: ejmf !ejected melt fraction
    real(SP),dimension(:),allocatable :: meltdist !its dimension should be the number of quasimc craters
 end type regodatatype
 

src/regolith/module_regolith.f90

@@ -95,13 +95,13 @@ 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)
+                 age_collector,xmints,xsfints,rsh,depthb,mixedregodata)
       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
+      type(regodatatype),intent(inout) :: newlayer, mixedregodata
       real(DP),intent(out) :: vmare,totseb
       real(SP),dimension(:),intent(inout) :: age_collector
       real(DP),intent(in)             :: xmints
@@ -111,13 +111,13 @@ end subroutine regolith_traverse_streamtube
 
    interface 
       subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer,vmare,totseb,& 
-                 age_collector,xmints,xsfints,vol) 
+                 age_collector,xmints,xsfints,vol,mixedregodata) 
       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
+      type(regodatatype),intent(inout)    :: newlayer, mixedregodata
       real(DP),intent(out)                :: vmare,totseb
       real(SP),dimension(:),intent(inout) :: age_collector
       real(DP),intent(in)                 :: xmints
@@ -128,13 +128,13 @@ 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)
+      totseb,age_collector,xmints,xsfints,depthb,mixedregodata)
       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
+      type(regodatatype),intent(inout) :: newlayer,mixedregodata
       real(DP),intent(inout) :: vmare,totseb
       real(SP),dimension(:),intent(inout) :: age_collector
       real(DP),intent(in)             :: xmints
@@ -143,11 +143,12 @@ end subroutine regolith_streamtube_lineseg
    end interface 
 
    interface 
-      subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector)
+      subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,mixedregodata)
       use module_globals 
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),intent(in) :: surfi
+      type(regodatatype),intent(inout) :: mixedregodata
       real(DP),intent(in) :: deltar
       real(DP),intent(inout) :: totmare,tots
       real(SP),dimension(:),intent(inout) :: age_collector

src/regolith/regolith_melt_glass.f90

@@ -130,6 +130,10 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
    newlayer%meltfrac  = 0.0_DP ! default value: no melt (zero)
    newlayer%thickness = ebh    ! default value: stream tube's volume = paraboloid shell's volume
    newlayer%comp      = 0.0_DP
+   newlayer%meltvolume = 0.0_DP
+   newlayer%totvolume = 0.0_DP
+   newlayer%ejm       = 0.0_DP
+   newlayer%ejmf      = 0.0_DP
    rints              = sqrt(rm**2 - (crater%imp/2.0)**2)
    cosvints           = min(max(eradi / (crater%imp + eradi), -1.0_DP), 1.0_DP)
    sinvints           = sqrt(1.0 - cosvints**2)
@@ -140,6 +144,8 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
    if (eradi <= rints) then
       volm1    = vst - volv1
       melt     = volm1
+      newlayer%meltvolume = melt
+      newlayer%totvolume = melt
       newlayer%meltfrac = 1.0
       xmints   = rints 
    else if (eradi > rints) then
@@ -164,14 +170,18 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
            xmints       = eradi * (1.0 - cos(thetaq)) * sin(thetaq)
            volm1        = regolith_streamtube_volume_func(eradi,0.0_DP,xmints,deltar)
            melt         = volm1 - volv1
+           newlayer%meltvolume = melt
+           newlayer%totvolume = vst-volv1
            newlayer%meltfrac = melt/(vst-volv1)
+           newlayer%ejm = melt
+           newlayer%ejmf = newlayer%meltfrac
 
    end if
 
    allocate(newlayer%meltdist((domain%rcnum)))
    newlayer%meltdist(:) = 0.0_SP
    if(domain%currentqmc) then
-      newlayer%meltdist(domain%nqmc) = newlayer%meltfrac
+      newlayer%meltdist(domain%nqmc) = newlayer%meltfrac !might change this to melt
    end if 
 
    n_age = max(ceiling(age / age_resolution), 1)

src/regolith/regolith_streamtube.f90

@@ -89,7 +89,7 @@ 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
+   type(regodatatype) :: newlayer, mixedregodata
 
    ! Constrain the tangital tube's volume with CTEM result
    real(DP)     :: k1,k2,k3,k4,c1,c2
@@ -222,7 +222,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
    ! Purpose 2: Once we have the size information of a stream tube, we can
    ! calculate the distal melt: the precursor of glass spherules within a
    ! stream tube. The result is contained in a linked list "newlayer".
-   call regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints, volm)
+   call regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,volm)
    ! if (eradc>rm) then
    !    write(*,*) 'eradc > rm!'
    !    write(*,*) ebh, exp(ejb(k)%thick)
@@ -254,8 +254,11 @@ 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)
+           age_collector,xmints,xsfints,vol,mixedregodata)
 
       newlayer%age(:) = newlayer%age(:) + age_collector(:)
 
@@ -264,6 +267,9 @@ 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%meltfrac = newlayer%meltvolume / newlayer%totvolume
 
    else
       rbody = sqrt(xints(2)**2 + yints(2)**2)
@@ -275,6 +281,9 @@ 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)
          totmare = totmare + vmare
@@ -292,8 +301,11 @@ 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)
+               totseb,age_collector,xmints,xsfints,rsh,depthb,mixedregodata)
             totmare = totmare + vmare
             tots = tots + totseb 
          end if
@@ -302,12 +314,15 @@ 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)
+      call regolith_streamtube_head(user,surf(xstpi,ystpi),deltar,totmare,tots,age_collector,mixedregodata)
 
       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%meltfrac = newlayer%meltvolume / newlayer%totvolume
 
   end if
 

src/regolith/regolith_streamtube_head.f90

@@ -18,13 +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)
+subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,mixedregodata)
    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) :: regodatatype
    real(DP),intent(in) :: deltar
    real(DP),intent(inout) :: totmare,tots
    real(SP),dimension(:),intent(inout) :: age_collector
@@ -42,6 +43,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
 
    ! melt collector
    real(DP) :: recyratio
+   real(DP) :: headmeltvol
 
    !current => surfi%regolayer
    !current = surfi%regolayer
@@ -63,6 +65,13 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
       totmare = totmare + vsgly * current(N)%comp
       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 (mixedregoda%meltfrac > 1.0_DP) then
+         write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
+      end if
    else ! head is not intersected with layers. 
 
    do
@@ -74,6 +83,13 @@ 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 (mixedregoda%meltfrac > 1.0_DP) then
+            write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
+         end if
          !current => current%next
          N = N - 1
          z = z + current(N)%thickness
@@ -84,6 +100,13 @@ 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 (mixedregoda%meltfrac > 1.0_DP) then
+            write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)"
+         end if
          exit
       end if
    end do