fix shortwave underflows and mechred divide by zero

columnphysics/icepack_mechred.F90

@@ -814,7 +814,7 @@ subroutine ridge_itd (ncat,        aice0,           &
          Gsum        ! Gsum(n) = sum of areas in categories 0 to n
 
       real (kind=dbl_kind) :: &
-         work        ! temporary work array
+         work        ! temporary work variable
 
       real (kind=dbl_kind) :: &
          hi      , & ! ice thickness for each cat (m)
@@ -827,9 +827,9 @@ subroutine ridge_itd (ncat,        aice0,           &
       ! Initialize
       !-----------------------------------------------------------------
 
-      Gsum   (-1) = c0        ! by definition
-!      Gsum   (0)  = c1         ! to avoid divzero below
+      Gsum   (-1:ncat) = c0  ! initialize
 
+      Gsum   (-1) = c0     ! by definition
       if (aice0 > puny) then
          Gsum(0) = aice0
       else
@@ -838,7 +838,6 @@ subroutine ridge_itd (ncat,        aice0,           &
       apartic(0)  = c0
 
       do n = 1, ncat
-         Gsum   (n) = c1    ! to avoid divzero below
          apartic(n) = c0
          hrmin  (n) = c0
          hrmax  (n) = c0
@@ -864,10 +863,12 @@ subroutine ridge_itd (ncat,        aice0,           &
 
       ! normalize
 
-      work = c1 / Gsum(ncat)
-      do n = 0, ncat
-         Gsum(n) = Gsum(n) * work
-      enddo
+      if (Gsum(ncat) > c0) then
+         work = c1 / Gsum(ncat)
+         do n = 0, ncat
+            Gsum(n) = Gsum(n) * work
+         enddo
+      endif
 
       !-----------------------------------------------------------------
       ! Compute the participation function apartic; this is analogous to

columnphysics/icepack_shortwave.F90

@@ -76,8 +76,8 @@ module icepack_shortwave
          hpmin  = 0.005_dbl_kind, & ! minimum allowed melt pond depth (m)
          hp0    = 0.200_dbl_kind    ! pond depth below which transition to bare ice
 
-      real (kind=dbl_kind) :: &
-         exp_min                    ! minimum exponential value
+      real (kind=dbl_kind), parameter :: & 
+         exp_argmax = c10    ! maximum argument of exponential
 
 !=======================================================================
 
@@ -907,18 +907,13 @@ subroutine run_dEdd(dt,       tr_aero,   &
       logical (kind=log_kind) :: &
          linitonly       ! local initonly value
 
-      real (kind=dbl_kind), parameter :: & 
-         argmax = c10    ! maximum argument of exponential
-
       character(len=*),parameter :: subname='(run_dEdd)'
 
       linitonly = .false.
       if (present(initonly)) then
          linitonly = initonly
       endif
 
-      exp_min = exp(-argmax)
-
       ! cosine of the zenith angle
       call compute_coszen (tlat,          tlon, &
                            calendar_type, days_per_year, &
@@ -3223,6 +3218,9 @@ subroutine solution_dEdd                                 &
          smr      , & ! accumulator for rdif gaussian integration
          smt          ! accumulator for tdif gaussian integration
 
+      real (kind=dbl_kind) :: &
+         exp_min                    ! minimum exponential value
+
       character(len=*),parameter :: subname='(solution_dEdd)'
 
       ! Delta-Eddington solution expressions
@@ -3311,7 +3309,8 @@ subroutine solution_dEdd                                 &
             ! non-refracted beam instead
             if( srftyp < 2 .and. k < kfrsnl ) mu0n = mu0
 
-            extins = max(exp_min, exp(-lm*ts))
+            exp_min = min(exp_argmax,lm*ts)
+            extins = exp(-exp_min)
             ne = n(ue,extins)
 
             ! first calculation of rdif, tdif using Delta-Eddington formulas
@@ -3320,7 +3319,8 @@ subroutine solution_dEdd                                 &
             tdif_a(k) = c4*ue/ne
 
             ! evaluate rdir,tdir for direct beam
-            trnlay(k) = max(exp_min, exp(-ts/mu0n))
+            exp_min = min(exp_argmax,ts/mu0n)
+            trnlay(k) = exp(-exp_min)
             alp = alpha(ws,mu0n,gs,lm)
             gam = agamm(ws,mu0n,gs,lm)
             apg = alp + gam
@@ -3341,7 +3341,8 @@ subroutine solution_dEdd                                 &
                mu  = gauspt(ng)
                gwt = gauswt(ng)
                swt = swt + mu*gwt
-               trn = max(exp_min, exp(-ts/mu))
+               exp_min = min(exp_argmax,ts/mu)
+               trn = exp(-exp_min)
                alp = alpha(ws,mu,gs,lm)
                gam = agamm(ws,mu,gs,lm)
                apg = alp + gam