Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history 
added checks in age calculations to prevent arithmetic exception due to SP tiny value
  • Loading branch information
Austin Blevins committed Jun 23, 2023
1 parent fe93d4f commit 32df869
Show file tree
Hide file tree
Showing 5 changed files with 108 additions and 6 deletions.
2 changes: 1 addition & 1 deletion src/regolith/module_regolith.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer
use module_globals
implicit none
type(usertype),intent(in) :: user
type(surftype),intent(in) :: surfi
type(surftype),intent(inout) :: surfi
real(DP),intent(in) :: deltar,ri,rip1,eradi
type(regodatatype),intent(inout) :: newlayer
real(DP),intent(inout) :: meltinejecta,totvol
Expand Down
27 changes: 26 additions & 1 deletion src/regolith/regolith_streamtube_head.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,12 +40,13 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
! head intersected with underlying layers, about 30% of volume difference.
real(DP) :: z,zstart,zend,zmin,zmax
real(DP) :: tothead,totmarehead,marehead,vhead,vsgly
integer(I4B) :: N,M
integer(I4B) :: N,M,i

! melt collector
real(DP) :: recyratio
real(DP) :: headmeltvol
real(DP) :: ratio
real(SP) :: limit

!current => surfi%regolayer
!current = surfi%regolayer
Expand All @@ -71,6 +72,14 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
if (ratio > 1) then
ratio = 1.0_DP
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(M)%age(i)<limit) then
current(M)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (current(M)%age(:) * ratio)
headmeltvol = current(M)%meltvolume * ratio
meltinejecta = meltinejecta + headmeltvol
Expand All @@ -90,6 +99,14 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
if (ratio > 1) then
ratio = 1.0_DP
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(N)%age(i)<limit) then
current(N)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio)
headmeltvol = current(N)%meltvolume * ratio
meltinejecta = meltinejecta + headmeltvol
Expand All @@ -108,6 +125,14 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector
if (ratio > 1) then
ratio = 1.0_DP
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(N)%age(i)<limit) then
current(N)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio)
headmeltvol = current(N)%meltvolume * ratio
meltinejecta = meltinejecta + headmeltvol
Expand Down
43 changes: 42 additions & 1 deletion src/regolith/regolith_streamtube_lineseg.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
type(regodatatype),dimension(:),allocatable :: current
real(DP) :: z,zstart,zend,rstart,rend,r
real(DP) :: vsgly,x,vseg, vsh
integer(I4B) :: N,M
integer(I4B) :: N,M,i
real(SP) :: limit

! Melt zone
real(DP) :: recyratio, xsh, rst, ratio
Expand Down Expand Up @@ -88,6 +89,14 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
if (ratio > 1) then
ratio = 1.0_DP
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(N-1)%age(i)<limit) then
current(N-1)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (current(N-1)%age(:) * ratio)
!vol = vol + sum(current(N-1)%age(:)) * recyratio
linmelt = current(N-1)%meltvolume * ratio
Expand Down Expand Up @@ -132,6 +141,14 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
vseg = regolith_streamtube_volume_func(eradi,max(xmints,rstart),rend,deltar)
vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rstart,rend)
recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(N)%age(i)<limit) then
current(N)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + current(N)%age(:) * ratio
ratio = max((vseg-vsh),0.0_DP) / current(N)%totvolume
if (ratio > 1) then
Expand All @@ -153,6 +170,14 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
if (ratio > 1) then
ratio = 1.0_DP
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(N)%age(i)<limit) then
current(N)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio)
!vol = vol + sum(current(N)%age(:)) * recyratio
linmelt = current(N)%meltvolume * ratio
Expand Down Expand Up @@ -182,6 +207,14 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
if (ratio > 1) then
ratio = 1.0_DP
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(N)%age(i)<limit) then
current(N)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio)
!vol = vol + sum(current(N)%age(:)) * recyratio
linmelt = current(N)%meltvolume * ratio
Expand All @@ -201,6 +234,14 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad
if (ratio > 1) then
ratio = 1.0_DP
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(N)%age(i)<limit) then
current(N)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (current(N)%age(:) * ratio)
linmelt = current(N)%meltvolume * ratio
meltinejecta = meltinejecta + linmelt
Expand Down
31 changes: 29 additions & 2 deletions src/regolith/regolith_subpixel_streamtube.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer

