fix mechred exp underflows (#169)

CICE-Consortium · Feb 13, 2018 · 3e6b077 · 3e6b077
1 parent 5e964c4
commit 3e6b077
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 6 deletions.
diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90
@@ -64,6 +64,7 @@ module icepack_mechred
                 icepack_step_ridge
 
       real (kind=dbl_kind), parameter :: & 
+         exp_argmax = 100.0_dbl_kind, &    ! maximum argument of exponential for underflow
          Cs = p25         , & ! fraction of shear energy contrbtng to ridging 
          fsnowrdg = p5    , & ! snow fraction that survives in ridging 
          Gstar  = p15     , & ! max value of G(h) that participates 
@@ -398,8 +399,6 @@ subroutine ridge_ice (dt,          ndtd,       &
                            aredistn,    vredistn)    
          if (icepack_warnings_aborted(subname)) return
 
-         if (icepack_warnings_aborted(subname)) return
-
       !-----------------------------------------------------------------
       ! Make sure the new area = 1.  If not (because the closing
       ! and opening rates were reduced above), prepare to ridge again
@@ -1174,6 +1173,8 @@ subroutine ridge_shift (ntrcr,       dt,              &
 
       real (kind=dbl_kind) :: &
          work       , & ! temporary variable
+         expL_arg   , & ! temporary exp arg values
+         expR_arg   , & ! temporary exp arg values
          closing_gross  ! rate at which area removed, not counting
                         ! area of new ridges
 
@@ -1452,8 +1453,10 @@ subroutine ridge_shift (ntrcr,       dt,              &
                      else
                         hL = max (hi1, hin_max(nr-1))
                         hR = hin_max(nr)
-                        expL = exp(-(hL-hi1)/hexp)
-                        expR = exp(-(hR-hi1)/hexp)
+                        expL_arg = min(((hL-hi1)/hexp),exp_argmax)
+                        expR_arg = min(((hR-hi1)/hexp),exp_argmax)
+                        expL = exp(-(expL_arg))
+                        expR = exp(-(expR_arg))
                         farea = expL - expR
                         fvol  = ((hL + hexp)*expL  &
                                     - (hR + hexp)*expR) / (hi1 + hexp)
@@ -1465,7 +1468,8 @@ subroutine ridge_shift (ntrcr,       dt,              &
                      hexp = hrexp(n)
 
                      hL = max (hi1, hin_max(nr-1))
-                     expL = exp(-(hL-hi1)/hexp)
+                     expL_arg = min(((hL-hi1)/hexp),exp_argmax)
+                     expL = exp(-(expL_arg))
                      farea = expL
                      fvol  = (hL + hexp)*expL / (hi1 + hexp)
 

diff --git a/columnphysics/version.txt b/columnphysics/version.txt
@@ -1 +1 @@
-Icepack 0.0.01
+Icepack 1.0.0
diff --git a/configuration/scripts/tests/travis_suite.ts b/configuration/scripts/tests/travis_suite.ts
@@ -1,6 +1,8 @@
 # Test         Grid    PEs        Sets    BFB-compare
 smoke          col     1x1        diag1,run1year
 smoke          col     1x1        debug,run1year
+smoke          col     1x1        debug,bgcISPOL
+smoke          col     1x1        debug,bgcNICE
 smoke          col     1x1        debug,run1year,thermo1
 smoke          col     1x1        debug,run1year,swccsm3
 smoke          col     1x1        debug,run1year,alt01
@@ -11,6 +13,8 @@ restart        col     1x1        diag1
 restart        col     1x1        pondcesm
 restart        col     1x1        pondlvl
 restart        col     1x1        pondtopo
+restart        col     1x1        bgcISPOL
+restart        col     1x1        bgcNICE
 restart        col     1x1        thermo1
 restart        col     1x1        swccsm3
 restart        col     1x1        alt01