trying to convert output to volume; gets error in writing

MintonGroup · Feb 23, 2023 · 1c768a5 · 1c768a5
1 parent b6ec18e
commit 1c768a5
Show file tree

Hide file tree

Showing 5 changed files with 35 additions and 29 deletions.
diff --git 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
@@ -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
@@ -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,42 +78,42 @@ 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
             ! thickness(k) = current%regodata%thickness
             ! 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)
    close(FREGO)
    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
@@ -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,17 +173,20 @@ 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
 
    end if
 
    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
@@ -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