removed meltfrac as a variable and calculate melt entirely through volume (currently untested)

python/ctem/ctem/driver.py

@@ -336,7 +336,7 @@ def process_output(self):
                 self.redirect_outputs(['tcum'], 'dist')
             self.redirect_outputs(['impmass'], 'misc')
             if (self.user['saverego'].upper() == 'T') :
-               self.redirect_outputs(['stack','rego','age','melt','comp','distfrac','ejm','ejmf','meltdist','meltfrac'], 'rego')
+               self.redirect_outputs(['stack','rego','age','melt','comp','ejm','meltdist'], 'rego')
 
 
 

src/globals/module_globals.f90

@@ -66,16 +66,13 @@ module module_globals
 type regodatatype 
    real(SP),dimension(MAXAGEBINS) :: age 
    real(DP) :: thickness
-   real(DP) :: meltfrac 
    real(DP) :: comp 
    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 + 1
    real(SP),dimension(:),allocatable :: distvol !its dimension should be the number of quasimc craters + 1
 end type regodatatype
 
@@ -315,9 +312,6 @@ module module_globals
 character(*),parameter :: RCFILE     = 'craterlist.dat'
 character(*),parameter :: MDFILE     = 'surface_meltdist.dat'
 character(*),parameter :: EJMFILE    = 'surface_ejm.dat'
-character(*),parameter :: EJMFFILE   = 'surface_ejmf.dat'
-character(*),parameter :: MELTFRACFILE   = 'surface_meltfrac.dat'
-character(*),parameter :: DISTFRACFILE = 'surface_distfrac.dat'
 
 ! Global variables 
 integer(I4B),parameter :: PBCLIM = 1             ! periodic boundary condition limit

src/init/init_regolith_stack.f90

@@ -42,17 +42,13 @@ subroutine init_regolith_stack(user,surf,domain)
    ! Initialize the grid space  
    !=======================================
    bedrock%thickness = user%trad
-   bedrock%meltfrac  = 0._DP 
    bedrock%comp      = 0._DP
    bedrock%age(:)    = 0.0_SP
-   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
    bedrock%totvolume = bedrock%thickness * user%pix * user%pix
    bedrock%ejm = 0.0_DP
-   bedrock%ejmf = 0.0_DP
 
    do yp = 1, user%gridsize
       do xp = 1, user%gridsize

src/io/io_read_regotrack.f90

@@ -35,15 +35,12 @@ subroutine io_read_regotrack(user,surf,domain)
    integer(I4B), parameter :: FAGE = 13
    integer(I4B), parameter :: FMD = 14
    integer(I4B), parameter :: FEJM = 15
-   integer(I4B), parameter :: FEJMF = 16
-   integer(I4B), parameter :: FMF = 17
-   integer(I4B), parameter :: FDF = 18
    ! real(DP),dimension(user%gridsize,user%gridsize) :: regotop,melt,comp,ejm,ejmf,meltfrac
    ! real(SP),dimension(user%gridsize,user%gridsize,domain%rcnum) :: meltdist, distfrac
    ! real(SP),dimension(user%gridsize,user%gridsize,MAXAGEBINS) :: age
    integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num 
-   real(DP),dimension(:),allocatable :: regotop,melt,comp,ejm,ejmf,meltfrac,thickness,meltvolume,agei
-   real(SP),dimension(:,:),allocatable :: age, meltdist, distfrac, distvol
+   real(DP),dimension(:),allocatable :: regotop,melt,comp,ejm,thickness,meltvolume,agei
+   real(SP),dimension(:,:),allocatable :: age, distvol
    !real(DP), dimension(:), allocatable :: regotopi,melti,compi,agei,dfi,ejmi,ejmfi,mdi,mfi
    type(regodatatype) :: newsurfi
    integer(I4B) :: ioerr,i,j,k,q,itmp,N
@@ -82,35 +79,18 @@ subroutine io_read_regotrack(user,surf,domain)
        stop
    end if   
 
