Skip to content

Commit

Permalink
Minor change to prevent a rare divzero
Browse files Browse the repository at this point in the history
In subroutine glissade_vertical_remap, there are checks that the total column thickness
('thck' and single-layer thickness 'hlyr' are greater than 0.
In rare cases (encountered during a 2-km Antarctic spin-up), we can have thck very slightly > 0,
but still get a divzero error when dividing by hlyr.

The fix is to require thck and hlyr to be greater than tiny(0.0d0).
This avoids the divzero for the run that crashed.
  • Loading branch information
whlipscomb committed Jan 8, 2020
1 parent beec486 commit d245c72
Showing 1 changed file with 8 additions and 5 deletions.
13 changes: 8 additions & 5 deletions libglissade/glissade_transport.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1974,9 +1974,10 @@ subroutine glissade_vertical_remap(nx, ny, &

!-----------------------------------------------------------------
! If thck > 0, do vertical remapping of tracers
! WHL - Use tiny instead of 0.0d0 to avoid a rare divzero error.
!-----------------------------------------------------------------

if (thck > 0.d0) then
if (thck > tiny(0.0d0)) then

!-----------------------------------------------------------------
! Determine vertical coordinate z1, given input layer thicknesses.
Expand Down Expand Up @@ -2127,11 +2128,13 @@ subroutine glissade_vertical_remap(nx, ny, &
enddo

! compute new tracer values
! Note: Since thck > 0, we should have hlyr > 0 for all k.
! But to be safe, allow for thck very slightly > 0 (e.g., 1.e-300) and hlyr = 0.0.

! WHL - Use tiny instead of 0.0d0 to avoid a rare divzero error.
! Such an error occurred during a 2-km Antarctic spin-up in Dec. 2019.
! Note: Since thck > tiny, we should have hlyr > tiny for all k.
! To be safe, however, allow for the possibility that thck > tiny but hlyr is not.

do k = 1, nlyr
if (hlyr(i,j,k) > 0.0d0) then
if (hlyr(i,j,k) > tiny(0.0d0)) then
trcr(i,j,:,k) = htsum(:,k) / hlyr(i,j,k)
else
trcr(i,j,:,k) = 0.0d0
Expand Down

0 comments on commit d245c72

Please sign in to comment.