! Arguments
type(usertype),intent(in) :: user
type(surftype),intent(in) :: surfi
type(surftype),intent(inout) :: surfi
real(DP),intent(in) :: deltar,ri,rip1,eradi
type(regodatatype),intent(inout) :: newlayer
real(DP),intent(out) :: meltinejecta, totvol
Expand All @@ -77,7 +77,8 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer
!type(regolisttype),pointer :: current
type(regodatatype),dimension(:),allocatable :: current
real(DP) :: z,zmax,zstart,zend,rlefti,rleftf,rrighti,rrightf,rc,vsgly,vsgly1,vsgly2,x,mvl,mvr
integer(I4B) :: N,M
integer(I4B) :: N,M,i
real(SP) :: limit, limit2

! Stream tube's distance from the edge of a melt zone
real(DP) :: zm, recyratio, xmints1, vseg, recyratio2, ratio, ratio2
Expand Down Expand Up @@ -108,6 +109,8 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer
vsh2 = 0.0_DP
ratio = 0.0_DP
ratio2 = 0.0_DP
limit = 0.0_SP
limit2 = 0.0_SP

! Two cases: subpixel is inside the first layer, and its volume is simply the landing ejecta blanket.
if (zend>=zmax) then
Expand All @@ -124,6 +127,14 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer
meltinejecta = surfi%regolayer(M)%meltvolume * ratio
distvol(:) = distvol(:) + (surfi%regolayer(M)%distvol(:) * ratio)
totvol = vseg
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(surfi%regolayer(M)%age(i)<limit) then
surfi%regolayer(M)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (surfi%regolayer(M)%age(:) * ratio)
!vol = vol + sum(age_collector(:))
! write(*,*) '1',eradi, xmints, xsfints, &
Expand Down Expand Up @@ -224,6 +235,22 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer
!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))
if (ratio2 > 0) then
limit2 = (TINY(1._SP) / ratio2) * 2
do i=1,size(age_collector)
if(current(N)%age(i)<limit2) then
current(N)%age(i) = 0
end if
end do
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(current(N)%age(i)<limit) then
current(N)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (current(N)%age(:)*(ratio)) + (current(N)%age*(ratio2))
else
distvol(:) = distvol(:) + current(N)%distvol(:)
Expand Down
11 changes: 10 additions & 1 deletion src/regolith/regolith_traverse_streamtube.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,8 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne
! Traversing a linked list
real(DP) :: zri,zrip1,cosi,coso,rzmax
real(DP) :: erad,z,zmin,zmax,thetast,vseg
integer(I4B) :: N
integer(I4B) :: N,i
real(SP) :: limit

real(DP),parameter :: a = 0.936457
real(DP),parameter :: b = 1.12368
Expand Down Expand Up @@ -95,6 +96,14 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne
if (ratio > 1) then
ratio = 1.0_DP
end if
if (ratio > 0) then
limit = (TINY(1._SP) / ratio) * 2
do i=1,size(age_collector)
if(surfi%regolayer(N)%age(i)<limit) then
surfi%regolayer(N)%age(i) = 0
end if
end do
end if
age_collector(:) = age_collector(:) + (surfi%regolayer(N)%age(:) * ratio)
meltinejecta = meltinejecta + (surfi%regolayer(N)%meltvolume * ratio)
distvol(:) = distvol(:) + (surfi%regolayer(N)%distvol(:) * ratio)
Expand Down

0 comments on commit 32df869

Please sign in to comment.