-   open(FDF,file=DISTFRACFILE,status='old',form='unformatted',iostat=ioerr)
-   if (ioerr/=0) then 
-       write(*,*) 'Error! Cannot read file ',trim(adjustl(DISTFRACFILE))
-       stop
-   end if  
-
    open(FEJM,file=EJMFILE,status='old',form='unformatted',iostat=ioerr)
    if (ioerr/=0) then 
        write(*,*) 'Error! Cannot read file ',trim(adjustl(EJMFILE))
        stop
    end if
-
-   open(FEJMF,file=EJMFFILE,status='old',form='unformatted',iostat=ioerr)
-   if (ioerr/=0) then 
-       write(*,*) 'Error! Cannot read file ',trim(adjustl(EJMFFILE))
-       stop
-   end if  
 
    open(FMD,file=MDFILE,status='old',form='unformatted',iostat=ioerr)
    if (ioerr/=0) then 
        write(*,*) 'Error! Cannot read file ',trim(adjustl(MDFILE))
        stop
    end if  
 
-   open(FMF,file=MELTFRACFILE,status='old',form='unformatted',iostat=ioerr)
-   if (ioerr/=0) then 
-       write(*,*) 'Error! Cannot read file ',trim(adjustl(MELTFRACFILE))
-       stop
-   end if  
 
    open(FAGE,file=AGEFILE,status='old',form='unformatted',iostat=ioerr)
    if (ioerr/=0) then
@@ -119,7 +99,6 @@ subroutine io_read_regotrack(user,surf,domain)
    end if
 
    ! Start pushing regolith thickness and melt fraction of each layer
-   allocate(newsurfi%meltdist(1+domain%rcnum))
    allocate(newsurfi%distvol(1+domain%rcnum))
 
    do j=1,user%gridsize
@@ -128,18 +107,14 @@ subroutine io_read_regotrack(user,surf,domain)
          !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(1+domain%rcnum,N),ejm(N),ejmf(N),&
-            meltfrac(N),meltdist(1+domain%rcnum,N))
+         allocate(meltvolume(N),thickness(N),comp(N),age(MAXAGEBINS,N),distvol(1+domain%rcnum,N),ejm(N))
 
          read(FMELT) meltvolume(:)
          read(FREGO) thickness(:)
          read(FCOMP) comp(:)
          read(FAGE) age(:,:)
          read(FMD) distvol(:,:)
          read(FEJM) ejm(:)
-         read(FEJMF) ejmf(:)
-         read(FMF) meltfrac(:)
-         read(FDF) meltdist(:,:)
 
          allocate(agei(MAXAGEBINS * stacks_num(i,j)))
 
@@ -156,33 +131,27 @@ subroutine io_read_regotrack(user,surf,domain)
          do k=1,max(stacks_num(i,j),1),1
             newsurfi%thickness = thickness(k)
             newsurfi%comp = comp(k)
-            newsurfi%meltfrac  = meltfrac(k)
             newsurfi%meltvolume = meltvolume(k)
             newsurfi%ejm = ejm(k)
-            newsurfi%ejmf = ejmf(k)
             newsurfi%totvolume = thickness(k) * user%gridsize * user%gridsize
             do q=1,MAXAGEBINS
                newsurfi%age(q) = agei(MAXAGEBINS*k-(MAXAGEBINS-q))
             end do
             do q=1,1+domain%rcnum
-               newsurfi%meltdist(q) = meltdist(q,k)
                newsurfi%distvol(q) = distvol(q,k)
             end do
             call util_push_array(surf(i,j)%regolayer,newsurfi)
          end do
 
-         deallocate(meltvolume,thickness,comp,age,distvol,ejm,ejmf,meltfrac,meltdist,agei)
+         deallocate(meltvolume,thickness,comp,age,distvol,ejm,agei)
 
       end do
    end do
    close(FMELT)
    close(FREGO)
    close(FCOMP)
    close(FAGE)
-   close(FDF)
    close(FEJM)
-   close(FEJMF)
    close(FMD)
-   close(FMF)
    return
 end subroutine io_read_regotrack

