Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
attempt to fix some bugs in correction
  • Loading branch information
Austin Blevins committed Mar 2, 2023
1 parent c43b005 commit 9d3c1db
Showing 1 changed file with 11 additions and 4 deletions.
15 changes: 11 additions & 4 deletions src/regolith/regolith_streamtube.f90
Original file line number Diff line number Diff line change
Expand Up @@ -339,23 +339,30 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
!Apply a correction factor to ensure conservation of volume


factor = (totvol - newlayer%ejm) / newlayer%totvolume
factor = newlayer%totvolume / (totvol - newlayer%ejm)
meltinejecta = meltinejecta * factor
distvol(:) = distvol(:) * factor
totvol = newlayer%totvolume - meltinejecta
if (newlayer%ejm > newlayer%totvolume) then !entire pixel is ejected melt
newlayer%ejm = newlayer%totvolume
newlayer%meltdist(:) = newlayer%distvol(:) * 1.0
newlayer%meltfrac = 1.0
newlayer%meltvolume = newlayer%ejm
else
meltinejecta = meltinejecta * factor
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%meltfrac = (meltinejecta + newlayer%ejm) / newlayer%totvolume
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
factor = newlayer%totvolume / newlayer%meltvolume
newlayer%meltvolume = newlayer%totvolume
distvol(:) = distvol(:) * factor
end if
newlayer%distvol(:) = factor * (newlayer%distvol(:) + distvol(:))
newlayer%distvol(:) = newlayer%distvol(:) + distvol(:)
newlayer%meltdist(:) = newlayer%distvol(:) / newlayer%totvolume
newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume
end if

call util_push_array(surf(xpi,ypi)%regolayer,newlayer)
Expand Down

0 comments on commit 9d3c1db

Please sign in to comment.