From 367d9ff8572d18ab5fd9ef81283c58a073589e47 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Mon, 27 Feb 2023 14:23:25 -0500 Subject: [PATCH] Fixed a few bugs in volume calculation --- src/regolith/regolith_streamtube_head.f90 | 15 +++++++------- src/regolith/regolith_streamtube_lineseg.f90 | 20 +++++++++---------- src/regolith/regolith_subpixel_streamtube.f90 | 17 ++++++++-------- src/regolith/regolith_traverse_streamtube.f90 | 4 ++-- 4 files changed, 29 insertions(+), 27 deletions(-) diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index 57e57c4e..3ed34b93 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -68,7 +68,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector age_collector(:) = age_collector(:) + current(M)%age(:) * recyratio headmeltvol = current(M)%meltfrac * vsgly * recyratio meltinejecta = meltinejecta + headmeltvol - distvol(:) = distvol(:) + (current(M)%meltdist(:)*vsgly*recyratio) + distvol(:) = distvol(:) + (current(M)%meltdist(:)*vsgly)!*recyratio) totvol = totvol + vsgly else ! head is not intersected with layers. @@ -81,10 +81,10 @@ 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 + headmeltvol = current(N)%meltfrac * tothead! * recyratio meltinejecta = meltinejecta + headmeltvol - distvol(:) = distvol(:) + (current(N)%meltdist(:)*tothead*recyratio) - totvol = totvol + tothead + distvol(:) = distvol(:) + (current(N)%meltdist(:)*(vhead*vratio))!*recyratio) + !totvol = totvol + tothead !current => current%next !N = N - 1 z = z + current(N-1)%thickness @@ -95,10 +95,11 @@ 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 * tothead*recyratio + headmeltvol = current(N)%meltfrac * (vsgly-tothead)!*recyratio meltinejecta = meltinejecta + headmeltvol - distvol(:) = distvol(:) + (current(N)%meltdist(:)*tothead*recyratio) - totvol = totvol + tothead + !distvol(:) = distvol(:) + (current(N)%meltdist(:)*tothead*recyratio) !should be +0 since recyratio=0 + distvol(:) = distvol(:) + (current(N)%meltdist(:)*(vsgly-tothead)) + totvol = totvol + vsgly exit end if end do diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90 index 2c690378..5dc1cb92 100644 --- a/src/regolith/regolith_streamtube_lineseg.f90 +++ b/src/regolith/regolith_streamtube_lineseg.f90 @@ -84,9 +84,9 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max(ebh_recyl,0.0_DP) / current(N-1)%thickness age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio vol = vol + sum(current(N-1)%age(:)) * recyratio - linmelt = current(N-1)%meltfrac * vsgly * recyratio + linmelt = current(N-1)%meltfrac * (vsgly-vsh)! * recyratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N-1)%meltdist(:)*vsgly*recyratio) + distvol(:) = distvol(:) + (current(N-1)%meltdist(:)*(vsgly-vsh))!*recyratio) totvol = totvol + vsgly else if (ri <= xmints .and. rip1 > xmints) then vseg = regolith_streamtube_volume_func(eradi,xmints,rip1,deltar) @@ -94,9 +94,9 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N-1)%thickness age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio vol = vol + sum(current(N-1)%age(:)) * recyratio - linmelt = current(N-1)%meltfrac * vseg * recyratio + linmelt = current(N-1)%meltfrac * (vseg-vsh)! * recyratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N-1)%meltdist(:)*vseg*recyratio) + distvol(:) = distvol(:) + (current(N-1)%meltdist(:)*(vseg-vsh))!*recyratio) totvol = totvol + vseg end if exit @@ -124,9 +124,9 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max((vseg-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 - linmelt = current(N)%meltfrac * vseg * recyratio + linmelt = current(N)%meltfrac * (vseg-vsh)! * recyratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N)%meltdist(:)*vseg*recyratio) + distvol(:) = distvol(:) + (current(N)%meltdist(:)*(vseg-vsh))!*recyratio) totvol = totvol + vseg end if @@ -137,9 +137,9 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max((vseg-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 - linmelt = current(N)%meltfrac * vseg * recyratio + linmelt = current(N)%meltfrac * (vseg-vsh)! * recyratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N)%meltdist(:)*vseg*recyratio) + distvol(:) = distvol(:) + (current(N)%meltdist(:)*(vseg-vsh))!*recyratio) totvol = totvol + vseg end if !current => current%next @@ -162,9 +162,9 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max((vseg-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 - linmelt = current(N)%meltfrac * vseg * recyratio + linmelt = current(N)%meltfrac * vseg! * recyratio meltinejecta = meltinejecta + linmelt - distvol(:) = distvol(:) + (current(N)%meltdist(:)*vseg*recyratio) + distvol(:) = distvol(:) + (current(N)%meltdist(:)*(vseg-vsh))!*recyratio) totvol = totvol + vseg end if diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index b8e16f16..3a668407 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -83,7 +83,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer real(DP) :: zm, recyratio, xmints1, vseg, recyratio2 ! Parameters for calculating shocked segment of stream tube - real(DP) :: x_up_sh, x_low_sh, vsh + real(DP) :: x_up_sh, x_low_sh, vsh, vsh2 ! The depth that a stream tube dips zmax = rip1/4.0_DP @@ -113,9 +113,9 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer 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 * vseg * recyratio - distvol(:) = distvol(:) + (surfi%regolayer(M)%meltdist(:)*vseg*recyratio) - totvol = vseg + meltinejecta = surfi%regolayer(M)%meltfrac * (vseg-vsh) * recyratio + distvol(:) = distvol(:) + (surfi%regolayer(M)%meltdist(:)*(vseg-vsh)*recyratio) + totvol = vseg - vsh age_collector(:) = age_collector(:) + surfi%regolayer(M)%age(:) * recyratio vol = vol + sum(age_collector(:)) ! write(*,*) '1',eradi, xmints, xsfints, & @@ -191,8 +191,9 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer recyratio = max((vsgly2 -vsh),0.0)/ (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio - mvr = vsgly2 * current(N)%meltfrac * recyratio + mvr = (vsgly2-vsh) * current(N)%meltfrac! * recyratio recyratio2 = recyratio + vsh2 = vsh end if if (rlefti > xmints) then @@ -200,13 +201,13 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer recyratio = max((vsgly1-vsh),0.0) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio - mvl = vsgly1 * current(N)%meltfrac*recyratio + mvl = (vsgly1-vsh) * current(N)%meltfrac!*recyratio end if !current => current%next meltinejecta = meltinejecta + mvl + mvr - distvol(:) = distvol(:) + (current(N)%meltdist(:)*((vsgly1*recyratio)+(vsgly2*recyratio2))) - !distvol = distvol + (current(N)%meltdist(:)*(vsgly1+vsgly2)) + !distvol(:) = distvol(:) + (current(N)%meltdist(:)*((vsgly1-vsh)*recyratio)+((vsgly2-vsh2)*recyratio2)) + distvol(:) = distvol(:) + (current(N)%meltdist(:)*((vsgly1-vsh)+(vsgly2-vsh2))) totvol = totvol + vsgly1 + vsgly2 !N = N - 1 z = z + current(N-1)%thickness diff --git a/src/regolith/regolith_traverse_streamtube.f90 b/src/regolith/regolith_traverse_streamtube.f90 index f8ab86ac..3ffc37a6 100644 --- a/src/regolith/regolith_traverse_streamtube.f90 +++ b/src/regolith/regolith_traverse_streamtube.f90 @@ -89,8 +89,8 @@ 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 * vseg*recyratio - distvol(:) = distvol(:) + (surfi%regolayer(N)%meltdist(:)*vseg*recyratio) + meltinejecta = meltinejecta + surfi%regolayer(N)%meltfrac*(vseg-vsh)!*recyratio + distvol(:) = distvol(:) + (surfi%regolayer(N)%meltdist(:)*(vseg-vsh))!*recyratio) totvol = totvol + vseg end if