src/io/io_write_regotrack.f90

@@ -35,14 +35,11 @@ subroutine io_write_regotrack(user,surf,domain)
    integer(I4B), parameter :: FAGE = 13
    integer(I4B), parameter :: FMD = 14
    integer(I4B), parameter :: FEJM = 15
-   integer(I4B), parameter :: FEJMF = 16
-   integer(I4B), parameter :: FMF = 17
-   integer(I4B), parameter :: FDF = 18
    !type(regolisttype),pointer :: current => null()
    type(regodatatype),dimension(:),allocatable :: current
    integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num
-   real(DP),dimension(:),allocatable :: meltfrac, thickness, comp, ejm, ejmf, meltvolume
-   real(SP),dimension(:,:),allocatable :: age, meltdist, distvol
+   real(DP),dimension(:),allocatable :: thickness, comp, ejm, meltvolume
+   real(SP),dimension(:,:),allocatable :: age, distvol
    integer(kind=8) :: recsize
    real(DP) :: dtmp
    real(SP) :: stmp
@@ -57,9 +54,6 @@ subroutine io_write_regotrack(user,surf,domain)
    open(FAGE,file=AGEFILE,status='replace',form='unformatted')
    open(FMD,file=MDFILE,status='replace',form='unformatted')
    open(FEJM,file=EJMFILE,status='replace',form='unformatted')
-   open(FMF,file=MELTFRACFILE,status='replace',form='unformatted')
-   open(FDF,file=DISTFRACFILE,status='replace',form='unformatted')
-   open(FEJMF,file=EJMFFILE,status='replace',form='unformatted')
 
    ! First pass to get stack numbers
    stacks_num(:,:) = 0
@@ -69,11 +63,6 @@ subroutine io_write_regotrack(user,surf,domain)
          allocate(current,source=surf(i,j)%regolayer)
          stacks_num(i,j) = size(current)
          deallocate(current)
-         ! do 
-         !    if (.not. associated(current)) exit ! We've reached the bottom of the linked list
-         !    stacks_num(i,j) = stacks_num(i,j) + 1
-         !    current => current%next
-         ! end do
       end do 
    end do
 
@@ -82,24 +71,16 @@ 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(1+domain%rcnum,N),ejm(N),ejmf(N),&
-            meltfrac(N),meltdist(1+domain%rcnum,N))
+         allocate(meltvolume(N),thickness(N),comp(N),age(MAXAGEBINS,N),distvol(1+domain%rcnum,N),ejm(N))
          allocate(current,source=surf(i,j)%regolayer)
          do k=1,N
-            meltfrac(k) = current(k)%meltfrac
-            ! thickness(k) = current%regodata%thickness
-            ! comp(k) = current%regodata%comp
-            ! age(:,k) = current%regodata%age(:)
-            ! current => current%next
             meltvolume(k) = current(k)%meltvolume
             thickness(k) = current(k)%thickness
             comp(k) = current(k)%comp
             age(:,k) = current(k)%age(:)
             !write(*,*) i, j
             distvol(:,k) = current(k)%distvol(:)
-            meltdist(:,k) = current(k)%meltdist(:)
             ejm(k) = current(k)%ejm
-            ejmf(k) = current(k)%ejmf
          end do
          deallocate(current)
          write(FMELT) meltvolume(:)
@@ -108,10 +89,7 @@ subroutine io_write_regotrack(user,surf,domain)
          write(FAGE) age(:,:)
          write(FMD) distvol(:,:)
          write(FEJM) ejm(:)
-         write(FEJMF) ejmf(:)
-         write(FMF) meltfrac(:)
-         write(FDF) meltdist(:,:)
-         deallocate(meltvolume,thickness,comp,age,distvol,ejm,ejmf,meltfrac,meltdist)
+         deallocate(meltvolume,thickness,comp,age,distvol,ejm)
       end do 
    end do
    close(FMELT)
@@ -120,9 +98,6 @@ subroutine io_write_regotrack(user,surf,domain)
    close(FAGE)
    close(FMD)
    close(FEJM)
