Skip to content

Commit

Permalink
Merge branch 'jonbob/seaice/fix-ir-zero-check' (PR #5252)
Browse files Browse the repository at this point in the history 
Fix floating-point exception in MPASSI incremental remap code

High-res (ne120pg2_r0125_oRRS18to6v3.WCYCL1950) B-case tests on summit
using gnu have been failing with a floating-point exception that points
to a line in the seaice incremental remap code:
       if (abs(massTracerProd) > 0.0_RKIND) then
@whannah1 reported the issue and has found that the error occurs when
massTracerProd is extremely small, so its inverse ends up being
extremely large. He also found that it can be avoided by replacing the
0.0 with a tiny number. However, after discussion this PR fixes the
issue by removing the reciprocal calculation and dividing by the mean0
(and massTracerProd) quantities instead -- @ambrad thank you for the
suggestion.

Fixes #5463

[non-BFB]
  • Loading branch information
@jonbob
jonbob committed Sep 22, 2023
2 parents 4311d16 + a578a04 commit a746149
Showing 1 changed file with 34 additions and 38 deletions.
72 changes: 34 additions & 38 deletions components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4689,7 +4689,7 @@ subroutine compute_barycenter_coordinates(&
cxxx, cxxy, cxyy, cyyy, &
cxxxx, cxxxy, cxxyy, cxyyy, cyyyy

real(kind=RKIND) :: reciprocal, massTracerProd
real(kind=RKIND) :: massTracerProd

! compute the barycenter coordinates, depending on how many tracer types were passed in
! Note: These formulas become more complex as the tracer hierarchy deepens. See the model documentation for details.
Expand All @@ -4706,18 +4706,18 @@ subroutine compute_barycenter_coordinates(&
c0 = center0
cx = xGrad0
cy = yGrad0

xBarycenter = 0.0_RKIND
yBarycenter = 0.0_RKIND
if (abs(mean0) > 0.0_RKIND) then
reciprocal = 1.0_RKIND / mean0
else
reciprocal = 0.0_RKIND
xBarycenter = (c0 * geomAvgCell % x (iCell) &
+ cx * geomAvgCell % xx(iCell) + cy * geomAvgCell % xy(iCell)) &
/ mean0
yBarycenter = (c0 * geomAvgCell % y (iCell) &
+ cx * geomAvgCell % xy(iCell) + cy * geomAvgCell % yy(iCell)) &
/ mean0
endif

xBarycenter = (c0 * geomAvgCell % x (iCell) &
+ cx * geomAvgCell % xx(iCell) + cy * geomAvgCell % xy(iCell)) &
* reciprocal
yBarycenter = (c0 * geomAvgCell % y (iCell) &
+ cx * geomAvgCell % xy(iCell) + cy * geomAvgCell % yy(iCell)) &
* reciprocal

elseif (.not.present(mean2)) then

Expand All @@ -4730,22 +4730,20 @@ subroutine compute_barycenter_coordinates(&
cxy = xGrad0 * yGrad1 + yGrad0 * xGrad1
cyy = yGrad0 * yGrad1

xBarycenter = 0.0_RKIND
yBarycenter = 0.0_RKIND
massTracerProd = mean0 * mean1
if (abs(massTracerProd) > 0.0_RKIND) then
reciprocal = 1.0_RKIND / massTracerProd
else
reciprocal = 0.0_RKIND
xBarycenter = (c0 * geomAvgCell % x (iCell) &
+ cx * geomAvgCell % xx (iCell) + cy * geomAvgCell % xy (iCell) &
+ cxx * geomAvgCell % xxx(iCell) + cxy * geomAvgCell % xxy(iCell) + cyy * geomAvgCell % xyy(iCell)) &
/ massTracerProd
yBarycenter = (c0 * geomAvgCell % y (iCell) &
+ cx * geomAvgCell % xy (iCell) + cy * geomAvgCell % yy (iCell) &
+ cxx * geomAvgCell % xxy(iCell) + cxy * geomAvgCell % xyy(iCell) + cyy * geomAvgCell % yyy(iCell)) &
/ massTracerProd
endif

xBarycenter = (c0 * geomAvgCell % x (iCell) &
+ cx * geomAvgCell % xx (iCell) + cy * geomAvgCell % xy (iCell) &
+ cxx * geomAvgCell % xxx(iCell) + cxy * geomAvgCell % xxy(iCell) + cyy * geomAvgCell % xyy(iCell)) &
* reciprocal
yBarycenter = (c0 * geomAvgCell % y (iCell) &
+ cx * geomAvgCell % xy (iCell) + cy * geomAvgCell % yy (iCell) &
+ cxx * geomAvgCell % xxy(iCell) + cxy * geomAvgCell % xyy(iCell) + cyy * geomAvgCell % yyy(iCell)) &
* reciprocal

else ! mass field, tracer1 and tracer2 are all present

! barycenter = center of mass*tracer1*tracer2
Expand All @@ -4762,26 +4760,24 @@ subroutine compute_barycenter_coordinates(&
cxyy = yGrad0 * yGrad1 * xGrad2 + yGrad0 * xGrad1 * yGrad2 + xGrad0 * yGrad1 * yGrad2
cyyy = yGrad0 * yGrad1 * yGrad2

xBarycenter = 0.0_RKIND
yBarycenter = 0.0_RKIND
massTracerProd = mean0 * mean1 * mean2
if (abs(massTracerProd) > 0.0_RKIND) then
reciprocal = 1.0_RKIND / massTracerProd
else
reciprocal = 0.0_RKIND
xBarycenter = (c0 * geomAvgCell % x (iCell) &
+ cx * geomAvgCell % xx (iCell) + cy * geomAvgCell % xy (iCell) &
+ cxx * geomAvgCell % xxx (iCell) + cxy * geomAvgCell % xxy (iCell) + cyy * geomAvgCell % xyy (iCell) &
+ cxxx * geomAvgCell % xxxx(iCell) + cxxy * geomAvgCell % xxxy(iCell) + cxyy * geomAvgCell % xxyy(iCell) &
+ cyyy * geomAvgCell % xyyy(iCell)) &
/ massTracerProd
yBarycenter = (c0 * geomAvgCell % y (iCell) &
+ cx * geomAvgCell % xy (iCell) + cy * geomAvgCell % yy (iCell) &
+ cxx * geomAvgCell % xxy (iCell) + cxy * geomAvgCell % xyy (iCell) + cyy * geomAvgCell % yyy (iCell) &
+ cxxx * geomAvgCell % xxxy(iCell) + cxxy * geomAvgCell % xxyy(iCell) + cxyy * geomAvgCell % xyyy(iCell) &
+ cyyy * geomAvgCell % yyyy(iCell)) &
/ massTracerProd
endif

xBarycenter = (c0 * geomAvgCell % x (iCell) &
+ cx * geomAvgCell % xx (iCell) + cy * geomAvgCell % xy (iCell) &
+ cxx * geomAvgCell % xxx (iCell) + cxy * geomAvgCell % xxy (iCell) + cyy * geomAvgCell % xyy (iCell) &
+ cxxx * geomAvgCell % xxxx(iCell) + cxxy * geomAvgCell % xxxy(iCell) + cxyy * geomAvgCell % xxyy(iCell) &
+ cyyy * geomAvgCell % xyyy(iCell)) &
* reciprocal
yBarycenter = (c0 * geomAvgCell % y (iCell) &
+ cx * geomAvgCell % xy (iCell) + cy * geomAvgCell % yy (iCell) &
+ cxx * geomAvgCell % xxy (iCell) + cxy * geomAvgCell % xyy (iCell) + cyy * geomAvgCell % yyy (iCell) &
+ cxxx * geomAvgCell % xxxy(iCell) + cxxy * geomAvgCell % xxyy(iCell) + cxyy * geomAvgCell % xyyy(iCell) &
+ cyyy * geomAvgCell % yyyy(iCell)) &
* reciprocal

endif

end subroutine compute_barycenter_coordinates
Expand Down

0 comments on commit a746149

Please sign in to comment.