apply mixing volume changes to allocatable_arrays branch

.gitignore

@@ -76,3 +76,5 @@ python/ctem/ctem.egg-info/
 *.i90
 
 INSTALL
+
+examples/global-lunar-bombardment/temp.in

examples/global-lunar-bombardment/qmctest2.in

@@ -0,0 +1,8 @@
+#Dcrat(m)   vel(m/s)   ang(deg)  xoffset(m)  yoffset(m)  t(Ga)
+2050000 18300.0 45.0 -2908888.89 -1607137.57 0.6
+658628 18300.0 45.0 570078.99 770213.11 0.5
+1012000 18300.0 45.0 -297849.189 1121963.97 0.2
+623000 18300.0 45.0 1064349.59 -473044.26 0.4
+809000 18300.0 45.0 1670883.00 509433.00 0.3
+639000 18300.0 45.0 -1874653.62 -609499.34 0.1
+233931 18300.0 45.0 -533541.38 1523910.2 0.05

examples/global-lunar-bombardment/qmctest3.in

@@ -0,0 +1,3 @@
+#Dcrat(m)   vel(m/s)   ang(deg)  xoffset(m)  yoffset(m)  t(Ga)
+886977 18300.0 45.0 1482572.48 -62471.4873 0.02
+10000 18300.0 45.0 70000.0 10000.0 0.006

src/regolith/regolith_mix.f90

@@ -49,8 +49,8 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain)
    newlayer%age(:)    = 0.0_DP
    allocate(newlayer%meltdist(1+domain%rcnum))
    allocate(newlayer%distvol(1+domain%rcnum))
-   newlayer%meltdist(:) = 0.0_SP
-   newlayer%distvol(:) = 0.0_SP
+   newlayer%meltdist(:) = 0.0_DP
+   newlayer%distvol(:) = 0.0_DP
    newlayer%ejm       = 0.0_DP
    newlayer%ejmf      = 0.0_DP
    newlayer%meltvolume = 0.0_DP
@@ -62,27 +62,15 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain)
    do i = N,1,-1
       newlayer%thickness = newlayer%thickness + poppedarray(i)%thickness
       newlayer%comp      = newlayer%comp + poppedarray(i)%thickness * poppedarray(i)%comp       
-      !newlayer%meltfrac  = newlayer%meltfrac + poppedarray(i)%thickness * poppedarray(i)%meltfrac
       newlayer%age(:)    = newlayer%age(:) + poppedarray(i)%age(:)
-      !newlayer%meltdist(:) = newlayer%meltdist(:) + poppedarray(i)%thickness * poppedarray(i)%meltdist(:)
-      newlayer%distvol(:) = newlayer%distvol(:) + poppedarray(i)%thickness * poppedarray(i)%distvol(:)
-      newlayer%ejm       = newlayer%ejm + poppedarray(i)%thickness * poppedarray(i)%ejm
-      !newlayer%ejmf      = newlayer%ejmf + poppedarray(i)%thickness * poppedarray(i)%ejmf
-      newlayer%meltvolume = newlayer%meltvolume + poppedarray(i)%thickness * poppedarray(i)%meltvolume
-      !newlayer%totvolume = newlayer%totvolume + poppedarray(i)%thickness * poppedarray(i)%totvolume
+      newlayer%distvol(:) = newlayer%distvol(:) + poppedarray(i)%distvol(:)
+      newlayer%ejm       = newlayer%ejm + poppedarray(i)%ejm
+      newlayer%meltvolume = newlayer%meltvolume + poppedarray(i)%meltvolume
    end do
 
    ! Get average values of composition and melt fraction
    newlayer%comp = newlayer%comp / newlayer%thickness 
-   !newlayer%meltfrac = newlayer%meltfrac / newlayer%thickness 
-   !newlayer%meltdist(:) = newlayer%meltdist(:) / newlayer%thickness
-   newlayer%distvol(:) = newlayer%distvol(:) / newlayer%thickness
-   newlayer%ejm = newlayer%ejm / newlayer%thickness
-   !newlayer%ejmf = newlayer%ejmf / newlayer%thickness
-   newlayer%meltvolume = newlayer%meltvolume / newlayer%thickness
-   !newlayer%totvolume = newlayer%totvolume / newlayer%thickness
 