-   close(FEJMF)
-   close(FMF)
-   close(FDF)
 
    recsize = sizeof(itmp) * user%gridsize * user%gridsize
    open(LUN,file=STACKNUMFILE,status='replace',form='unformatted',recl=recsize,access='direct')

src/regolith/regolith_interior.f90

@@ -41,7 +41,6 @@ subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltshee
     !Executable code
 
     hmeltsheet = vmeltsheet / (nmeltsheet*user%pix*user%pix)
-    allocate(newlayer%meltdist(1+domain%rcnum))
     allocate(newlayer%distvol(1+domain%rcnum))
 
     inc = incval
@@ -66,26 +65,20 @@ subroutine regolith_interior(user,surf,crater,domain,incval,nmeltsheet,vmeltshee
             deallocate(poppedarray)
 
             !fill top layer with melt sheet of given thickness hmeltsheet
-            newlayer%meltfrac = 1.0_DP
             newlayer%ejm = 0.0_DP
-            newlayer%ejmf = 0.0_DP
             newlayer%thickness = hmeltsheet
             newlayer%meltvolume = vmeltsheet / nmeltsheet
             newlayer%totvolume = newlayer%meltvolume
-            newlayer%meltdist(:) = 0.0_SP
             newlayer%distvol(:) = 0.0_SP
             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
     end do
 
-    deallocate(newlayer%meltdist)
     deallocate(newlayer%distvol)
 
     return

src/regolith/regolith_melt_glass.f90

@@ -127,13 +127,11 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
 
    vst = (0.25_DP * PI * deltar**2 * a**2 * eradi / b *(tan(b)-b)) + sqrt(2.0_DP)/2.0_DP*PI*deltar**3
 
-   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 = newlayer%thickness * user%pix * user%pix
    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)
@@ -146,9 +144,7 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
       melt     = volm1
       newlayer%meltvolume = melt
       !newlayer%totvolume = volm1
-      newlayer%meltfrac = 1.0_DP
       newlayer%ejm = melt
-      newlayer%ejmf = 1.0_DP
       xmints   = rints 
    else if (eradi > rints) then
            depthb = crater%imp / 2.0
@@ -174,18 +170,13 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
            melt         = volm1 - volv1
            newlayer%meltvolume = melt
            newlayer%totvolume = newlayer%thickness * user%pix * user%pix
-           newlayer%meltfrac = melt / newlayer%totvolume
            newlayer%ejm = melt
-           newlayer%ejmf = newlayer%meltfrac
 
    end if
 
-   allocate(newlayer%meltdist((1+domain%rcnum)))
-   newlayer%meltdist(:) = 0.0_SP
    allocate(newlayer%distvol((1+domain%rcnum)))
    newlayer%distvol(:) = 0.0_SP
    if(domain%currentqmc) then
-      newlayer%meltdist(domain%nqmc) = newlayer%meltfrac
       newlayer%distvol(domain%nqmc) = newlayer%meltvolume
    end if 
 

src/regolith/regolith_mix.f90

@@ -45,14 +45,10 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain)
 
    newlayer%thickness = 0.0_DP
    newlayer%comp      = 0.0_DP
-   newlayer%meltfrac  = 0.0_DP
    newlayer%age(:)    = 0.0_DP
-   allocate(newlayer%meltdist(1+domain%rcnum))
    allocate(newlayer%distvol(1+domain%rcnum))
-   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
    newlayer%totvolume = 0.0_DP
 
@@ -72,9 +68,6 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain)
    newlayer%comp = newlayer%comp / newlayer%thickness 
 
    newlayer%totvolume = newlayer%thickness * user%pix * user%pix
-   newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume
-   newlayer%meltdist(:) = newlayer%distvol(:) / newlayer%totvolume
-   newlayer%ejmf = newlayer%ejm / newlayer%totvolume
 
    call util_push_array(surfi%regolayer, newlayer)
    !call util_destroy_list(poppedlist_top)