Skip to content

Commit

Permalink
Remove reciprocal approach that had been causing issues
Browse files Browse the repository at this point in the history
  • Loading branch information
@jonbob
jonbob committed Feb 22, 2023
1 parent 2b5e2b6 commit a578a04
Showing 1 changed file with 37 additions and 44 deletions.
81 changes: 37 additions & 44 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 @@ -4663,9 +4663,6 @@ subroutine compute_barycenter_coordinates(&
mean1, center1, xGrad1, yGrad1, &
mean2, center2, xGrad2, yGrad2)

use seaice_constants, only: &
seaicePuny

type(geometric_avg_cell_type), pointer :: &
geomAvgCell !< Input: derived type holding geometric averages

Expand All @@ -4692,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 @@ -4709,18 +4706,18 @@ subroutine compute_barycenter_coordinates(&
c0 = center0
cx = xGrad0
cy = yGrad0
if (abs(mean0) > seaicePuny) then
reciprocal = 1.0_RKIND / mean0
else
reciprocal = 0.0_RKIND

xBarycenter = 0.0_RKIND
yBarycenter = 0.0_RKIND
if (abs(mean0) > 0.0_RKIND) then
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 @@ -4733,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) > seaicePuny) then
reciprocal = 1.0_RKIND / massTracerProd
else
reciprocal = 0.0_RKIND
if (abs(massTracerProd) > 0.0_RKIND) then
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 @@ -4765,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) > seaicePuny) then
reciprocal = 1.0_RKIND / massTracerProd
else
reciprocal = 0.0_RKIND
if (abs(massTracerProd) > 0.0_RKIND) then
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 a578a04

Please sign in to comment.