From 24755dda595f9ca72856b41dbd1ef27612894dbb Mon Sep 17 00:00:00 2001 From: apcraig Date: Mon, 5 Feb 2018 19:03:03 +0000 Subject: [PATCH] fix shortwave underflows and mechred divide by zero --- columnphysics/icepack_mechred.F90 | 17 +++++++++-------- columnphysics/icepack_shortwave.F90 | 21 +++++++++++---------- 2 files changed, 20 insertions(+), 18 deletions(-) diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90 index df13b0f9f..5c148814f 100644 --- a/columnphysics/icepack_mechred.F90 +++ b/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 diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 5c566cf24..8ba56d6ac 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/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,9 +907,6 @@ 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. @@ -917,8 +914,6 @@ subroutine run_dEdd(dt, tr_aero, & 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