Local melt is now a part of the meltdist and distvol arrays

src/globals/module_globals.f90

@@ -75,8 +75,8 @@ module module_globals
    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
-   real(SP),dimension(:),allocatable :: distvol !its dimension should be the number of quasimc craters
+   real(SP),dimension(:),allocatable :: meltdist !its dimension should be the number of quasimc craters + 1
+   real(SP),dimension(:),allocatable :: distvol !its dimension should be the number of quasimc craters + 1
 end type regodatatype
 
 type regolisttype

src/init/init_regolith_stack.f90

@@ -45,8 +45,8 @@ subroutine init_regolith_stack(user,surf,domain)
    bedrock%meltfrac  = 0._DP 
    bedrock%comp      = 0._DP
    bedrock%age(:)    = 0.0_SP
-   allocate(bedrock%meltdist(domain%rcnum))
-   allocate(bedrock%distvol(domain%rcnum))
+   allocate(bedrock%meltdist(1+domain%rcnum))
+   allocate(bedrock%distvol(1+domain%rcnum))
    bedrock%meltdist(:) = 0.0_SP
    bedrock%distvol(:) = 0.0_SP
    bedrock%meltvolume = 0.0_DP

src/io/io_read_regotrack.f90

@@ -120,17 +120,17 @@ subroutine io_read_regotrack(user,surf,domain)
 
    ! Start pushing regolith thickness and melt fraction of each layer in FILO manner
 
-   allocate(newsurfi%meltdist(domain%rcnum))
-   allocate(newsurfi%distvol(domain%rcnum))
+   allocate(newsurfi%meltdist(1+domain%rcnum))
+   allocate(newsurfi%distvol(1+domain%rcnum))
 
    do j=1,user%gridsize
       do i=1,user%gridsize
 
          !call util_init_list(surf(i,j)%regolayer,initstat)
          call util_init_array(user,surf(i,j)%regolayer,domain,initstat)
          N = stacks_num(i,j)
-         allocate(meltvolume(N),thickness(N),comp(N),age(MAXAGEBINS,N),distvol(domain%rcnum,N),ejm(N),ejmf(N),&
-            meltfrac(N),meltdist(domain%rcnum,N))
+         allocate(meltvolume(N),thickness(N),comp(N),age(MAXAGEBINS,N),distvol(1+domain%rcnum,N),ejm(N),ejmf(N),&
+            meltfrac(N),meltdist(1+domain%rcnum,N))
 
          read(FMELT) meltvolume(:)
          read(FREGO) thickness(:)
@@ -167,7 +167,7 @@ subroutine io_read_regotrack(user,surf,domain)
             do q=1,MAXAGEBINS
                newsurfi%age(q) = agei(MAXAGEBINS*k-(MAXAGEBINS-q))
             end do
-            do q=1,domain%rcnum
+            do q=1,1+domain%rcnum
                newsurfi%meltdist(q) = meltdist(q,k)
                newsurfi%distvol(q) = distvol(q,k)
             end do

src/io/io_write_regotrack.f90

@@ -82,8 +82,8 @@ subroutine io_write_regotrack(user,surf,domain)
       do i=1,user%gridsize
          !current => surf(i,j)%regolayer
          N = stacks_num(i,j)
-         allocate(meltvolume(N),thickness(N),comp(N),age(MAXAGEBINS,N),distvol(domain%rcnum,N),ejm(N),ejmf(N),&
-            meltfrac(N),meltdist(domain%rcnum,N))
+         allocate(meltvolume(N),thickness(N),comp(N),age(MAXAGEBINS,N),distvol(1+domain%rcnum,N),ejm(N),ejmf(N),&
+            meltfrac(N),meltdist(1+domain%rcnum,N))
          allocate(current,source=surf(i,j)%regolayer)
          do k=1,N
             meltfrac(k) = current(k)%meltfrac

src/regolith/regolith_interior.f90

@@ -41,8 +41,8 @@ subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltshee
     !Executable code
 
     hmeltsheet = vmeltsheet / (nmeltsheet*user%pix*user%pix)
