Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix shortwave underflows and mechred divide by zero #127

Merged
merged 1 commit into from
Feb 6, 2018
Merged

fix shortwave underflows and mechred divide by zero #127

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 9 additions & 8 deletions columnphysics/icepack_mechred.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
21 changes: 11 additions & 10 deletions columnphysics/icepack_shortwave.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

!=======================================================================

Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down