From 1c768a5ad1f26ae5a96f5272da6ecc05e9f3fd2c Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Thu, 23 Feb 2023 12:35:46 -0500 Subject: [PATCH] trying to convert output to volume; gets error in writing --- src/ejecta/ejecta_emplace.f90 | 10 +++++----- src/globals/module_globals.f90 | 1 + src/io/io_write_regotrack.f90 | 30 ++++++++++++++-------------- src/regolith/regolith_melt_glass.f90 | 13 +++++++----- src/regolith/regolith_streamtube.f90 | 10 ++++++---- 5 files changed, 35 insertions(+), 29 deletions(-) diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index ce8436e7..ae4a986f 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -287,11 +287,11 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r end do !close(74) !!$OMP END PARALLEL DO - if(user%doregotrack .and. user%testflag) then - write(*,*) 'Ejected Melt: ', vmelt - write(*,*) 'Total Melt: ', totmelt - write(*,*) 'ejected / total melt:', vmelt/totmelt - end if + ! if(user%doregotrack .and. user%testflag) then + ! write(*,*) 'Ejected Melt: ', vmelt + ! write(*,*) 'Total Melt: ', totmelt + ! write(*,*) 'ejected / total melt:', vmelt/totmelt + ! end if vmeltsheet = totmelt - vmelt ejbmass = sum(cumulative_elchange) diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 36e72477..a0393321 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -77,6 +77,7 @@ module module_globals 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 end type regodatatype type regolisttype diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90 index 7cd17398..b9be4794 100644 --- a/src/io/io_write_regotrack.f90 +++ b/src/io/io_write_regotrack.f90 @@ -39,8 +39,8 @@ subroutine io_write_regotrack(user,surf,domain) !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 - real(SP),dimension(:,:),allocatable :: age, meltdist + real(DP),dimension(:),allocatable :: meltfrac, thickness, comp, ejm, ejmf, meltvolume + real(SP),dimension(:,:),allocatable :: age, meltdist, distvol integer(kind=8) :: recsize real(DP) :: dtmp real(SP) :: stmp @@ -54,8 +54,8 @@ subroutine io_write_regotrack(user,surf,domain) open(FCOMP,file=COMPFILE,status='replace',form='unformatted') 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(FEJMF,file=EJMFFILE,status='replace',form='unformatted') + open(FEJM,file=EJMFILE,status='replace',form='unformatted') + !open(FEJMF,file=EJMFFILE,status='replace',form='unformatted') ! First pass to get stack numbers stacks_num(:,:) = 0 @@ -78,7 +78,7 @@ subroutine io_write_regotrack(user,surf,domain) do i=1,user%gridsize !current => surf(i,j)%regolayer N = stacks_num(i,j) - allocate(meltfrac(N),thickness(N),comp(N),age(MAXAGEBINS,N),meltdist(domain%rcnum,N),ejm(N),ejmf(N)) + allocate(meltvolume(N),thickness(N),comp(N),age(MAXAGEBINS,N),distvol(domain%rcnum,N),ejm(N),ejmf(N)) allocate(current,source=surf(i,j)%regolayer) do k=1,N ! meltfrac(k) = current%regodata%meltfrac @@ -86,25 +86,25 @@ subroutine io_write_regotrack(user,surf,domain) ! comp(k) = current%regodata%comp ! age(:,k) = current%regodata%age(:) ! current => current%next - meltfrac(k) = current(k)%meltfrac + meltvolume(k) = current(k)%meltvolume thickness(k) = current(k)%thickness comp(k) = current(k)%comp age(:,k) = current(k)%age(:) !write(*,*) i, j - meltdist(:,k) = current(k)%meltdist(:) + distvol(:,k) = current(k)%distvol(:) ejm(k) = current(k)%ejm - ejmf(k) = current(k)%ejmf + !ejmf(k) = current(k)%ejmf end do deallocate(current) - write(FMELT) meltfrac(:) + write(FMELT) meltvolume(:) write(FREGO) thickness(:) write(FCOMP) comp(:) write(FAGE) age(:,:) - write(FMD) meltdist(:,:) - !write(FEJM) ejm(:) - write(FEJMF) ejmf(:) - deallocate(meltfrac,thickness,comp,age,meltdist,ejm,ejmf) + write(FMD) distvol(:,:) + write(FEJM) ejm(:) + !write(FEJMF) ejmf(:) + deallocate(meltfrac,thickness,comp,age,distvol,ejm,ejmf) end do end do close(FMELT) @@ -112,8 +112,8 @@ subroutine io_write_regotrack(user,surf,domain) close(FCOMP) close(FAGE) close(FMD) - !close(FEJM) - close(FEJMF) + close(FEJM) + !close(FEJMF) 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_melt_glass.f90 b/src/regolith/regolith_melt_glass.f90 index 5caa511a..b7944a26 100644 --- a/src/regolith/regolith_melt_glass.f90 +++ b/src/regolith/regolith_melt_glass.f90 @@ -131,7 +131,7 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad 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%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) @@ -145,7 +145,7 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad volm1 = vst - volv1 melt = volm1 newlayer%meltvolume = melt - newlayer%totvolume = volm1 + !newlayer%totvolume = volm1 newlayer%meltfrac = 1.0_DP newlayer%ejm = melt newlayer%ejmf = 1.0_DP @@ -173,8 +173,8 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad 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%totvolume = newlayer%thickness * user%pix * user%pix + newlayer%meltfrac = melt / newlayer%totvolume newlayer%ejm = melt newlayer%ejmf = newlayer%meltfrac @@ -182,8 +182,11 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad allocate(newlayer%meltdist((domain%rcnum))) newlayer%meltdist(:) = 0.0_SP + allocate(newlayer%distvol((domain%rcnum))) + newlayer%distvol(:) = 0.0_SP if(domain%currentqmc) then - newlayer%meltdist(domain%nqmc) = newlayer%meltfrac !might change this to melt + newlayer%meltdist(domain%nqmc) = newlayer%meltfrac + newlayer%distvol(domain%nqmc) = newlayer%meltvolume end if n_age = max(ceiling(age / age_resolution), 1) diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index 48ecb89c..7eea441f 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -273,11 +273,12 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%meltvolume = newlayer%meltvolume + meltinejecta newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume distvol(:) = distvol(:) + (newlayer%ejm*newlayer%meltdist(:)) + newlayer%distvol(:) = distvol(:) newlayer%meltdist(:) = distvol(:) / newlayer%totvolume - 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 + ! 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) @@ -326,6 +327,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%meltvolume = newlayer%meltvolume + meltinejecta newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume distvol(:) = distvol(:) + (newlayer%ejm*newlayer%meltdist(:)) + newlayer%distvol(:) = distvol(:) newlayer%meltdist(:) = distvol(:) / newlayer%totvolume end if if (newlayer%meltfrac > 1.0_DP) then