Skip to content

Commit

Permalink
Fixed a rare logic error in the f_ground computation
Browse files Browse the repository at this point in the history
When computing f_ground for cells or quadrants, the number of floating corners (nfloat)
is based on how many of the four corners have f_flotation > 0.  Arbitrarily, corners with
f_flotation = 0.0 are considered to be grounded rather than floating.

This creates a potential problem for a cell where 1 <= nfloat <= 3, but all the corners
that are nominally grounded have f_flotation = 0.  (This can happen in rare cases, because
of another piece of logic that caps f_flotation at values slightly greater than or less than 0,
so that small positive values can exactly cancel with small negative values.)
Then we have a log(0) in the integral formula, leading to f_ground = NaN.

The one-line fix is to treat any cell or quadrant with minval(f_flotation) = 0 as if nfloat = 4,
implying f_ground = 0.

This fix is answer-changing in a good way; it gives finite answers in these rare cases, instead of NaNs.

The fix was introduced partway through the ISMIP6 Antarctic production runs, so that these runs could finish.
  • Loading branch information
whlipscomb committed Dec 24, 2019
1 parent c16ef4f commit 132a4c3
Showing 1 changed file with 5 additions and 1 deletion.
6 changes: 5 additions & 1 deletion libglissade/glissade_grounding_line.F90
Original file line number Diff line number Diff line change
Expand Up @@ -971,12 +971,16 @@ subroutine compute_grounded_fraction(i, j, q, rank, &

! Given nfloat, compute f_ground for each vertex
! First the easy cases...
! Note: Cells with f_flotation = 0 are arbitrarily identified as grounded rather than floating.
! This can be problematic if 1<= nfloat <= 3, but all the corners that are identified as grounded
! have f_flotation = 0 (i.e., minval(f_flotation) = 0). Then we should set f_ground = 0 (as if we had nfloat = 4).

if (nfloat == 0) then

f_ground = 1.0d0 ! fully grounded

elseif (nfloat == 4) then
!! elseif (nfloat == 4) then ! need the extra logic in the line below
elseif (nfloat == 4 .or. minval(f_flotation) == 0.0d0) then

f_ground = 0.0d0 ! fully floating

Expand Down

0 comments on commit 132a4c3

Please sign in to comment.