initial shortwave module changes

components/mpas-seaice/src/column/ice_shortwave.F90

@@ -60,6 +60,9 @@ module ice_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), parameter :: &
+         exp_argmax = c10    ! maximum argument of exponential
+
       real (kind=dbl_kind) :: &
          exp_min                    ! minimum exponential value
 
@@ -926,16 +929,11 @@ 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
-
       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, &
@@ -2397,13 +2395,13 @@ subroutine compute_dEdd (nilyr,    nslyr,    klev,  klevp,  &
 
                  ! get grain size index:
                  ! works for 25 < snw_rds < 1625 um:
-                 if (tmp_gs < 125) then
-                   tmp1 = tmp_gs/50
+                 if (tmp_gs < 125.0_dbl_kind) then
+                   tmp1 = tmp_gs/50.0_dbl_kind
                    k_bcini(k) = nint(tmp1)
-                 elseif (tmp_gs < 175) then
+                 elseif (tmp_gs < 175.0_dbl_kind) then
                    k_bcini(k) = 2
                  else
-                   tmp1 = (tmp_gs/250)+2
+                   tmp1 = (tmp_gs/250.0_dbl_kind)+c2
                    k_bcini(k) = nint(tmp1)
                  endif
               else                  ! use the largest snow grain size for ice
@@ -3347,7 +3345,10 @@ subroutine solution_dEdd                                 &
          tdr      , & ! tdir for gaussian integration
          smr      , & ! accumulator for rdif gaussian integration
          smt          ! accumulator for tdif gaussian integration
-
+
+      real (kind=dbl_kind) :: &
+         exp_min                    ! minimum exponential value
+
       ! Delta-Eddington solution expressions
       alpha(w,uu,gg,e) = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
       agamm(w,uu,gg,e) = p5*w*((c1 + c3*gg*(c1-w)*uu*uu)/(c1-e*e*uu*uu))
@@ -3385,7 +3386,7 @@ subroutine solution_dEdd                                 &
       ! value below the fresnel level, i.e. the cosine solar zenith 
       ! angle below the fresnel level for the refracted solar beam:
       mu0nij = sqrt(c1-((c1-mu0**2)/(refindx*refindx)))
- 
+
       ! compute level of fresnel refraction
       ! if ponded sea ice, fresnel level is the top of the pond.
       kfrsnl = 0
@@ -3434,7 +3435,9 @@ subroutine solution_dEdd                                 &
             ! non-refracted beam instead
             if( srftyp < 2 .and. k < kfrsnl ) mu0n = mu0
 
-            extins = max(exp_min, exp(-lm*ts))
+            !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
@@ -3443,7 +3446,9 @@ subroutine solution_dEdd                                 &
             tdif_a(k) = c4*ue/ne
 
             ! evaluate rdir,tdir for direct beam
-            trnlay(k) = max(exp_min, exp(-ts/mu0n))
+            !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
@@ -3464,7 +3469,9 @@ subroutine solution_dEdd                                 &
                mu  = gauspt(ng)
                gwt = gauswt(ng)
                swt = swt + mu*gwt
-               trn = max(exp_min, exp(-ts/mu))
+               !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
@@ -4649,13 +4656,13 @@ subroutine compute_dEdd_5bd (nilyr,    nslyr,    klev,  klevp,  &
 
                  ! get grain size index:
                  ! works for 25 < snw_rds < 1625 um:
-                 if (tmp_gs < 125) then
-                   tmp1 = tmp_gs/50
+                 if (tmp_gs < 125._dbl_kind) then
+                   tmp1 = tmp_gs/50._dbl_kind
                    k_bcini(k) = nint(tmp1)
-                 elseif (tmp_gs < 175) then
+                 elseif (tmp_gs < 175._dbl_kind) then
                    k_bcini(k) = 2
                  else
-                   tmp1 = (tmp_gs/250)+2
+                   tmp1 = (tmp_gs/250._dbl_kind) + c2
                    k_bcini(k) = nint(tmp1)
                  endif
               else                  ! use the largest snow grain size for ice