From 29e4b76fe0a2a5ab5802de9e00081dd189362e55 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 13 Jul 2022 16:56:34 -0400 Subject: [PATCH] Fixed underflow bug --- src/regolith/regolith_depth_model.f90 | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/regolith/regolith_depth_model.f90 b/src/regolith/regolith_depth_model.f90 index 8c98bc80..8f37eff8 100644 --- a/src/regolith/regolith_depth_model.f90 +++ b/src/regolith/regolith_depth_model.f90 @@ -19,6 +19,7 @@ subroutine regolith_depth_model(user,domain,finterval,nflux,p) use module_globals use module_util + use, intrinsic :: ieee_arithmetic use module_regolith, EXCEPT_THIS_ONE => regolith_depth_model implicit none @@ -35,7 +36,9 @@ subroutine regolith_depth_model(user,domain,finterval,nflux,p) integer(I4B) :: i, j real(DP) :: h, dr, f, fmin, fmax real(DP) :: ntotsubcrat + logical :: underflow + if (ieee_support_underflow_control(psum)) call ieee_set_underflow_mode(gradual=.true.) ! Smallest crater size in sub-pixel crater regime (regolith scaling column in "nflux") ! Do it in terms of number of craters ntotsubcrat = 0._DP @@ -93,7 +96,11 @@ subroutine regolith_depth_model(user,domain,finterval,nflux,p) f = 0.5_DP * dr * ( fmin + fmax ) psum = psum + f end do - p(2,i) = 1.0_DP - exp(-1.0_DP * PI * psum * t) + if (abs(psum) > epsilon(psum)) then + p(2,i) = 1.0_DP - exp(-1.0_DP * PI * psum * t) + else + p(2,i) = 1.0_DP + end if end do return