From 2b5e2b6f5251458b4b8cea7193eaf5260a6e0bfe Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Wed, 26 Oct 2022 10:24:38 -0700 Subject: [PATCH 1/2] Change checks from 0.0 to puny in compute_barycenter_coordinates --- .../src/shared/mpas_seaice_advection_incremental_remap.F | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F b/components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F index 49c2a314e3d2..0cb4a268a271 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F @@ -4663,6 +4663,9 @@ 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 @@ -4706,7 +4709,7 @@ subroutine compute_barycenter_coordinates(& c0 = center0 cx = xGrad0 cy = yGrad0 - if (abs(mean0) > 0.0_RKIND) then + if (abs(mean0) > seaicePuny) then reciprocal = 1.0_RKIND / mean0 else reciprocal = 0.0_RKIND @@ -4731,7 +4734,7 @@ subroutine compute_barycenter_coordinates(& cyy = yGrad0 * yGrad1 massTracerProd = mean0 * mean1 - if (abs(massTracerProd) > 0.0_RKIND) then + if (abs(massTracerProd) > seaicePuny) then reciprocal = 1.0_RKIND / massTracerProd else reciprocal = 0.0_RKIND @@ -4763,7 +4766,7 @@ subroutine compute_barycenter_coordinates(& cyyy = yGrad0 * yGrad1 * yGrad2 massTracerProd = mean0 * mean1 * mean2 - if (abs(massTracerProd) > 0.0_RKIND) then + if (abs(massTracerProd) > seaicePuny) then reciprocal = 1.0_RKIND / massTracerProd else reciprocal = 0.0_RKIND From a578a04e28c472e5393e0713df25af26c1c9eedd Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Wed, 22 Feb 2023 16:56:23 -0600 Subject: [PATCH 2/2] Remove reciprocal approach that had been causing issues --- .../mpas_seaice_advection_incremental_remap.F | 81 +++++++++---------- 1 file changed, 37 insertions(+), 44 deletions(-) diff --git a/components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F b/components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F index 0cb4a268a271..d57d1ddd29b3 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_advection_incremental_remap.F @@ -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 @@ -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. @@ -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 @@ -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 @@ -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