-    allocate(newlayer%meltdist(domain%rcnum))
-    allocate(newlayer%distvol(domain%rcnum))
+    allocate(newlayer%meltdist(1+domain%rcnum))
+    allocate(newlayer%distvol(1+domain%rcnum))
 
     inc = incval
 
@@ -77,6 +77,9 @@ subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltshee
             if(domain%currentqmc) then
                 newlayer%meltdist(domain%nqmc) = newlayer%meltfrac
                 newlayer%distvol(domain%nqmc) = newlayer%meltvolume
+            else
+                newlayer%meltdist(1+domain%rcnum) = newlayer%meltfrac
+                newlayer%distvol(1+domain%rcnum) = newlayer%meltvolume
             end if
             call util_push_array(surf(xpi,ypi)%regolayer,newlayer)
         end do

src/regolith/regolith_melt_glass.f90

@@ -180,9 +180,9 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
 
    end if
 
-   allocate(newlayer%meltdist((domain%rcnum)))
+   allocate(newlayer%meltdist((1+domain%rcnum)))
    newlayer%meltdist(:) = 0.0_SP
-   allocate(newlayer%distvol((domain%rcnum)))
+   allocate(newlayer%distvol((1+domain%rcnum)))
    newlayer%distvol(:) = 0.0_SP
    if(domain%currentqmc) then
       newlayer%meltdist(domain%nqmc) = newlayer%meltfrac

src/regolith/regolith_mix.f90

@@ -47,8 +47,8 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain)
    newlayer%comp      = 0.0_DP
    newlayer%meltfrac  = 0.0_DP
    newlayer%age(:)    = 0.0_DP
-   allocate(newlayer%meltdist(domain%rcnum))
-   allocate(newlayer%distvol(domain%rcnum))
+   allocate(newlayer%meltdist(1+domain%rcnum))
+   allocate(newlayer%distvol(1+domain%rcnum))
    newlayer%meltdist(:) = 0.0_SP
    newlayer%distvol(:) = 0.0_SP
    newlayer%ejm       = 0.0_DP

src/regolith/regolith_streamtube.f90

@@ -239,7 +239,7 @@ 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))
+   allocate(distvol(1+domain%rcnum))
    distvol(:) = 0.0_SP
 
    ! if (eradc<=user%testimp) then
@@ -355,7 +355,18 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
       end if
       newlayer%distvol(:) = newlayer%distvol(:) + distvol(:)
       newlayer%meltdist(:) = newlayer%distvol(:) / newlayer%totvolume
+
+      newlayer%distvol(1+domain%rcnum) = newlayer%meltvolume - sum(newlayer%distvol(1:domain%rcnum))
+      if (newlayer%distvol(1+domain%rcnum) < 0.0_DP) then !pixel consists entirely of QMC melt
+         newlayer%distvol(1+domain%rcnum) = 0.0_SP
+      end if
+      if (sum(newlayer%distvol) > newlayer%totvolume) then
+         factor = newlayer%totvolume / sum(newlayer%distvol)
+         newlayer%distvol(:) = newlayer%distvol(:) * factor
+      end if
+      newlayer%meltvolume = sum(newlayer%distvol)
       newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume
+
    end if
 
   call util_push_array(surf(xpi,ypi)%regolayer,newlayer)

src/util/util_init_array.f90

@@ -56,9 +56,9 @@ subroutine util_init_array(user,regolayer,domain,initstat)
         regolayer(1)%age(:)   = 0.0_SP
         regolayer(1)%ejm = 0.0_DP
         regolayer(1)%ejmf = 0.0_DP
-        allocate(regolayer(1)%meltdist(domain%rcnum))
+        allocate(regolayer(1)%meltdist(1+domain%rcnum))
         regolayer(1)%meltdist(:) = 0.0_SP
-        allocate(regolayer(1)%distvol(domain%rcnum))
+        allocate(regolayer(1)%distvol(1+domain%rcnum))
         regolayer(1)%distvol(:) = 0.0_SP
         regolayer(1)%meltvolume = 0.0_DP
         regolayer(1)%totvolume = regolayer(1)%thickness * user%pix * user%pix