From 6afe99de2eb8bbc81fe62b927a13b8cf49e276c2 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Mon, 12 Feb 2018 14:36:09 -0800 Subject: [PATCH] fix underflow in brine (#158) --- columnphysics/icepack_brine.F90 | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/columnphysics/icepack_brine.F90 b/columnphysics/icepack_brine.F90 index 86e114a19..9b5c85d80 100644 --- a/columnphysics/icepack_brine.F90 +++ b/columnphysics/icepack_brine.F90 @@ -518,12 +518,12 @@ subroutine update_hbrine (meltt, & dhrunoff = -dhS_top*aice0 hbrocn = max(c0,hbrocn - dhrunoff) exp_arg = darcy_coeff/bphi_min*dt -! tcx tcraig avoids underflows but is not bit-for-bit -! if (exp_arg > exp_argmax) then -! hbrocn_new = c0 -! else +! tcraig avoids underflows + if (exp_arg > exp_argmax) then + hbrocn_new = c0 + else hbrocn_new = hbrocn*exp(-exp_arg) -! endif + endif hbr = max(hbrmin, h_ocn + hbrocn_new) hbrocn_new = hbr-h_ocn darcy_V = -SIGN((hbrocn-hbrocn_new)/dt*bphi_min, hbrocn) @@ -532,7 +532,13 @@ subroutine update_hbrine (meltt, & dh_top_chl = dh_top_chl - darcy_V_chl*dt/bphi_min + dhrunoff dh_direct = dhrunoff elseif (hbrocn < c0 .AND. hbr > thinS) then - hbrocn_new = hbrocn*exp(-darcy_coeff/bphi_min*dt) + exp_arg = darcy_coeff/bphi_min*dt +! tcraig avoids underflows + if (exp_arg > exp_argmax) then + hbrocn_new = c0 + else + hbrocn_new = hbrocn*exp(-exp_arg) + endif dhflood = max(c0,hbrocn_new - hbrocn)*aice0 hbr = max(hbrmin, h_ocn + hbrocn_new) darcy_V = -SIGN((hbrocn-hbrocn_new + dhflood)/dt*bphi_min, hbrocn)