-   !test forcing melt fraction to equal melt volume / total volume
    newlayer%totvolume = newlayer%thickness * user%pix * user%pix
    newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume
    newlayer%meltdist(:) = newlayer%distvol(:) / newlayer%totvolume
@@ -91,5 +79,12 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain)
    call util_push_array(surfi%regolayer, newlayer)
    !call util_destroy_list(poppedlist_top)
 
+
+   ! do i = N,1,-1
+   !    if (abs(surfi%regolayer(i)%meltvolume - sum(surfi%regolayer(i)%distvol) > 1e-5)) then
+   !    write(*,*) "melt array =/= melt value!", domain%nqmc, domain%currentqmc, abs(surfi%regolayer(i)%meltvolume - sum(surfi%regolayer(i)%distvol))
+   !    end if
+   ! end do
+
    return
 end subroutine regolith_mix

src/util/util_traverse_pop_array.f90

@@ -65,19 +65,29 @@ subroutine util_traverse_pop_array(user,regolayer,traverse_depth,poppedarray)
 
 
     allocate(poppedarray,source=regolayer(maxi:N))
+    allocate(oldregodata,source=regolayer(maxi:maxi))
 
     !for #1 element of poppedarray, shrink thickness by whatever was lefr over. In corresponding maxi of regolayer, also need to change that.
 
     poppedarray(1)%thickness = poppedarray(1)%thickness - depth_diff
     regolayer(maxi)%thickness = regolayer(maxi)%thickness - poppedarray(1)%thickness
+    poppedarray(1)%meltvolume = (poppedarray(1)%thickness/oldregodata(1)%thickness) * oldregodata(1)%meltvolume
+    poppedarray(1)%distvol(:) = (poppedarray(1)%thickness/oldregodata(1)%thickness) * oldregodata(1)%distvol(:)
+    poppedarray(1)%ejm = (poppedarray(1)%thickness/oldregodata(1)%thickness) * oldregodata(1)%ejm
+    regolayer(maxi)%meltvolume = (regolayer(maxi)%thickness/oldregodata(1)%thickness) * oldregodata(1)%meltvolume
+    regolayer(maxi)%distvol(:) = (regolayer(maxi)%thickness/oldregodata(1)%thickness) * oldregodata(1)%distvol(:)
+    regolayer(maxi)%ejm = (regolayer(maxi)%thickness/oldregodata(1)%thickness) * oldregodata(1)%ejm
+
     poppedarray(1)%totvolume = poppedarray(1)%thickness * user%pix * user%pix
-    poppedarray(1)%meltvolume = poppedarray(1)%meltfrac * poppedarray(1)%totvolume
-    poppedarray(1)%distvol(:) = poppedarray(1)%meltdist(:) * poppedarray(1)%totvolume
-    poppedarray(1)%ejm = poppedarray(1)%ejmf * poppedarray(1)%totvolume
+    if(poppedarray(1)%totvolume > 0) then
+        poppedarray(1)%meltfrac = poppedarray(1)%meltvolume / poppedarray(1)%totvolume
+        poppedarray(1)%ejmf = poppedarray(1)%ejm / poppedarray(1)%totvolume
+    end if
     regolayer(maxi)%totvolume = regolayer(maxi)%thickness * user%pix * user%pix
-    regolayer(maxi)%distvol(:) = regolayer(maxi)%meltdist(:) * regolayer(maxi)%totvolume
-    regolayer(maxi)%meltvolume = regolayer(maxi)%meltfrac * regolayer(maxi)%totvolume
-    regolayer(maxi)%ejm = regolayer(maxi)%ejmf * regolayer(maxi)%totvolume
+    regolayer(maxi)%meltfrac = regolayer(maxi)%meltvolume / regolayer(maxi)%totvolume
+    regolayer(maxi)%ejmf = regolayer(maxi)%ejm / regolayer(maxi)%totvolume
+
+    deallocate(oldregodata)
 
     ! copy regolayer from 1 to maxi to temp variable, then deallocate regolayer, then movealloc templayer onto regolayer <--may need temp array
     allocate(oldregodata,source=regolayer(1:maxi))