From e2c3e9080228ce9e9532e097487e625d4a3cab16 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Fri, 9 Jun 2023 18:07:31 -0400 Subject: [PATCH] removed meltfrac as a variable and calculate melt entirely through volume (currently untested) --- python/ctem/ctem/driver.py | 2 +- src/globals/module_globals.f90 | 6 --- src/init/init_regolith_stack.f90 | 4 -- src/io/io_read_regotrack.f90 | 39 ++-------------- src/io/io_write_regotrack.f90 | 33 ++----------- src/regolith/regolith_interior.f90 | 7 --- src/regolith/regolith_melt_glass.f90 | 9 ---- src/regolith/regolith_mix.f90 | 7 --- src/regolith/regolith_streamtube.f90 | 17 +------ src/regolith/regolith_streamtube_head.f90 | 26 ++++++++--- src/regolith/regolith_streamtube_lineseg.f90 | 46 +++++++++++++------ src/regolith/regolith_subpixel_streamtube.f90 | 31 ++++++++++--- src/regolith/regolith_transport.f90 | 2 +- src/regolith/regolith_traverse_streamtube.f90 | 10 ++-- src/util/util_init_array.f90 | 4 -- src/util/util_traverse_pop_array.f90 | 6 --- 16 files changed, 93 insertions(+), 156 deletions(-) diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 32c17cb3..0c9bd611 100644 --- a/python/ctem/ctem/driver.py +++ b/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') diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index f1851922..3dda6811 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -66,7 +66,6 @@ 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. @@ -74,8 +73,6 @@ module module_globals 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 diff --git a/src/init/init_regolith_stack.f90 b/src/init/init_regolith_stack.f90 index c9a5dd60..27649e33 100644 --- a/src/init/init_regolith_stack.f90 +++ b/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 diff --git a/src/io/io_read_regotrack.f90 b/src/io/io_read_regotrack.f90 index 9a08df9d..fc199eaf 100644 --- a/src/io/io_read_regotrack.f90 +++ b/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,23 +79,11 @@ 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 @@ -106,11 +91,6 @@ subroutine io_read_regotrack(user,surf,domain) 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,8 +107,7 @@ 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(:) @@ -137,9 +115,6 @@ subroutine io_read_regotrack(user,surf,domain) 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,22 +131,19 @@ 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 @@ -179,10 +151,7 @@ subroutine io_read_regotrack(user,surf,domain) close(FREGO) close(FCOMP) close(FAGE) - close(FDF) close(FEJM) - close(FEJMF) close(FMD) - close(FMF) return end subroutine io_read_regotrack diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90 index 2cc183ff..4e11d099 100644 --- a/src/io/io_write_regotrack.f90 +++ b/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') diff --git a/src/regolith/regolith_interior.f90 b/src/regolith/regolith_interior.f90 index 217fb1ae..ef4046ad 100644 --- a/src/regolith/regolith_interior.f90 +++ b/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 diff --git a/src/regolith/regolith_melt_glass.f90 b/src/regolith/regolith_melt_glass.f90 index b5280ac6..df47d1c4 100644 --- a/src/regolith/regolith_melt_glass.f90 +++ b/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 diff --git a/src/regolith/regolith_mix.f90 b/src/regolith/regolith_mix.f90 index 7dec2759..38aed0a1 100644 --- a/src/regolith/regolith_mix.f90 +++ b/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) diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index 385dad24..1fe3584c 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -270,15 +270,6 @@ 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 + meltinejecta - !newlayer%meltfrac = newlayer%meltvolume / totvol - !distvol(:) = distvol(:) + (newlayer%ejm*newlayer%meltdist(:)) - !newlayer%distvol(:) = newlayer%distvol(:) + distvol(:) - !newlayer%meltdist(:) = newlayer%distvol(:) / totvol - ! if (newlayer%meltfrac > 1.0_DP) then - ! write(*,*) "Melt fraction >1! (Subpixel)", xpi,ypi,crater%timestamp,crater%fcrat,crater%xlpx,crater%ylpx,& - ! newlayer%meltvolume, newlayer%totvolume, newlayer%ejm, newlayer%ejmf, totvol - ! end if else rbody = sqrt(xints(2)**2 + yints(2)**2) @@ -338,23 +329,18 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, !totvol = newlayer%totvolume - meltinejecta if (newlayer%ejm > newlayer%totvolume) then !entire pixel is ejected melt newlayer%ejm = newlayer%totvolume - newlayer%meltfrac = 1.0_DP newlayer%meltvolume = newlayer%ejm - newlayer%ejmf = 1.0_DP else if (meltinejecta + newlayer%ejm > newlayer%totvolume) then !entire pixel is melt, but not all of it is ejected meltinejecta = newlayer%totvolume - newlayer%ejm end if newlayer%meltvolume = meltinejecta + newlayer%ejm - newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume - if (newlayer%meltfrac > 1.0_DP) then !edge case caused by floating point math could result in melt fraction slightly higher than 1 - newlayer%meltfrac = 1.0_DP + if (newlayer%meltvolume > newlayer%totvolume) then !edge case caused by floating point math could result in melt fraction slightly higher than 1 factor = newlayer%totvolume / newlayer%meltvolume newlayer%meltvolume = newlayer%totvolume distvol(:) = distvol(:) * factor 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 @@ -365,7 +351,6 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%distvol(:) = newlayer%distvol(:) * factor end if newlayer%meltvolume = sum(newlayer%distvol) - newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume end if diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index 03daf09b..9460260c 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -45,6 +45,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector ! melt collector real(DP) :: recyratio real(DP) :: headmeltvol + real(DP) :: ratio !current => surfi%regolayer !current = surfi%regolayer @@ -66,9 +67,13 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector totmare = totmare + vsgly * current(M)%comp recyratio = vsgly / (user%pix**2) /current(M)%thickness age_collector(:) = age_collector(:) + current(M)%age(:) * recyratio - headmeltvol = current(M)%meltfrac * vsgly! * recyratio + ratio = vsgly / current(M)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if + headmeltvol = current(M)%meltvolume * ratio meltinejecta = meltinejecta + headmeltvol - distvol(:) = distvol(:) + (current(M)%meltdist(:)*vsgly)!*recyratio) + distvol(:) = distvol(:) + (current(M)%distvol(:) * ratio) totvol = totvol + vsgly else ! head is not intersected with layers. @@ -81,9 +86,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 - headmeltvol = current(N)%meltfrac * tothead! * recyratio + ratio = (vhead * vratio) / current(N)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if + headmeltvol = current(N)%meltvolume * ratio meltinejecta = meltinejecta + headmeltvol - distvol(:) = distvol(:) + (current(N)%meltdist(:)*(vhead*vratio))!*recyratio) + distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) totvol = totvol + tothead !current => current%next !N = N - 1 @@ -95,10 +104,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 - headmeltvol = current(N)%meltfrac * (vsgly-tothead)!*recyratio + ratio = (vsgly-tothead) / current(N)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if + headmeltvol = current(N)%meltvolume * ratio meltinejecta = meltinejecta + headmeltvol - !distvol(:) = distvol(:) + (current(N)%meltdist(:)*tothead*recyratio) !should be +0 since recyratio=0 - distvol(:) = distvol(:) + (current(N)%meltdist(:)*(vsgly-tothead)) + distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) totvol = totvol + vsgly exit end if diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90 index 03ad8b1e..a834e8c4 100644 --- a/src/regolith/regolith_streamtube_lineseg.f90 +++ b/src/regolith/regolith_streamtube_lineseg.f90 @@ -43,7 +43,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad integer(I4B) :: N,M ! Melt zone - real(DP) :: recyratio, xsh, rst + real(DP) :: recyratio, xsh, rst, ratio real(DP) :: theta1, theta2, r1, r2, vol real(DP) :: linmelt ! Shock damaged zone @@ -81,23 +81,31 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad if (ri > xmints) then vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) - ebh_recyl = (1.0 - newlayer%meltfrac) * newlayer%thickness - vsh / user%pix**2 + ebh_recyl = (1.0_DP - (newlayer%meltvolume/newlayer%totvolume)) * newlayer%thickness - vsh / user%pix**2 recyratio = max(ebh_recyl,0.0_DP) / current(N-1)%thickness age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio + ratio = max(vsgly-vsh,0.0_DP) / current(N-1)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if vol = vol + sum(current(N-1)%age(:)) * recyratio - linmelt = current(N-1)%meltfrac * max(vsgly-vsh,0.0_DP)! * recyratio + linmelt = current(N-1)%meltvolume * ratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N-1)%meltdist(:)*max(vsgly-vsh,0.0_DP))!*recyratio) + distvol(:) = distvol(:) + (current(N-1)%distvol(:) * ratio) totvol = totvol + vsgly else if (ri <= xmints .and. rip1 > xmints) then vseg = regolith_streamtube_volume_func(eradi,xmints,rip1,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N-1)%thickness age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio + ratio = max(vseg-vsh,0.0_DP) / current(N-1)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if vol = vol + sum(current(N-1)%age(:)) * recyratio - linmelt = current(N-1)%meltfrac * max(vseg-vsh,0.0_DP)! * recyratio + linmelt = current(N-1)%meltvolume * ratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N-1)%meltdist(:)*max(vseg-vsh,0.0_DP))!*recyratio) + distvol(:) = distvol(:) + (current(N-1)%distvol(:) * ratio) totvol = totvol + vseg end if exit @@ -124,10 +132,14 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rstart,rend) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio + ratio = max((vseg-vsh),0.0_DP) / current(N)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if vol = vol + sum(current(N)%age(:)) * recyratio - linmelt = current(N)%meltfrac * max(vseg-vsh,0.0_DP)! * recyratio + linmelt = current(N)%meltvolume * ratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N)%meltdist(:)*max(vseg-vsh,0.0_DP))!*recyratio) + distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) totvol = totvol + vseg end if @@ -137,10 +149,14 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rend,rstart) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio + ratio = max((vseg-vsh),0.0_DP) / current(N)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if vol = vol + sum(current(N)%age(:)) * recyratio - linmelt = current(N)%meltfrac * max(vseg-vsh,0.0_DP)! * recyratio + linmelt = current(N)%meltvolume * ratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N)%meltdist(:)*max(vseg-vsh,0.0_DP))!*recyratio) + distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) totvol = totvol + vseg end if !current => current%next @@ -162,10 +178,11 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio + ratio = max((vseg-vsh),0.0_DP) / current(N)%totvolume vol = vol + sum(current(N)%age(:)) * recyratio - linmelt = current(N)%meltfrac * max(vseg-vsh,0.0_DP)! * recyratio + linmelt = current(N)%meltvolume * ratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N)%meltdist(:)*max(vseg-vsh,0.0_DP))!*recyratio) + distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) totvol = totvol + vseg end if @@ -176,9 +193,10 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vsgly = regolith_streamtube_volume_func(eradi,ri,rip1,deltar) vmare = vmare + (vsgly - totseb) * current(N)%comp totseb = vsgly - linmelt = current(N)%meltfrac * vsgly + ratio = vsgly / current(N)%totvolume + linmelt = current(N)%meltvolume * ratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N)%meltdist(:)*vsgly) + distvol(:) = distvol(:) + (current(N)%distvol(:) * ratio) totvol = totvol + vsgly exit end if diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index 84317355..ca271fac 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -80,7 +80,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer integer(I4B) :: N,M ! Stream tube's distance from the edge of a melt zone - real(DP) :: zm, recyratio, xmints1, vseg, recyratio2 + real(DP) :: zm, recyratio, xmints1, vseg, recyratio2, ratio, ratio2 ! Parameters for calculating shocked segment of stream tube real(DP) :: x_up_sh, x_low_sh, vsh, vsh2 @@ -114,9 +114,13 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer if (eradi>xmints) then vseg = regolith_streamtube_volume_func(eradi,xmints,eradi,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,0.0_DP,eradi) - recyratio = max(vseg-vsh,0.0 )/ (user%pix**2) / (surfi%regolayer(M)%thickness) - meltinejecta = surfi%regolayer(M)%meltfrac * max((vseg-vsh),0.0_DP)! * recyratio - distvol(:) = distvol(:) + (surfi%regolayer(M)%meltdist(:)*max((vseg-vsh),0.0_DP))!*recyratio) + recyratio = max(vseg-vsh,0.0_DP) / (user%pix**2) / (surfi%regolayer(M)%thickness) + ratio = max(vseg-vsh,0.0_DP) / surfi%regolayer(M)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if + meltinejecta = surfi%regolayer(M)%meltvolume * ratio + distvol(:) = distvol(:) + (surfi%regolayer(M)%distvol(:) * ratio) totvol = vseg age_collector(:) = age_collector(:) + surfi%regolayer(M)%age(:) * recyratio vol = vol + sum(age_collector(:)) @@ -193,7 +197,11 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer recyratio = max((vsgly2 -vsh),0.0_DP)/ (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio - mvr = max((vsgly2 -vsh),0.0_DP) * current(N)%meltfrac! * recyratio + ratio2 = max((vsgly2-vsh),0.0_DP) / current(N)%totvolume + if (ratio2 > 1) then + ratio2 = 1.0_DP + end if + mvr = ratio2 * current(N)%meltvolume recyratio2 = recyratio vsh2 = vsh end if @@ -203,13 +211,22 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer recyratio = max((vsgly1-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio - mvl = max((vsgly1-vsh),0.0_DP) * current(N)%meltfrac!*recyratio + ratio = max((vsgly1-vsh),0.0_DP) / current(N)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if + mvl = ratio * current(N)%meltvolume end if !current => current%next meltinejecta = meltinejecta + mvl + mvr !distvol(:) = distvol(:) + (current(N)%meltdist(:)*((vsgly1-vsh)*recyratio)+((vsgly2-vsh2)*recyratio2)) - distvol(:) = distvol(:) + (current(N)%meltdist(:)*(max((vsgly1-vsh),0.0_DP)+max((vsgly2-vsh2),0.0_DP))) + !distvol(:) = distvol(:) + (current(N)%meltdist(:)*(max((vsgly1-vsh),0.0_DP)+max((vsgly2-vsh2),0.0_DP))) + if (ratio + ratio2 < 1.0_DP) then + distvol(:) = distvol(:) + (current(N)%distvol(:)*(ratio)) + (current(N)%distvol(:)*(ratio2)) + else + distvol(:) = distvol(:) + current(N)%distvol(:) + end if totvol = totvol + vsgly1 + vsgly2 !N = N - 1 z = z + current(N-1)%thickness diff --git a/src/regolith/regolith_transport.f90 b/src/regolith/regolith_transport.f90 index e25f833e..01ca20a1 100644 --- a/src/regolith/regolith_transport.f90 +++ b/src/regolith/regolith_transport.f90 @@ -64,7 +64,7 @@ subroutine regolith_transport(user,surfi,crater,domain,ejb,ejtble,lrad,ebh,newla melt = ejb(k)%meltfrac - ((ejb(k)%meltfrac - ejb(k+1)%meltfrac) * frac) end if - newlayer%meltfrac = melt + newlayer%meltvolume = melt * newlayer%totvolume call util_push_array(surfi%regolayer,newlayer) diff --git a/src/regolith/regolith_traverse_streamtube.f90 b/src/regolith/regolith_traverse_streamtube.f90 index d62af82c..cdceb053 100644 --- a/src/regolith/regolith_traverse_streamtube.f90 +++ b/src/regolith/regolith_traverse_streamtube.f90 @@ -47,7 +47,7 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne real(DP),parameter :: b = 1.12368 ! Melt zone - real(DP) :: recyratio + real(DP) :: recyratio, ratio ! Shock zone real(DP) :: vsh @@ -91,8 +91,12 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) recyratio = max(vseg-vsh,0.0_DP) / (user%pix**2) / surfi%regolayer(N)%thickness age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio - meltinejecta = meltinejecta + surfi%regolayer(N)%meltfrac*max(vseg-vsh,0.0_DP)!*recyratio - distvol(:) = distvol(:) + (surfi%regolayer(N)%meltdist(:)*max(vseg-vsh,0.0_DP))!*recyratio) + ratio = max(vseg-vsh,0.0_DP) / surfi%regolayer(N)%totvolume + if (ratio > 1) then + ratio = 1.0_DP + end if + meltinejecta = meltinejecta + (surfi%regolayer(N)%meltvolume * ratio) + distvol(:) = distvol(:) + (surfi%regolayer(N)%distvol(:) * ratio) totvol = totvol + vseg end if diff --git a/src/util/util_init_array.f90 b/src/util/util_init_array.f90 index b441c197..812efaf8 100644 --- a/src/util/util_init_array.f90 +++ b/src/util/util_init_array.f90 @@ -51,13 +51,9 @@ subroutine util_init_array(user,regolayer,domain,initstat) initstat = .true. regolayer(1)%thickness = sqrt(VBIG) ! This generates a buffer layer that the model should never reach if the run is structured properly regolayer(1)%comp = 0.0_DP - regolayer(1)%meltfrac = 0.0_DP regolayer(1)%porosity = 0.0_DP regolayer(1)%age(:) = 0.0_SP regolayer(1)%ejm = 0.0_DP - regolayer(1)%ejmf = 0.0_DP - allocate(regolayer(1)%meltdist(1+domain%rcnum)) - regolayer(1)%meltdist(:) = 0.0_SP allocate(regolayer(1)%distvol(1+domain%rcnum)) regolayer(1)%distvol(:) = 0.0_SP regolayer(1)%meltvolume = 0.0_DP diff --git a/src/util/util_traverse_pop_array.f90 b/src/util/util_traverse_pop_array.f90 index fc670800..e6e650f4 100644 --- a/src/util/util_traverse_pop_array.f90 +++ b/src/util/util_traverse_pop_array.f90 @@ -79,13 +79,7 @@ subroutine util_traverse_pop_array(user,regolayer,traverse_depth,poppedarray) regolayer(maxi)%ejm = (regolayer(maxi)%thickness/oldregodata(1)%thickness) * oldregodata(1)%ejm poppedarray(1)%totvolume = poppedarray(1)%thickness * user%pix * user%pix - 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)%meltfrac = regolayer(maxi)%meltvolume / regolayer(maxi)%totvolume - regolayer(maxi)%ejmf = regolayer(maxi)%ejm / regolayer(maxi)%totvolume deallocate(oldregodata)