From d245c72f3ececbf54d9314161cf86a881b13e6b3 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Wed, 8 Jan 2020 13:43:02 -0700 Subject: [PATCH] Minor change to prevent a rare divzero 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. --- libglissade/glissade_transport.F90 | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/libglissade/glissade_transport.F90 b/libglissade/glissade_transport.F90 index c3892539..0cf5ac4b 100644 --- a/libglissade/glissade_transport.F90 +++ b/libglissade/glissade_transport.F90 @@ -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. @@ -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