From 87a3788d42dce366298a807ca40d9a5faa1ec257 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 25 Jul 2023 17:14:15 +0000 Subject: [PATCH 01/18] Modified doc to specify that l_fixed_area is T for C-grid --- doc/source/science_guide/sg_horiztrans.rst | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/doc/source/science_guide/sg_horiztrans.rst b/doc/source/science_guide/sg_horiztrans.rst index 10b668755..4ccf00e9b 100644 --- a/doc/source/science_guide/sg_horiztrans.rst +++ b/doc/source/science_guide/sg_horiztrans.rst @@ -477,9 +477,16 @@ Remote Sensing Center (Norway), who applied an earlier version of the CICE remapping scheme to an ocean model. The implementation in CICE is somewhat more general, allowing for departure regions lying on both sides of a cell edge. The extra triangle is constrained to lie in one -but not both of the grid cells that share the edge. Since this option -has yet to be fully tested in CICE, the current default is -`l\_fixed\_area` = false. +but not both of the grid cells that share the edge. + +The default value for the B grid is `l\_fixed\_area` = false. However, +idealized tests with the C grid have shown that prognostic fields such +as sea ice concentration exhibit a checkerboard pattern with +`l\_fixed\_area` = false. The logical `l\_fixed\_area` is therefore set +to true when using the C grid. The edge areas `edgearea\_e` and `edgearea\_n` +are in this case calculated with the C grid velocity components :math:`uvelE` +and :math:`vvelN`. + We made one other change in the scheme of :cite:`Dukowicz00` for locating triangles. In their paper, departure points are defined by From c66baf4636fb13cc165726dc5f18bddf0971215d Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 25 Jul 2023 18:00:37 +0000 Subject: [PATCH 02/18] Initial modifs to calc areafact based on linear interpolation of left and rigth values --- .../cicedyn/dynamics/ice_transport_remap.F90 | 267 ++++++++++++++---- 1 file changed, 210 insertions(+), 57 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index eb0dd17cf..0824e1dfa 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -1767,9 +1767,11 @@ subroutine locate_triangles (nx_block, ny_block, & real (kind=dbl_kind), dimension(nx_block,ny_block) :: & dx, dy , & ! scaled departure points - areafac_c , & ! area scale factor at center of edge areafac_l , & ! area scale factor at left corner - areafac_r ! area scale factor at right corner + areafac_r , & ! area scale factor at right corner + areafac_c , & ! (areafac_l + areafac_r)/2 + areafacS , & ! Sum of areafac_l and areafac_r + areafacD ! Difference of areafac_l and areafac_r real (kind=dbl_kind) :: & xcl, ycl , & ! coordinates of left corner point @@ -1790,6 +1792,7 @@ subroutine locate_triangles (nx_block, ny_block, & area1, area2 , & ! temporary triangle areas area3, area4 , & ! area_c , & ! center polygon area + areatp, xpm , & ! to calculate areafact puny , & ! w1, w2 ! work variables @@ -1859,9 +1862,11 @@ subroutine locate_triangles (nx_block, ny_block, & if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - areafac_c(:,:) = c0 areafac_l(:,:) = c0 areafac_r(:,:) = c0 + areafac_c(:,:) = c0 + areafacS(:,:) = c0 + areafacD(:,:) = c0 do ng = 1, ngroups do j = 1, ny_block do i = 1, nx_block @@ -1912,6 +1917,8 @@ subroutine locate_triangles (nx_block, ny_block, & areafac_l(i,j) = dxu(i-1,j)*dyu(i-1,j) areafac_r(i,j) = dxu(i ,j)*dyu(i ,j) areafac_c(i,j) = p5*(areafac_l(i,j) + areafac_r(i,j)) + areafacS (i,j) = areafac_r(i,j) + areafac_l(i,j) + areafacD (i,j) = areafac_r(i,j) - areafac_l(i,j) enddo enddo @@ -1945,7 +1952,9 @@ subroutine locate_triangles (nx_block, ny_block, & do i = ib, ie areafac_l(i,j) = dxu(i,j )*dyu(i,j ) areafac_r(i,j) = dxu(i,j-1)*dyu(i,j-1) - areafac_c(i,j) = p5 * (areafac_l(i,j) + areafac_r(i,j)) + areafac_c(i,j) = p5*(areafac_l(i,j) + areafac_r(i,j)) + areafacS (i,j) = areafac_r(i,j) + areafac_l(i,j) + areafacD (i,j) = areafac_r(i,j) - areafac_l(i,j) enddo enddo @@ -2070,7 +2079,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl - areafact(i,j,ng) = -areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_l(i,j) elseif (yil < c0 .and. xdl < xcl .and. ydl < c0) then @@ -2085,7 +2097,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yil iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl - areafact(i,j,ng) = areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_l(i,j) elseif (yil < c0 .and. xdl < xcl .and. ydl >= c0) then @@ -2100,7 +2115,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl - areafact(i,j,ng) = areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_l(i,j) ! BL1 (group 3) @@ -2113,7 +2131,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yil iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl - areafact(i,j,ng) = areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_l(i,j) elseif (yil > c0 .and. xdl < xcl .and. ydl < c0) then @@ -2128,7 +2149,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl - areafact(i,j,ng) = -areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_l(i,j) ! BL2 (group 1) @@ -2141,7 +2165,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl - areafact(i,j,ng) = -areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_l(i,j) endif ! TL and BL triangles @@ -2163,7 +2190,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yir iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr - areafact(i,j,ng) = -areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_r(i,j) elseif (yir < c0 .and. xdr >= xcr .and. ydr < c0) then @@ -2178,7 +2208,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br - areafact(i,j,ng) = areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_r(i,j) elseif (yir < c0 .and. xdr >= xcr .and. ydr >= c0) then @@ -2193,7 +2226,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr - areafact(i,j,ng) = areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_r(i,j) ! BR1 (group 3) @@ -2206,7 +2242,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br - areafact(i,j,ng) = areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_r(i,j) elseif (yir > c0 .and. xdr >= xcr .and. ydr < c0) then @@ -2221,7 +2260,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yir iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr - areafact(i,j,ng) = -areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_r(i,j) ! BR2 (group 2) @@ -2234,7 +2276,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br - areafact(i,j,ng) = -areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_r(i,j) endif ! TR and BR triangles @@ -2303,7 +2348,7 @@ subroutine locate_triangles (nx_block, ny_block, & area_c = edgearea(i,j) - area1 - area2 - area3 ! shift midpoint so that the area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) & + w1 = c2*area_c/areafac_c(i,j) & ! jlem could be improved... + (xdr-xcl)*ydl + (xcr-xdl)*ydr w2 = (xdr-xdl)**2 + (ydr-ydl)**2 w1 = w1/w2 @@ -2338,11 +2383,11 @@ subroutine locate_triangles (nx_block, ny_block, & ydm = p5 * ydr ! compute area of triangle adjacent to left corner - area4 = p5 * (xcl - xic) * ydl * areafac_l(i,j) + area4 = p5 * (xcl - xic) * ydl * areafac_l(i,j) ! jlem could be improved... area_c = edgearea(i,j) - area1 - area2 - area3 - area4 ! shift midpoint so that area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr + w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr ! jlem could be improved... w2 = (xdr-xic)**2 + ydr**2 w1 = w1/w2 xdm = xdm + ydr*w1 @@ -2366,11 +2411,11 @@ subroutine locate_triangles (nx_block, ny_block, & xdm = p5 * (xicr + xdl) ydm = p5 * ydl - area4 = p5 * (xic - xcr) * ydr * areafac_r(i,j) + area4 = p5 * (xic - xcr) * ydr * areafac_r(i,j) ! jlem could be improved... area_c = edgearea(i,j) - area1 - area2 - area3 - area4 ! shift midpoint so that area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl + w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl ! jlem could be improved... w2 = (xic-xdl)**2 + ydl**2 w1 = w1/w2 xdm = xdm - ydl*w1 @@ -2410,7 +2455,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! TC2a (group 5) @@ -2423,7 +2471,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! TC3a (group 6) @@ -2436,7 +2487,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl >= c0 .and. ydr >= c0 .and. ydm < c0) then ! rare @@ -2451,7 +2505,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! TC2b (group 5) @@ -2464,7 +2521,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! BC3b (group 6) @@ -2477,7 +2537,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr < c0 .and. ydm < c0) then @@ -2492,7 +2555,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ycr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) ! BC2a (group 5) @@ -2505,7 +2571,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) ! BC3a (group 6) @@ -2518,7 +2587,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr < c0 .and. ydm >= c0) then ! rare @@ -2533,7 +2605,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) ! BC2b (group 5) @@ -2546,7 +2621,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) ! TC3b (group 6) @@ -2559,7 +2637,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! Now consider cases where the two DPs lie in different grid cells ! For these cases, one triangle is given the area factor associated @@ -2579,7 +2660,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! BC2b (group 5) @@ -2592,7 +2676,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_r(i,j) ! TC3b (group 6) @@ -2605,7 +2692,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 & .and. ydm < c0 ) then ! less common @@ -2621,7 +2711,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! BC2b (group 5) @@ -2634,7 +2727,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_r(i,j) ! BC3b (group 6) @@ -2647,7 +2743,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 & .and. ydm < c0) then @@ -2663,7 +2762,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_l(i,j) ! BC2b (group 5) @@ -2676,7 +2778,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) ! BC3b (group 6) @@ -2689,7 +2794,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 & .and. ydm >= c0) then ! less common @@ -2705,7 +2813,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_l(i,j) ! BC2b (group 5) @@ -2718,7 +2829,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) ! TC3b (group 6) @@ -2731,7 +2845,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 & .and. ydm >= c0) then @@ -2747,7 +2864,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_l(i,j) ! TC2b (group 5) @@ -2760,7 +2880,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! TC3b (group 6) @@ -2773,7 +2896,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 & .and. ydm < c0) then ! less common @@ -2789,7 +2915,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_l(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_l(i,j) ! TC2b (group 5) @@ -2802,7 +2931,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) ! BC3b (group 6) @@ -2815,7 +2947,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 & .and. ydm < c0) then @@ -2831,7 +2966,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) ! TC2b (group 5) @@ -2844,7 +2982,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_r(i,j) ! BC3b (group 6) @@ -2857,7 +2998,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 & .and. ydm >= c0) then ! less common @@ -2873,7 +3017,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = areatp + !areafact(i,j,ng) = areafac_c(i,j) ! TC2b (group 5) @@ -2886,7 +3033,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_r(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_r(i,j) ! TC3b (group 6) @@ -2899,7 +3049,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) + xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) + areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areafact(i,j,ng) = -areatp + !areafact(i,j,ng) = -areafac_c(i,j) endif ! TC and BC triangles From 33317f3f4f5974a1f9624ddb5c4c498d7816b6d7 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 25 Jul 2023 18:49:06 +0000 Subject: [PATCH 03/18] put back l_fixed_area = .true. for C-grid --- cicecore/cicedyn/dynamics/ice_transport_remap.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 0824e1dfa..6cc0a263d 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -317,7 +317,7 @@ subroutine init_remap !------------------------------------------------------------------- if (grid_ice == 'CD' .or. grid_ice == 'C') then - l_fixed_area = .false. !jlem temporary + l_fixed_area = .true. else l_fixed_area = .false. endif From 9df6bc7eb327f5d8efa0f06a668c0d2029bd215a Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Wed, 26 Jul 2023 13:15:08 +0000 Subject: [PATCH 04/18] added temporary comments for PR review --- cicecore/cicedyn/dynamics/ice_transport_remap.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 6cc0a263d..9aad9ed8b 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -2348,7 +2348,7 @@ subroutine locate_triangles (nx_block, ny_block, & area_c = edgearea(i,j) - area1 - area2 - area3 ! shift midpoint so that the area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) & ! jlem could be improved... + w1 = c2*area_c/areafac_c(i,j) & ! jlem should this be modified? + (xdr-xcl)*ydl + (xcr-xdl)*ydr w2 = (xdr-xdl)**2 + (ydr-ydl)**2 w1 = w1/w2 @@ -2383,11 +2383,11 @@ subroutine locate_triangles (nx_block, ny_block, & ydm = p5 * ydr ! compute area of triangle adjacent to left corner - area4 = p5 * (xcl - xic) * ydl * areafac_l(i,j) ! jlem could be improved... + area4 = p5 * (xcl - xic) * ydl * areafac_l(i,j) ! jlem should this be modified? area_c = edgearea(i,j) - area1 - area2 - area3 - area4 ! shift midpoint so that area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr ! jlem could be improved... + w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr ! jlem should this be modified? w2 = (xdr-xic)**2 + ydr**2 w1 = w1/w2 xdm = xdm + ydr*w1 @@ -2411,11 +2411,11 @@ subroutine locate_triangles (nx_block, ny_block, & xdm = p5 * (xicr + xdl) ydm = p5 * ydl - area4 = p5 * (xic - xcr) * ydr * areafac_r(i,j) ! jlem could be improved... + area4 = p5 * (xic - xcr) * ydr * areafac_r(i,j) ! jlem should this be modified? area_c = edgearea(i,j) - area1 - area2 - area3 - area4 ! shift midpoint so that area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl ! jlem could be improved... + w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl ! jlem should this be modified? w2 = (xic-xdl)**2 + ydl**2 w1 = w1/w2 xdm = xdm - ydl*w1 From c4b573f25873e351e06766a167cc47c23e3c6af3 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Wed, 9 Aug 2023 16:23:41 +0000 Subject: [PATCH 05/18] Modified areafac calc for case 1 and case 2 --- .../cicedyn/dynamics/ice_transport_remap.F90 | 219 +++++++++--------- 1 file changed, 113 insertions(+), 106 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 9aad9ed8b..a2427f639 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -32,7 +32,7 @@ module ice_transport_remap use ice_blocks, only: nx_block, ny_block use ice_calendar, only: istep1 use ice_communicate, only: my_task - use ice_constants, only: c0, c1, c2, c12, p333, p4, p5, p6, & + use ice_constants, only: c0, c1, c2, c12, p25, p333, p4, p5, p6, & eps13, eps16, & field_loc_center, field_type_scalar, & field_loc_NEcorner, field_type_vector @@ -1770,7 +1770,6 @@ subroutine locate_triangles (nx_block, ny_block, & areafac_l , & ! area scale factor at left corner areafac_r , & ! area scale factor at right corner areafac_c , & ! (areafac_l + areafac_r)/2 - areafacS , & ! Sum of areafac_l and areafac_r areafacD ! Difference of areafac_l and areafac_r real (kind=dbl_kind) :: & @@ -1792,7 +1791,9 @@ subroutine locate_triangles (nx_block, ny_block, & area1, area2 , & ! temporary triangle areas area3, area4 , & ! area_c , & ! center polygon area - areatp, xpm , & ! to calculate areafact + areatp, , & ! temporary area calculation + areatplone , & ! temporary area calculation for lone triangle + xpm , & ! mean value of x points puny , & ! w1, w2 ! work variables @@ -1865,7 +1866,6 @@ subroutine locate_triangles (nx_block, ny_block, & areafac_l(:,:) = c0 areafac_r(:,:) = c0 areafac_c(:,:) = c0 - areafacS(:,:) = c0 areafacD(:,:) = c0 do ng = 1, ngroups do j = 1, ny_block @@ -1917,7 +1917,6 @@ subroutine locate_triangles (nx_block, ny_block, & areafac_l(i,j) = dxu(i-1,j)*dyu(i-1,j) areafac_r(i,j) = dxu(i ,j)*dyu(i ,j) areafac_c(i,j) = p5*(areafac_l(i,j) + areafac_r(i,j)) - areafacS (i,j) = areafac_r(i,j) + areafac_l(i,j) areafacD (i,j) = areafac_r(i,j) - areafac_l(i,j) enddo enddo @@ -1953,7 +1952,6 @@ subroutine locate_triangles (nx_block, ny_block, & areafac_l(i,j) = dxu(i,j )*dyu(i,j ) areafac_r(i,j) = dxu(i,j-1)*dyu(i,j-1) areafac_c(i,j) = p5*(areafac_l(i,j) + areafac_r(i,j)) - areafacS (i,j) = areafac_r(i,j) + areafac_l(i,j) areafacD (i,j) = areafac_r(i,j) - areafac_l(i,j) enddo enddo @@ -2080,7 +2078,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_l(i,j) @@ -2098,7 +2096,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_l(i,j) @@ -2116,7 +2114,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_l(i,j) @@ -2132,7 +2130,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_l(i,j) @@ -2150,7 +2148,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_l(i,j) @@ -2166,7 +2164,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_l(i,j) @@ -2191,7 +2189,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_r(i,j) @@ -2209,7 +2207,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_r(i,j) @@ -2227,7 +2225,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_r(i,j) @@ -2243,7 +2241,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_r(i,j) @@ -2261,7 +2259,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_r(i,j) @@ -2277,7 +2275,7 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_r(i,j) @@ -2347,8 +2345,12 @@ subroutine locate_triangles (nx_block, ny_block, & ! compute required area of central departure region area_c = edgearea(i,j) - area1 - area2 - area3 + ! compute areafac for quadrilateral cl-cr-dr-dl + xpm = p25*(xcl + xcr + xdr + xdl) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! shift midpoint so that the area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) & ! jlem should this be modified? + w1 = c2*area_c/areatp & + (xdr-xcl)*ydl + (xcr-xdl)*ydr w2 = (xdr-xdl)**2 + (ydr-ydl)**2 w1 = w1/w2 @@ -2382,12 +2384,18 @@ subroutine locate_triangles (nx_block, ny_block, & xdm = p5 * (xdr + xicl) ydm = p5 * ydr - ! compute area of triangle adjacent to left corner - area4 = p5 * (xcl - xic) * ydl * areafac_l(i,j) ! jlem should this be modified? + ! compute area of (lone) triangle adjacent to left corner + xpm = p33*(xcl + xdl + xic) + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + area4 = p5 * (xcl - xic) * ydl * areatplone area_c = edgearea(i,j) - area1 - area2 - area3 - area4 + ! compute areafac for right triangle ic-cr-dr + xpm = p33*(xic + xcr + xdr) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! shift midpoint so that area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr ! jlem should this be modified? + w1 = c2*area_c/areatp + (xcr-xic)*ydr w2 = (xdr-xic)**2 + ydr**2 w1 = w1/w2 xdm = xdm + ydr*w1 @@ -2411,11 +2419,18 @@ subroutine locate_triangles (nx_block, ny_block, & xdm = p5 * (xicr + xdl) ydm = p5 * ydl - area4 = p5 * (xic - xcr) * ydr * areafac_r(i,j) ! jlem should this be modified? + ! compute area of (lone) triangle adjacent to right corner + xpm = p33*(xic + xcr + xdr) + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + area4 = p5 * (xic - xcr) * ydr * areatplone area_c = edgearea(i,j) - area1 - area2 - area3 - area4 + ! compute areafac for left triangle cl-dl-ic + xpm = p33*(xcl + xdl + xic) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! shift midpoint so that area of remaining triangles = area_c - w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl ! jlem should this be modified? + w1 = c2*area_c/areatp + (xic-xcl)*ydl w2 = (xic-xdl)**2 + ydl**2 w1 = w1/w2 xdm = xdm - ydl*w1 @@ -2444,6 +2459,10 @@ subroutine locate_triangles (nx_block, ny_block, & if (ydl >= c0 .and. ydr >= c0 .and. ydm >= c0) then + ! compute areafac for quadrilateral cl-cr-dr-dl + xpm = p25*(xcl + xcr + xdr + xdl) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! TC1a (group 4) ng = 4 @@ -2455,8 +2474,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) @@ -2471,8 +2488,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) @@ -2487,13 +2502,15 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl >= c0 .and. ydr >= c0 .and. ydm < c0) then ! rare + ! compute areafac for quadrilateral cl-cr-dr-dl + xpm = p25*(xcl + xcr + xdr + xdl) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! TC1b (group 4) ng = 4 @@ -2505,8 +2522,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) @@ -2521,8 +2536,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) @@ -2537,13 +2550,15 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr < c0 .and. ydm < c0) then + ! compute areafac for quadrilateral cl-cr-dr-dl + xpm = p25*(xcl + xcr + xdr + xdl) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! BC1a (group 4) ng = 4 @@ -2555,8 +2570,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ycr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) @@ -2571,8 +2584,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) @@ -2587,13 +2598,15 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr < c0 .and. ydm >= c0) then ! rare + ! compute areafac for quadrilateral cl-cr-dr-dl + xpm = p25*(xcl + xcr + xdr + xdl) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! BC1b (group 4) ng = 4 @@ -2605,8 +2618,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) @@ -2621,8 +2632,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) @@ -2637,8 +2646,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) @@ -2649,6 +2656,10 @@ subroutine locate_triangles (nx_block, ny_block, & elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 & .and. ydm >= c0) then + ! compute areafac for left triangle cl-dl-ic + xpm = p33*(xcl + xdl + xic) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! TC1b (group 4) ng = 4 @@ -2660,12 +2671,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) - ! BC2b (group 5) + ! BC2b (group 5) lone triangle ng = 5 xp (i,j,1,ng) = xcr @@ -2677,8 +2686,8 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) - areafact(i,j,ng) = areatp + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + areafact(i,j,ng) = areatplone !areafact(i,j,ng) = areafac_r(i,j) ! TC3b (group 6) @@ -2692,14 +2701,16 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 & .and. ydm < c0 ) then ! less common + ! compute areafac for left triangle cl-dl-ic + xpm = p33*(xcl + xdl + xic) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! TC1b (group 4) ng = 4 @@ -2711,12 +2722,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) - ! BC2b (group 5) + ! BC2b (group 5) lone triangle ng = 5 xp (i,j,1,ng) = xcr @@ -2728,8 +2737,8 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) - areafact(i,j,ng) = areatp + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + areafact(i,j,ng) = areatplone !areafact(i,j,ng) = areafac_r(i,j) ! BC3b (group 6) @@ -2743,15 +2752,17 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 & .and. ydm < c0) then - ! TC1b (group 4) + ! compute areafac for right triangle ic-cr-dr + xpm = p33*(xic + xcr + xdr) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + + ! TC1b (group 4) lone triangle ng = 4 xp (i,j,1,ng) = xcl @@ -2763,8 +2774,8 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) - areafact(i,j,ng) = -areatp + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + areafact(i,j,ng) = -areatplone !areafact(i,j,ng) = -areafac_l(i,j) ! BC2b (group 5) @@ -2778,8 +2789,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) @@ -2794,15 +2803,17 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 & .and. ydm >= c0) then ! less common - ! TC1b (group 4) + ! compute areafac for right triangle ic-cr-dr + xpm = p33*(xic + xcr + xdr) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + + ! TC1b (group 4) lone triangle ng = 4 xp (i,j,1,ng) = xcl @@ -2814,8 +2825,8 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) - areafact(i,j,ng) = -areatp + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + areafact(i,j,ng) = -areatplone !areafact(i,j,ng) = -areafac_l(i,j) ! BC2b (group 5) @@ -2829,8 +2840,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) @@ -2845,15 +2854,17 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 & .and. ydm >= c0) then - ! BC1b (group 4) + ! compute areafac for right triangle ic-cr-dr + xpm = p33*(xic + xcr + xdr) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + + ! BC1b (group 4) lone triangle ng = 4 xp (i,j,1,ng) = xcl @@ -2865,8 +2876,8 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) - areafact(i,j,ng) = areatp + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + areafact(i,j,ng) = areatplone !areafact(i,j,ng) = areafac_l(i,j) ! TC2b (group 5) @@ -2880,8 +2891,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) @@ -2896,15 +2905,17 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 & .and. ydm < c0) then ! less common - ! BC1b (group 4) + ! compute areafac for right triangle ic-cr-dr + xpm = p33*(xic + xcr + xdr) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + + ! BC1b (group 4) lone triangle ng = 4 xp (i,j,1,ng) = xcl @@ -2916,8 +2927,8 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) - areafact(i,j,ng) = areatp + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + areafact(i,j,ng) = areatplone !areafact(i,j,ng) = areafac_l(i,j) ! TC2b (group 5) @@ -2931,8 +2942,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) @@ -2947,14 +2956,16 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 & .and. ydm < c0) then + ! compute areafac for left triangle cl-dl-ic + xpm = p33*(xcl + xdl + xic) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! BC1b (group 4) ng = 4 @@ -2966,12 +2977,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - ! TC2b (group 5) + ! TC2b (group 5) lone triangle ng = 5 xp (i,j,1,ng) = xcr @@ -2983,8 +2992,8 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) - areafact(i,j,ng) = -areatp + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + areafact(i,j,ng) = -areatplone !areafact(i,j,ng) = -areafac_r(i,j) ! BC3b (group 6) @@ -2998,14 +3007,16 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 & .and. ydm >= c0) then ! less common + ! compute areafac for left triangle cl-dl-ic + xpm = p33*(xcl + xdl + xic) + areatp = areafacD(i,j)*xpm + areafac_c(i,j) + ! BC1b (group 4) ng = 4 @@ -3017,12 +3028,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - ! TC2b (group 5) + ! TC2b (group 5) lone triangle ng = 5 xp (i,j,1,ng) = xcr @@ -3034,8 +3043,8 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) - areafact(i,j,ng) = -areatp + areatplone = areafacD(i,j)*xpm + areafac_c(i,j) + areafact(i,j,ng) = -areatplone !areafact(i,j,ng) = -areafac_r(i,j) ! TC3b (group 6) @@ -3049,8 +3058,6 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + p5*areafacS(i,j) areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) @@ -3094,7 +3101,7 @@ subroutine locate_triangles (nx_block, ny_block, & (xp(i,j,3,ng)-xp(i,j,1,ng)) ) & * areafact(i,j,ng) - if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then + if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then ! jf ici triarea(i,j,ng) = c0 else icells(ng) = icells(ng) + 1 @@ -3113,7 +3120,7 @@ subroutine locate_triangles (nx_block, ny_block, & do ij = 1, icellsd i = indxid(ij) j = indxjd(ij) - if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then + if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then ! jf ici write(nu_diag,*) '' write(nu_diag,*) 'Areas do not add up: m, i, j, edge =', & my_task, i, j, trim(edge) From dd6e12a49b2c0df1b1e6768ce95549798586e79e Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Wed, 9 Aug 2023 20:30:26 +0000 Subject: [PATCH 06/18] Corrected minor compilation issues --- .../cicedyn/dynamics/ice_transport_remap.F90 | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index a2427f639..bef542624 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -1791,7 +1791,7 @@ subroutine locate_triangles (nx_block, ny_block, & area1, area2 , & ! temporary triangle areas area3, area4 , & ! area_c , & ! center polygon area - areatp, , & ! temporary area calculation + areatp , & ! temporary area calculation areatplone , & ! temporary area calculation for lone triangle xpm , & ! mean value of x points puny , & ! @@ -2385,13 +2385,13 @@ subroutine locate_triangles (nx_block, ny_block, & ydm = p5 * ydr ! compute area of (lone) triangle adjacent to left corner - xpm = p33*(xcl + xdl + xic) + xpm = p333*(xcl + xdl + xic) areatplone = areafacD(i,j)*xpm + areafac_c(i,j) area4 = p5 * (xcl - xic) * ydl * areatplone area_c = edgearea(i,j) - area1 - area2 - area3 - area4 ! compute areafac for right triangle ic-cr-dr - xpm = p33*(xic + xcr + xdr) + xpm = p333*(xic + xcr + xdr) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! shift midpoint so that area of remaining triangles = area_c @@ -2420,13 +2420,13 @@ subroutine locate_triangles (nx_block, ny_block, & ydm = p5 * ydl ! compute area of (lone) triangle adjacent to right corner - xpm = p33*(xic + xcr + xdr) + xpm = p333*(xic + xcr + xdr) areatplone = areafacD(i,j)*xpm + areafac_c(i,j) area4 = p5 * (xic - xcr) * ydr * areatplone area_c = edgearea(i,j) - area1 - area2 - area3 - area4 ! compute areafac for left triangle cl-dl-ic - xpm = p33*(xcl + xdl + xic) + xpm = p333*(xcl + xdl + xic) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! shift midpoint so that area of remaining triangles = area_c @@ -2657,7 +2657,7 @@ subroutine locate_triangles (nx_block, ny_block, & .and. ydm >= c0) then ! compute areafac for left triangle cl-dl-ic - xpm = p33*(xcl + xdl + xic) + xpm = p333*(xcl + xdl + xic) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! TC1b (group 4) @@ -2708,7 +2708,7 @@ subroutine locate_triangles (nx_block, ny_block, & .and. ydm < c0 ) then ! less common ! compute areafac for left triangle cl-dl-ic - xpm = p33*(xcl + xdl + xic) + xpm = p333*(xcl + xdl + xic) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! TC1b (group 4) @@ -2759,7 +2759,7 @@ subroutine locate_triangles (nx_block, ny_block, & .and. ydm < c0) then ! compute areafac for right triangle ic-cr-dr - xpm = p33*(xic + xcr + xdr) + xpm = p333*(xic + xcr + xdr) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! TC1b (group 4) lone triangle @@ -2810,7 +2810,7 @@ subroutine locate_triangles (nx_block, ny_block, & .and. ydm >= c0) then ! less common ! compute areafac for right triangle ic-cr-dr - xpm = p33*(xic + xcr + xdr) + xpm = p333*(xic + xcr + xdr) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! TC1b (group 4) lone triangle @@ -2861,7 +2861,7 @@ subroutine locate_triangles (nx_block, ny_block, & .and. ydm >= c0) then ! compute areafac for right triangle ic-cr-dr - xpm = p33*(xic + xcr + xdr) + xpm = p333*(xic + xcr + xdr) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! BC1b (group 4) lone triangle @@ -2912,7 +2912,7 @@ subroutine locate_triangles (nx_block, ny_block, & .and. ydm < c0) then ! less common ! compute areafac for right triangle ic-cr-dr - xpm = p33*(xic + xcr + xdr) + xpm = p333*(xic + xcr + xdr) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! BC1b (group 4) lone triangle @@ -2963,7 +2963,7 @@ subroutine locate_triangles (nx_block, ny_block, & .and. ydm < c0) then ! compute areafac for left triangle cl-dl-ic - xpm = p33*(xcl + xdl + xic) + xpm = p333*(xcl + xdl + xic) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! BC1b (group 4) @@ -3014,7 +3014,7 @@ subroutine locate_triangles (nx_block, ny_block, & .and. ydm >= c0) then ! less common ! compute areafac for left triangle cl-dl-ic - xpm = p33*(xcl + xdl + xic) + xpm = p333*(xcl + xdl + xic) areatp = areafacD(i,j)*xpm + areafac_c(i,j) ! BC1b (group 4) From 978fdf382f09254f9c759175297d2db99089d44f Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Thu, 10 Aug 2023 15:24:57 +0000 Subject: [PATCH 07/18] Corrected conditions for case 1 to make sure areas add up --- cicecore/cicedyn/dynamics/ice_transport_remap.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index bef542624..4c6f4309b 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -2553,7 +2553,7 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl < c0 .and. ydr < c0 .and. ydm < c0) then + elseif (ydl <= c0 .and. ydr <= c0 .and. ydm < c0) then ! compute areafac for quadrilateral cl-cr-dr-dl xpm = p25*(xcl + xcr + xdr + xdl) @@ -2601,7 +2601,7 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl < c0 .and. ydr < c0 .and. ydm >= c0) then ! rare + elseif (ydl <= c0 .and. ydr <= c0 .and. ydm >= c0) then ! rare ! compute areafac for quadrilateral cl-cr-dr-dl xpm = p25*(xcl + xcr + xdr + xdl) From 8fe30953eaebc6079047566dd0ce2eee96f08e21 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Thu, 10 Aug 2023 15:31:13 +0000 Subject: [PATCH 08/18] Small modif in l_fixed_area section to ensure only one condition is true --- cicecore/cicedyn/dynamics/ice_transport_remap.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 4c6f4309b..c59e82155 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -2375,7 +2375,7 @@ subroutine locate_triangles (nx_block, ny_block, & endif yicr = c0 - elseif (xic < c0) then ! fix ICL = IC + elseif (xic < c0 .and. xic > xcl) then ! fix ICL = IC ! Bill please check this xicl = xic yicl = yic @@ -2410,7 +2410,7 @@ subroutine locate_triangles (nx_block, ny_block, & endif yicr = c0 - elseif (xic >= c0) then ! fix ICR = IR + elseif (xic >= c0 .and. xic < xcr) then ! fix ICR = IR ! Bill please check this xicr = xic yicr = yic @@ -3101,7 +3101,7 @@ subroutine locate_triangles (nx_block, ny_block, & (xp(i,j,3,ng)-xp(i,j,1,ng)) ) & * areafact(i,j,ng) - if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then ! jf ici + if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then ! Bill please check this triarea(i,j,ng) = c0 else icells(ng) = icells(ng) + 1 @@ -3120,7 +3120,7 @@ subroutine locate_triangles (nx_block, ny_block, & do ij = 1, icellsd i = indxid(ij) j = indxjd(ij) - if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then ! jf ici + if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then ! Bill please check this write(nu_diag,*) '' write(nu_diag,*) 'Areas do not add up: m, i, j, edge =', & my_task, i, j, trim(edge) From 0ef0ac0df9ee9ca14705abbac235f2e2d74df662 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Thu, 10 Aug 2023 19:19:34 +0000 Subject: [PATCH 09/18] Modified conditions in locate triangle to be consistent with previous changes for case 1 --- .../cicedyn/dynamics/ice_transport_remap.F90 | 36 +++++++++---------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index c59e82155..3c78acc84 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -2553,7 +2553,7 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl <= c0 .and. ydr <= c0 .and. ydm < c0) then + elseif (ydl <= c0 .and. ydr <= c0 .and. ydm <= c0) then ! Bill I also changed ydm < to <= ! compute areafac for quadrilateral cl-cr-dr-dl xpm = p25*(xcl + xcr + xdr + xdl) @@ -2601,7 +2601,7 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl <= c0 .and. ydr <= c0 .and. ydm >= c0) then ! rare + elseif (ydl <= c0 .and. ydr <= c0 .and. ydm > c0) then ! rare ! Bill I also changed ydm >= to > ! compute areafac for quadrilateral cl-cr-dr-dl xpm = p25*(xcl + xcr + xdr + xdl) @@ -2653,8 +2653,8 @@ subroutine locate_triangles (nx_block, ny_block, & ! For these cases, one triangle is given the area factor associated ! with the adjacent corner, to avoid rare negative masses on curved grids. - elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 & - .and. ydm >= c0) then + elseif (ydl > c0 .and. ydr < c0 .and. xic >= c0 & + .and. ydm >= c0) then ! compute areafac for left triangle cl-dl-ic xpm = p333*(xcl + xdl + xic) @@ -2704,8 +2704,8 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) - elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 & - .and. ydm < c0 ) then ! less common + elseif (ydl > c0 .and. ydr < c0 .and. xic >= c0 & + .and. ydm < c0 ) then ! less common ! compute areafac for left triangle cl-dl-ic xpm = p333*(xcl + xdl + xic) @@ -2755,8 +2755,8 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 & - .and. ydm < c0) then + elseif (ydl > c0 .and. ydr < c0 .and. xic < c0 & + .and. ydm < c0) then ! compute areafac for right triangle ic-cr-dr xpm = p333*(xic + xcr + xdr) @@ -2806,8 +2806,8 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 & - .and. ydm >= c0) then ! less common + elseif (ydl > c0 .and. ydr < c0 .and. xic < c0 & + .and. ydm >= c0) then ! less common ! compute areafac for right triangle ic-cr-dr xpm = p333*(xic + xcr + xdr) @@ -2857,8 +2857,8 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) - elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 & - .and. ydm >= c0) then + elseif (ydl < c0 .and. ydr > c0 .and. xic < c0 & + .and. ydm >= c0) then ! compute areafac for right triangle ic-cr-dr xpm = p333*(xic + xcr + xdr) @@ -2908,8 +2908,8 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = -areatp !areafact(i,j,ng) = -areafac_c(i,j) - elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 & - .and. ydm < c0) then ! less common + elseif (ydl < c0 .and. ydr > c0 .and. xic < c0 & + .and. ydm < c0) then ! less common ! compute areafac for right triangle ic-cr-dr xpm = p333*(xic + xcr + xdr) @@ -2959,8 +2959,8 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 & - .and. ydm < c0) then + elseif (ydl < c0 .and. ydr > c0 .and. xic >= c0 & + .and. ydm < c0) then ! compute areafac for left triangle cl-dl-ic xpm = p333*(xcl + xdl + xic) @@ -3010,8 +3010,8 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = areatp !areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 & - .and. ydm >= c0) then ! less common + elseif (ydl < c0 .and. ydr > c0 .and. xic >= c0 & + .and. ydm >= c0) then ! less common ! compute areafac for left triangle cl-dl-ic xpm = p333*(xcl + xdl + xic) From 69a6374b0785131ac1c61ff83aea1ee9e53a9e57 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Wed, 23 Aug 2023 15:08:32 +0000 Subject: [PATCH 10/18] Use other edge areafac_c for TL, BL, TR and BR triangles --- .../cicedyn/dynamics/ice_transport_remap.F90 | 344 +++++++----------- 1 file changed, 122 insertions(+), 222 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 3c78acc84..59cd46e55 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -63,7 +63,7 @@ module ice_transport_remap ! if false, area flux is determined internally ! and is passed out - logical (kind=log_kind), parameter :: bugcheck = .false. + logical (kind=log_kind), parameter :: bugcheck = .true. !======================================================================= ! Here is some information about how the incremental remapping scheme @@ -1756,7 +1756,11 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_tr, jshift_tr , & ! i,j indices of TR cell relative to edge ishift_br, jshift_br , & ! i,j indices of BR cell relative to edge ishift_tc, jshift_tc , & ! i,j indices of TC cell relative to edge - ishift_bc, jshift_bc ! i,j indices of BC cell relative to edge + ishift_bc, jshift_bc , & ! i,j indices of BC cell relative to edge + ise_tl, jse_tl , & ! i,j of TL other edge relative to edge + ise_bl, jse_bl , & ! i,j of BL other edge relative to edge + ise_tr, jse_tr , & ! i,j of TR other edge relative to edge + ise_br, jse_br ! i,j of BR other edge relative to edge integer (kind=int_kind) :: & icellsd ! number of cells where departure area > 0. @@ -1767,12 +1771,12 @@ subroutine locate_triangles (nx_block, ny_block, & real (kind=dbl_kind), dimension(nx_block,ny_block) :: & dx, dy , & ! scaled departure points - areafac_l , & ! area scale factor at left corner - areafac_r , & ! area scale factor at right corner areafac_c , & ! (areafac_l + areafac_r)/2 - areafacD ! Difference of areafac_l and areafac_r + areafac_ce ! areafac_c on other edge real (kind=dbl_kind) :: & + areafac_l , & ! area scale factor at left corner + areafac_r , & ! area scale factor at right corner xcl, ycl , & ! coordinates of left corner point ! (relative to midpoint of edge) xdl, ydl , & ! left departure point @@ -1791,9 +1795,6 @@ subroutine locate_triangles (nx_block, ny_block, & area1, area2 , & ! temporary triangle areas area3, area4 , & ! area_c , & ! center polygon area - areatp , & ! temporary area calculation - areatplone , & ! temporary area calculation for lone triangle - xpm , & ! mean value of x points puny , & ! w1, w2 ! work variables @@ -1863,10 +1864,9 @@ subroutine locate_triangles (nx_block, ny_block, & if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - areafac_l(:,:) = c0 - areafac_r(:,:) = c0 - areafac_c(:,:) = c0 - areafacD(:,:) = c0 + areafac_c(:,:) = c0 + areafac_ce(:,:) = c0 + do ng = 1, ngroups do j = 1, ny_block do i = 1, nx_block @@ -1910,14 +1910,34 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_bc = 0 jshift_bc = 0 + ! index shifts for neighbor east edges + + ise_tl = -1 + jse_tl = 1 + ise_bl = -1 + jse_bl = 0 + ise_tr = 0 + jse_tr = 1 + ise_br = 0 + jse_br = 0 + ! area scale factor do j = jb, je do i = ib, ie - areafac_l(i,j) = dxu(i-1,j)*dyu(i-1,j) - areafac_r(i,j) = dxu(i ,j)*dyu(i ,j) - areafac_c(i,j) = p5*(areafac_l(i,j) + areafac_r(i,j)) - areafacD (i,j) = areafac_r(i,j) - areafac_l(i,j) + areafac_l = dxu(i-1,j)*dyu(i-1,j) + areafac_r = dxu(i ,j)*dyu(i ,j) + areafac_c(i,j) = p5*( areafac_l + areafac_r ) + enddo + enddo + + ! area scale factor for other edge (east) + + do j = jb, je+1 + do i = ib-1, ie + areafac_l = dxu(i,j )*dyu(i,j ) + areafac_r = dxu(i,j-1)*dyu(i,j-1) + areafac_ce(i,j) = p5*( areafac_l + areafac_r ) enddo enddo @@ -1945,14 +1965,34 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_bc = 0 jshift_bc = 0 + ! index shifts for neighbor north edges + + ise_tl = 1 + jse_tl = 0 + ise_bl = 0 + jse_bl = 0 + ise_tr = 1 + jse_tr = -1 + ise_br = 0 + jse_br = -1 + ! area scale factors do j = jb, je do i = ib, ie - areafac_l(i,j) = dxu(i,j )*dyu(i,j ) - areafac_r(i,j) = dxu(i,j-1)*dyu(i,j-1) - areafac_c(i,j) = p5*(areafac_l(i,j) + areafac_r(i,j)) - areafacD (i,j) = areafac_r(i,j) - areafac_l(i,j) + areafac_l = dxu(i,j )*dyu(i,j ) + areafac_r = dxu(i,j-1)*dyu(i,j-1) + areafac_c(i,j) = p5*( areafac_l + areafac_r ) + enddo + enddo + + ! area scale factor for ther edge (north) + + do j = jb-1, je + do i = ib, ie+1 + areafac_l = dxu(i-1,j)*dyu(i-1,j) + areafac_r = dxu(i ,j)*dyu(i ,j) + areafac_ce(i,j) = p5*( areafac_l + areafac_r ) enddo enddo @@ -2077,9 +2117,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatp + areafact(i,j,ng) = -areafac_ce(i+ise_tl,j+jse_tl) !areafact(i,j,ng) = -areafac_l(i,j) elseif (yil < c0 .and. xdl < xcl .and. ydl < c0) then @@ -2095,9 +2133,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yil iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatp + areafact(i,j,ng) = areafac_ce(i+ise_bl,j+jse_bl) !areafact(i,j,ng) = areafac_l(i,j) elseif (yil < c0 .and. xdl < xcl .and. ydl >= c0) then @@ -2113,9 +2149,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatp + areafact(i,j,ng) = areafac_ce(i+ise_tl,j+jse_tl) !areafact(i,j,ng) = areafac_l(i,j) ! BL1 (group 3) @@ -2129,9 +2163,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yil iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatp + areafact(i,j,ng) = areafac_ce(i+ise_bl,j+jse_bl) !areafact(i,j,ng) = areafac_l(i,j) elseif (yil > c0 .and. xdl < xcl .and. ydl < c0) then @@ -2147,9 +2179,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatp + areafact(i,j,ng) = -areafac_ce(i+ise_tl,j+jse_tl) !areafact(i,j,ng) = -areafac_l(i,j) ! BL2 (group 1) @@ -2163,9 +2193,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatp + areafact(i,j,ng) = -areafac_ce(i+ise_bl,j+jse_bl) !areafact(i,j,ng) = -areafac_l(i,j) endif ! TL and BL triangles @@ -2188,9 +2216,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yir iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatp + areafact(i,j,ng) = -areafac_ce(i+ise_tr,j+jse_tr) !areafact(i,j,ng) = -areafac_r(i,j) elseif (yir < c0 .and. xdr >= xcr .and. ydr < c0) then @@ -2206,9 +2232,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatp + areafact(i,j,ng) = areafac_ce(i+ise_br,j+jse_br) !areafact(i,j,ng) = areafac_r(i,j) elseif (yir < c0 .and. xdr >= xcr .and. ydr >= c0) then @@ -2224,9 +2248,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatp + areafact(i,j,ng) = areafac_ce(i+ise_tr,j+jse_tr) !areafact(i,j,ng) = areafac_r(i,j) ! BR1 (group 3) @@ -2240,9 +2262,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatp + areafact(i,j,ng) = areafac_ce(i+ise_br,j+jse_br) !areafact(i,j,ng) = areafac_r(i,j) elseif (yir > c0 .and. xdr >= xcr .and. ydr < c0) then @@ -2258,9 +2278,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yir iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatp + areafact(i,j,ng) = -areafac_ce(i+ise_tr,j+jse_tr) !areafact(i,j,ng) = -areafac_r(i,j) ! BR2 (group 2) @@ -2274,9 +2292,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatp + areafact(i,j,ng) = -areafac_ce(i+ise_br,j+jse_br) !areafact(i,j,ng) = -areafac_r(i,j) endif ! TR and BR triangles @@ -2345,12 +2361,8 @@ subroutine locate_triangles (nx_block, ny_block, & ! compute required area of central departure region area_c = edgearea(i,j) - area1 - area2 - area3 - ! compute areafac for quadrilateral cl-cr-dr-dl - xpm = p25*(xcl + xcr + xdr + xdl) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! shift midpoint so that the area of remaining triangles = area_c - w1 = c2*area_c/areatp & + w1 = c2*area_c/areafac_c(i,j) & + (xdr-xcl)*ydl + (xcr-xdl)*ydr w2 = (xdr-xdl)**2 + (ydr-ydl)**2 w1 = w1/w2 @@ -2375,7 +2387,7 @@ subroutine locate_triangles (nx_block, ny_block, & endif yicr = c0 - elseif (xic < c0 .and. xic > xcl) then ! fix ICL = IC ! Bill please check this + elseif (xic < c0 .and. xic > xcl) then ! fix ICL = IC xicl = xic yicl = yic @@ -2385,17 +2397,11 @@ subroutine locate_triangles (nx_block, ny_block, & ydm = p5 * ydr ! compute area of (lone) triangle adjacent to left corner - xpm = p333*(xcl + xdl + xic) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - area4 = p5 * (xcl - xic) * ydl * areatplone + area4 = p5 * (xcl - xic) * ydl * areafac_c(i,j) ! Bill check this area_c = edgearea(i,j) - area1 - area2 - area3 - area4 - ! compute areafac for right triangle ic-cr-dr - xpm = p333*(xic + xcr + xdr) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! shift midpoint so that area of remaining triangles = area_c - w1 = c2*area_c/areatp + (xcr-xic)*ydr + w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr w2 = (xdr-xic)**2 + ydr**2 w1 = w1/w2 xdm = xdm + ydr*w1 @@ -2410,7 +2416,7 @@ subroutine locate_triangles (nx_block, ny_block, & endif yicr = c0 - elseif (xic >= c0 .and. xic < xcr) then ! fix ICR = IR ! Bill please check this + elseif (xic >= c0 .and. xic < xcr) then ! fix ICR = IR xicr = xic yicr = yic @@ -2420,17 +2426,11 @@ subroutine locate_triangles (nx_block, ny_block, & ydm = p5 * ydl ! compute area of (lone) triangle adjacent to right corner - xpm = p333*(xic + xcr + xdr) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - area4 = p5 * (xic - xcr) * ydr * areatplone + area4 = p5 * (xic - xcr) * ydr * areafac_c(i,j) ! Bill check this area_c = edgearea(i,j) - area1 - area2 - area3 - area4 - ! compute areafac for left triangle cl-dl-ic - xpm = p333*(xcl + xdl + xic) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! shift midpoint so that area of remaining triangles = area_c - w1 = c2*area_c/areatp + (xic-xcl)*ydl + w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl w2 = (xic-xdl)**2 + ydl**2 w1 = w1/w2 xdm = xdm - ydl*w1 @@ -2459,10 +2459,6 @@ subroutine locate_triangles (nx_block, ny_block, & if (ydl >= c0 .and. ydr >= c0 .and. ydm >= c0) then - ! compute areafac for quadrilateral cl-cr-dr-dl - xpm = p25*(xcl + xcr + xdr + xdl) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! TC1a (group 4) ng = 4 @@ -2474,9 +2470,8 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) - + areafact(i,j,ng) = -areafac_c(i,j) + ! TC2a (group 5) ng = 5 @@ -2488,9 +2483,8 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) - + areafact(i,j,ng) = -areafac_c(i,j) + ! TC3a (group 6) ng = 6 @@ -2502,15 +2496,10 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl >= c0 .and. ydr >= c0 .and. ydm < c0) then ! rare - ! compute areafac for quadrilateral cl-cr-dr-dl - xpm = p25*(xcl + xcr + xdr + xdl) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! TC1b (group 4) ng = 4 @@ -2522,8 +2511,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! TC2b (group 5) @@ -2536,8 +2524,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! BC3b (group 6) @@ -2550,14 +2537,9 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) - - elseif (ydl <= c0 .and. ydr <= c0 .and. ydm <= c0) then ! Bill I also changed ydm < to <= + areafact(i,j,ng) = areafac_c(i,j) - ! compute areafac for quadrilateral cl-cr-dr-dl - xpm = p25*(xcl + xcr + xdr + xdl) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) + elseif (ydl <= c0 .and. ydr <= c0 .and. ydm <= c0) then ! BC1a (group 4) @@ -2570,8 +2552,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ycr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! BC2a (group 5) @@ -2584,8 +2565,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! BC3a (group 6) @@ -2598,14 +2578,9 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) - elseif (ydl <= c0 .and. ydr <= c0 .and. ydm > c0) then ! rare ! Bill I also changed ydm >= to > - - ! compute areafac for quadrilateral cl-cr-dr-dl - xpm = p25*(xcl + xcr + xdr + xdl) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) + elseif (ydl <= c0 .and. ydr <= c0 .and. ydm > c0) then ! rare ! BC1b (group 4) @@ -2618,8 +2593,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! BC2b (group 5) @@ -2632,8 +2606,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! TC3b (group 6) @@ -2646,20 +2619,15 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! Now consider cases where the two DPs lie in different grid cells ! For these cases, one triangle is given the area factor associated - ! with the adjacent corner, to avoid rare negative masses on curved grids. + ! with the adjacent corner, to avoid rare negative masses on curved grids. elseif (ydl > c0 .and. ydr < c0 .and. xic >= c0 & .and. ydm >= c0) then - ! compute areafac for left triangle cl-dl-ic - xpm = p333*(xcl + xdl + xic) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! TC1b (group 4) ng = 4 @@ -2671,8 +2639,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! BC2b (group 5) lone triangle @@ -2685,10 +2652,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatplone - !areafact(i,j,ng) = areafac_r(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! Bill check this and see comments above ! TC3b (group 6) @@ -2701,16 +2665,11 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl > c0 .and. ydr < c0 .and. xic >= c0 & .and. ydm < c0 ) then ! less common - ! compute areafac for left triangle cl-dl-ic - xpm = p333*(xcl + xdl + xic) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! TC1b (group 4) ng = 4 @@ -2722,8 +2681,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! BC2b (group 5) lone triangle @@ -2736,10 +2694,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatplone - !areafact(i,j,ng) = areafac_r(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! Bill check this ! BC3b (group 6) @@ -2752,16 +2707,11 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) elseif (ydl > c0 .and. ydr < c0 .and. xic < c0 & .and. ydm < c0) then - ! compute areafac for right triangle ic-cr-dr - xpm = p333*(xic + xcr + xdr) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! TC1b (group 4) lone triangle ng = 4 @@ -2773,10 +2723,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatplone - !areafact(i,j,ng) = -areafac_l(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! Bill check this ! BC2b (group 5) @@ -2789,8 +2736,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! BC3b (group 6) @@ -2803,16 +2749,11 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) elseif (ydl > c0 .and. ydr < c0 .and. xic < c0 & .and. ydm >= c0) then ! less common - ! compute areafac for right triangle ic-cr-dr - xpm = p333*(xic + xcr + xdr) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! TC1b (group 4) lone triangle ng = 4 @@ -2824,10 +2765,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatplone - !areafact(i,j,ng) = -areafac_l(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! Bill check this ! BC2b (group 5) @@ -2840,8 +2778,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! TC3b (group 6) @@ -2854,16 +2791,11 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl < c0 .and. ydr > c0 .and. xic < c0 & .and. ydm >= c0) then - ! compute areafac for right triangle ic-cr-dr - xpm = p333*(xic + xcr + xdr) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! BC1b (group 4) lone triangle ng = 4 @@ -2875,10 +2807,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatplone - !areafact(i,j,ng) = areafac_l(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! Bill check this ! TC2b (group 5) @@ -2891,8 +2820,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! TC3b (group 6) @@ -2905,16 +2833,11 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) elseif (ydl < c0 .and. ydr > c0 .and. xic < c0 & .and. ydm < c0) then ! less common - ! compute areafac for right triangle ic-cr-dr - xpm = p333*(xic + xcr + xdr) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! BC1b (group 4) lone triangle ng = 4 @@ -2926,10 +2849,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = areatplone - !areafact(i,j,ng) = areafac_l(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! Bill check this ! TC2b (group 5) @@ -2942,8 +2862,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! BC3b (group 6) @@ -2956,16 +2875,11 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr > c0 .and. xic >= c0 & .and. ydm < c0) then - ! compute areafac for left triangle cl-dl-ic - xpm = p333*(xcl + xdl + xic) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! BC1b (group 4) ng = 4 @@ -2977,8 +2891,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! TC2b (group 5) lone triangle @@ -2991,10 +2904,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatplone - !areafact(i,j,ng) = -areafac_r(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! Bill check this ! BC3b (group 6) @@ -3007,16 +2917,11 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) elseif (ydl < c0 .and. ydr > c0 .and. xic >= c0 & .and. ydm >= c0) then ! less common - ! compute areafac for left triangle cl-dl-ic - xpm = p333*(xcl + xdl + xic) - areatp = areafacD(i,j)*xpm + areafac_c(i,j) - ! BC1b (group 4) ng = 4 @@ -3028,8 +2933,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areatp - !areafact(i,j,ng) = areafac_c(i,j) + areafact(i,j,ng) = areafac_c(i,j) ! TC2b (group 5) lone triangle @@ -3042,10 +2946,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - xpm= p333*( xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng) ) - areatplone = areafacD(i,j)*xpm + areafac_c(i,j) - areafact(i,j,ng) = -areatplone - !areafact(i,j,ng) = -areafac_r(i,j) + areafact(i,j,ng) = -areafac_c(i,j) ! Bill check this ! TC3b (group 6) @@ -3058,8 +2959,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydm iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areatp - !areafact(i,j,ng) = -areafac_c(i,j) + areafact(i,j,ng) = -areafac_c(i,j) endif ! TC and BC triangles @@ -3101,7 +3001,7 @@ subroutine locate_triangles (nx_block, ny_block, & (xp(i,j,3,ng)-xp(i,j,1,ng)) ) & * areafact(i,j,ng) - if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then ! Bill please check this + if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then triarea(i,j,ng) = c0 else icells(ng) = icells(ng) + 1 @@ -3120,7 +3020,7 @@ subroutine locate_triangles (nx_block, ny_block, & do ij = 1, icellsd i = indxid(ij) j = indxjd(ij) - if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then ! Bill please check this + if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then write(nu_diag,*) '' write(nu_diag,*) 'Areas do not add up: m, i, j, edge =', & my_task, i, j, trim(edge) From 7c287c3ce9a59dc3abc40a30b390db7d868b598e Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Wed, 23 Aug 2023 17:31:55 +0000 Subject: [PATCH 11/18] Some comments removed --- .../cicedyn/dynamics/ice_transport_remap.F90 | 42 ++++++------------- 1 file changed, 13 insertions(+), 29 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 59cd46e55..61deb2487 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -32,7 +32,7 @@ module ice_transport_remap use ice_blocks, only: nx_block, ny_block use ice_calendar, only: istep1 use ice_communicate, only: my_task - use ice_constants, only: c0, c1, c2, c12, p25, p333, p4, p5, p6, & + use ice_constants, only: c0, c1, c2, c12, p333, p4, p5, p6, & eps13, eps16, & field_loc_center, field_type_scalar, & field_loc_NEcorner, field_type_vector @@ -1986,7 +1986,7 @@ subroutine locate_triangles (nx_block, ny_block, & enddo enddo - ! area scale factor for ther edge (north) + ! area scale factor for other edge (north) do j = jb-1, je do i = ib, ie+1 @@ -2118,7 +2118,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl areafact(i,j,ng) = -areafac_ce(i+ise_tl,j+jse_tl) - !areafact(i,j,ng) = -areafac_l(i,j) elseif (yil < c0 .and. xdl < xcl .and. ydl < c0) then @@ -2134,7 +2133,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl areafact(i,j,ng) = areafac_ce(i+ise_bl,j+jse_bl) - !areafact(i,j,ng) = areafac_l(i,j) elseif (yil < c0 .and. xdl < xcl .and. ydl >= c0) then @@ -2150,7 +2148,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl areafact(i,j,ng) = areafac_ce(i+ise_tl,j+jse_tl) - !areafact(i,j,ng) = areafac_l(i,j) ! BL1 (group 3) @@ -2164,7 +2161,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl areafact(i,j,ng) = areafac_ce(i+ise_bl,j+jse_bl) - !areafact(i,j,ng) = areafac_l(i,j) elseif (yil > c0 .and. xdl < xcl .and. ydl < c0) then @@ -2180,7 +2176,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl areafact(i,j,ng) = -areafac_ce(i+ise_tl,j+jse_tl) - !areafact(i,j,ng) = -areafac_l(i,j) ! BL2 (group 1) @@ -2194,7 +2189,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl areafact(i,j,ng) = -areafac_ce(i+ise_bl,j+jse_bl) - !areafact(i,j,ng) = -areafac_l(i,j) endif ! TL and BL triangles @@ -2217,7 +2211,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr areafact(i,j,ng) = -areafac_ce(i+ise_tr,j+jse_tr) - !areafact(i,j,ng) = -areafac_r(i,j) elseif (yir < c0 .and. xdr >= xcr .and. ydr < c0) then @@ -2233,7 +2226,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br areafact(i,j,ng) = areafac_ce(i+ise_br,j+jse_br) - !areafact(i,j,ng) = areafac_r(i,j) elseif (yir < c0 .and. xdr >= xcr .and. ydr >= c0) then @@ -2249,7 +2241,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr areafact(i,j,ng) = areafac_ce(i+ise_tr,j+jse_tr) - !areafact(i,j,ng) = areafac_r(i,j) ! BR1 (group 3) @@ -2263,7 +2254,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br areafact(i,j,ng) = areafac_ce(i+ise_br,j+jse_br) - !areafact(i,j,ng) = areafac_r(i,j) elseif (yir > c0 .and. xdr >= xcr .and. ydr < c0) then @@ -2279,7 +2269,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr areafact(i,j,ng) = -areafac_ce(i+ise_tr,j+jse_tr) - !areafact(i,j,ng) = -areafac_r(i,j) ! BR2 (group 2) @@ -2293,7 +2282,6 @@ subroutine locate_triangles (nx_block, ny_block, & iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br areafact(i,j,ng) = -areafac_ce(i+ise_br,j+jse_br) - !areafact(i,j,ng) = -areafac_r(i,j) endif ! TR and BR triangles @@ -2349,9 +2337,7 @@ subroutine locate_triangles (nx_block, ny_block, & ! region so that the sum of all triangle areas is equal to the ! prescribed value. ! If two triangles are in one grid cell and one is in the other, - ! then compute the area of the lone triangle using an area factor - ! corresponding to the adjacent corner. This is necessary to prevent - ! negative masses in some rare cases on curved grids. Then adjust + ! then compute the area of the lone triangle. Then adjust ! the area of the remaining two-triangle region so that the sum of ! all triangle areas has the prescribed value. !----------------------------------------------------------- @@ -2397,7 +2383,7 @@ subroutine locate_triangles (nx_block, ny_block, & ydm = p5 * ydr ! compute area of (lone) triangle adjacent to left corner - area4 = p5 * (xcl - xic) * ydl * areafac_c(i,j) ! Bill check this + area4 = p5 * (xcl - xic) * ydl * areafac_c(i,j) area_c = edgearea(i,j) - area1 - area2 - area3 - area4 ! shift midpoint so that area of remaining triangles = area_c @@ -2426,7 +2412,7 @@ subroutine locate_triangles (nx_block, ny_block, & ydm = p5 * ydl ! compute area of (lone) triangle adjacent to right corner - area4 = p5 * (xic - xcr) * ydr * areafac_c(i,j) ! Bill check this + area4 = p5 * (xic - xcr) * ydr * areafac_c(i,j) area_c = edgearea(i,j) - area1 - area2 - area3 - area4 ! shift midpoint so that area of remaining triangles = area_c @@ -2622,8 +2608,6 @@ subroutine locate_triangles (nx_block, ny_block, & areafact(i,j,ng) = -areafac_c(i,j) ! Now consider cases where the two DPs lie in different grid cells - ! For these cases, one triangle is given the area factor associated - ! with the adjacent corner, to avoid rare negative masses on curved grids. elseif (ydl > c0 .and. ydr < c0 .and. xic >= c0 & .and. ydm >= c0) then @@ -2652,7 +2636,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) ! Bill check this and see comments above + areafact(i,j,ng) = areafac_c(i,j) ! TC3b (group 6) @@ -2694,7 +2678,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) ! Bill check this + areafact(i,j,ng) = areafac_c(i,j) ! BC3b (group 6) @@ -2723,7 +2707,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) ! Bill check this + areafact(i,j,ng) = -areafac_c(i,j) ! BC2b (group 5) @@ -2765,7 +2749,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) ! Bill check this + areafact(i,j,ng) = -areafac_c(i,j) ! BC2b (group 5) @@ -2807,7 +2791,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) ! Bill check this + areafact(i,j,ng) = areafac_c(i,j) ! TC2b (group 5) @@ -2849,7 +2833,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicl iflux (i,j,ng) = i + ishift_bc jflux (i,j,ng) = j + jshift_bc - areafact(i,j,ng) = areafac_c(i,j) ! Bill check this + areafact(i,j,ng) = areafac_c(i,j) ! TC2b (group 5) @@ -2904,7 +2888,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) ! Bill check this + areafact(i,j,ng) = -areafac_c(i,j) ! BC3b (group 6) @@ -2946,7 +2930,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yicr iflux (i,j,ng) = i + ishift_tc jflux (i,j,ng) = j + jshift_tc - areafact(i,j,ng) = -areafac_c(i,j) ! Bill check this + areafact(i,j,ng) = -areafac_c(i,j) ! TC3b (group 6) From f3770d352e5ebe37cee6607e8fbcf1a0c6a4eace Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 29 Aug 2023 19:01:05 +0000 Subject: [PATCH 12/18] Fixed out of bounds areafac_ce and now use earea and narea --- .../cicedyn/dynamics/ice_transport_remap.F90 | 47 +++++++++---------- 1 file changed, 21 insertions(+), 26 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 61deb2487..f28a2218e 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -356,7 +356,7 @@ subroutine horizontal_remap (dt, ntrace, & use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_remap use ice_blocks, only: block, get_block, nghost use ice_grid, only: HTE, HTN, dxu, dyu, & - tarear, hm, & + earea, narea, tarear, hm, & xav, yav, xxav, yyav ! xyav, xxxav, xxyav, xyyav, yyyav use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound @@ -727,6 +727,7 @@ subroutine horizontal_remap (dt, ntrace, & indxing(:,:), indxjng(:,:), & dpx (:,:,iblk), dpy(:,:,iblk), & dxu (:,:,iblk), dyu(:,:,iblk), & + earea (:,:,iblk), narea (:,:,iblk), & xp (:,:,:,:), yp (:,:,:,:), & iflux, jflux, & triarea, edgearea_e(:,:)) @@ -786,6 +787,7 @@ subroutine horizontal_remap (dt, ntrace, & indxing(:,:), indxjng(:,:), & dpx (:,:,iblk), dpy (:,:,iblk), & dxu (:,:,iblk), dyu (:,:,iblk), & + earea (:,:,iblk), narea (:,:,iblk), & xp (:,:,:,:), yp(:,:,:,:), & iflux, jflux, & triarea, edgearea_n(:,:)) @@ -1705,6 +1707,7 @@ subroutine locate_triangles (nx_block, ny_block, & indxi, indxj, & dpx, dpy, & dxu, dyu, & + earea, narea, & xp, yp, & iflux, jflux, & triarea, edgearea) @@ -1721,7 +1724,9 @@ subroutine locate_triangles (nx_block, ny_block, & dpx , & ! x coordinates of departure points at cell corners dpy , & ! y coordinates of departure points at cell corners dxu , & ! E-W dimension of U-cell (m) - dyu ! N-S dimension of U-cell (m) + dyu , & ! N-S dimension of U-cell (m) + earea , & ! area of E-cell + narea ! area of N-cell real (kind=dbl_kind), dimension (nx_block,ny_block,0:nvert,ngroups), intent(out) :: & xp, yp ! coordinates of triangle vertices @@ -1771,12 +1776,10 @@ subroutine locate_triangles (nx_block, ny_block, & real (kind=dbl_kind), dimension(nx_block,ny_block) :: & dx, dy , & ! scaled departure points - areafac_c , & ! (areafac_l + areafac_r)/2 - areafac_ce ! areafac_c on other edge + areafac_c , & ! earea or narea + areafac_ce ! areafac_c on other edge (narea or earea) real (kind=dbl_kind) :: & - areafac_l , & ! area scale factor at left corner - areafac_r , & ! area scale factor at right corner xcl, ycl , & ! coordinates of left corner point ! (relative to midpoint of edge) xdl, ydl , & ! left departure point @@ -1923,21 +1926,17 @@ subroutine locate_triangles (nx_block, ny_block, & ! area scale factor - do j = jb, je - do i = ib, ie - areafac_l = dxu(i-1,j)*dyu(i-1,j) - areafac_r = dxu(i ,j)*dyu(i ,j) - areafac_c(i,j) = p5*( areafac_l + areafac_r ) + do j = jlo-1, jhi + do i = ilo, ihi + areafac_c(i,j) = narea(i,j) enddo enddo ! area scale factor for other edge (east) - do j = jb, je+1 - do i = ib-1, ie - areafac_l = dxu(i,j )*dyu(i,j ) - areafac_r = dxu(i,j-1)*dyu(i,j-1) - areafac_ce(i,j) = p5*( areafac_l + areafac_r ) + do j = jlo-1, jhi+1 + do i = ilo-1, ihi + areafac_ce(i,j) = earea(i,j) enddo enddo @@ -1978,21 +1977,17 @@ subroutine locate_triangles (nx_block, ny_block, & ! area scale factors - do j = jb, je - do i = ib, ie - areafac_l = dxu(i,j )*dyu(i,j ) - areafac_r = dxu(i,j-1)*dyu(i,j-1) - areafac_c(i,j) = p5*( areafac_l + areafac_r ) + do j = jlo, jhi + do i = ilo-1, ihi + areafac_c(i,j) = earea(i,j) enddo enddo ! area scale factor for other edge (north) - do j = jb-1, je - do i = ib, ie+1 - areafac_l = dxu(i-1,j)*dyu(i-1,j) - areafac_r = dxu(i ,j)*dyu(i ,j) - areafac_ce(i,j) = p5*( areafac_l + areafac_r ) + do j = jlo-1, jhi + do i = ilo-1, ihi+1 + areafac_ce(i,j) = narea(i,j) enddo enddo From ad89a2acbe233c35311ef126123928e225ae2888 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 29 Aug 2023 19:54:18 +0000 Subject: [PATCH 13/18] Replaced ib,ie,jb,je in locate_triangle using ilo,ihi,jlo,jhi --- .../cicedyn/dynamics/ice_transport_remap.F90 | 39 ++++++++----------- 1 file changed, 16 insertions(+), 23 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index f28a2218e..8ae227f24 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -1753,7 +1753,7 @@ subroutine locate_triangles (nx_block, ny_block, & integer (kind=int_kind) :: & i, j, ij, ic , & ! horizontal indices - ib, ie, jb, je , & ! limits for loops over edges + ib, jb , & ! limits for loops for bugcheck ng, nv , & ! triangle indices ishift , jshift , & ! differences between neighbor cells ishift_tl, jshift_tl , & ! i,j indices of TL cell relative to edge @@ -1891,13 +1891,6 @@ subroutine locate_triangles (nx_block, ny_block, & if (trim(edge) == 'north') then - ! loop size - - ib = ilo - ie = ihi - jb = jlo - nghost ! lowest j index is a ghost cell - je = jhi - ! index shifts for neighbor cells ishift_tl = -1 @@ -1942,13 +1935,6 @@ subroutine locate_triangles (nx_block, ny_block, & else ! east edge - ! loop size - - ib = ilo - nghost ! lowest i index is a ghost cell - ie = ihi - jb = jlo - je = jhi - ! index shifts for neighbor cells ishift_tl = 1 @@ -1999,8 +1985,8 @@ subroutine locate_triangles (nx_block, ny_block, & icellsd = 0 if (trim(edge) == 'north') then - do j = jb, je - do i = ib, ie + do j = jlo-1, jhi + do i = ilo, ihi if (dpx(i-1,j)/=c0 .or. dpy(i-1,j)/=c0 & .or. & dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then @@ -2011,8 +1997,8 @@ subroutine locate_triangles (nx_block, ny_block, & enddo enddo else ! east edge - do j = jb, je - do i = ib, ie + do j = jlo, jhi + do i = ilo-1, ihi if (dpx(i,j-1)/=c0 .or. dpy(i,j-1)/=c0 & .or. & dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then @@ -2028,8 +2014,8 @@ subroutine locate_triangles (nx_block, ny_block, & ! Scale the departure points !------------------------------------------------------------------- - do j = 1, je - do i = 1, ie + do j = 1, jhi + do i = 1, ihi dx(i,j) = dpx(i,j) / dxu(i,j) dy(i,j) = dpy(i,j) / dyu(i,j) enddo @@ -3061,10 +3047,17 @@ subroutine locate_triangles (nx_block, ny_block, & endif if (bugcheck) then + if (trim(edge) == 'north') then + ib = ilo + jb = jlo-1 + else ! east edge + ib = ilo-1 + jb = jlo + endif do ng = 1, ngroups do nv = 1, nvert - do j = jb, je - do i = ib, ie + do j = jb, jhi + do i = ib, ihi if (abs(triarea(i,j,ng)) > puny) then if (abs(xp(i,j,nv,ng)) > p5+puny) then write(nu_diag,*) '' From 0b9334be1b2a8f05de7a19a74e4a27fde9360c17 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Fri, 8 Sep 2023 16:33:36 +0000 Subject: [PATCH 14/18] Modified areafac for TL1, BL2, TR1 and BR2 for area flux consistency --- .../cicedyn/dynamics/ice_transport_remap.F90 | 24 +++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 8ae227f24..0e7d405b6 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -1762,6 +1762,8 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_br, jshift_br , & ! i,j indices of BR cell relative to edge ishift_tc, jshift_tc , & ! i,j indices of TC cell relative to edge ishift_bc, jshift_bc , & ! i,j indices of BC cell relative to edge + is_l, js_l , & ! index shifts for TL1 and BL2 for tri area consistency + is_r, js_r , & ! index shifts for TR1 and BR2 for tri area consistency ise_tl, jse_tl , & ! i,j of TL other edge relative to edge ise_bl, jse_bl , & ! i,j of BL other edge relative to edge ise_tr, jse_tr , & ! i,j of TR other edge relative to edge @@ -1906,6 +1908,13 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_bc = 0 jshift_bc = 0 + ! index shifts for TL1, BL2, TR1 and BR2 for triangle area consistency + + is_l = -1 + js_l = 0 + is_r = +1 + js_r = 0 + ! index shifts for neighbor east edges ise_tl = -1 @@ -1950,6 +1959,13 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_bc = 0 jshift_bc = 0 + ! index shifts for TL1, BL2, TR1 and BR2 for triangle area consistency + + is_l = 0 + js_l = +1 + is_r = 0 + js_r = -1 + ! index shifts for neighbor north edges ise_tl = 1 @@ -2128,7 +2144,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_tl jflux (i,j,ng) = j + jshift_tl - areafact(i,j,ng) = areafac_ce(i+ise_tl,j+jse_tl) + areafact(i,j,ng) = areafac_c(i+is_l,j+js_l) ! BL1 (group 3) @@ -2169,7 +2185,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydl iflux (i,j,ng) = i + ishift_bl jflux (i,j,ng) = j + jshift_bl - areafact(i,j,ng) = -areafac_ce(i+ise_bl,j+jse_bl) + areafact(i,j,ng) = -areafac_c(i+is_l,j+js_l) endif ! TL and BL triangles @@ -2221,7 +2237,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = ydr iflux (i,j,ng) = i + ishift_tr jflux (i,j,ng) = j + jshift_tr - areafact(i,j,ng) = areafac_ce(i+ise_tr,j+jse_tr) + areafact(i,j,ng) = areafac_c(i+is_r,j+js_r) ! BR1 (group 3) @@ -2262,7 +2278,7 @@ subroutine locate_triangles (nx_block, ny_block, & yp (i,j,3,ng) = yic iflux (i,j,ng) = i + ishift_br jflux (i,j,ng) = j + jshift_br - areafact(i,j,ng) = -areafac_ce(i+ise_br,j+jse_br) + areafact(i,j,ng) = -areafac_c(i+is_r,j+js_r) endif ! TR and BR triangles From ec88a5516e2760bddd6899126769a100cf46d9d4 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Fri, 8 Sep 2023 19:52:52 +0000 Subject: [PATCH 15/18] Cosmetic changes --- cicecore/cicedyn/dynamics/ice_transport_remap.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 0e7d405b6..00ed41ec7 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -1762,8 +1762,8 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_br, jshift_br , & ! i,j indices of BR cell relative to edge ishift_tc, jshift_tc , & ! i,j indices of TC cell relative to edge ishift_bc, jshift_bc , & ! i,j indices of BC cell relative to edge - is_l, js_l , & ! index shifts for TL1 and BL2 for tri area consistency - is_r, js_r , & ! index shifts for TR1 and BR2 for tri area consistency + is_l, js_l , & ! i,j shifts for TL1, BL2 for area consistency + is_r, js_r , & ! i,j shifts for TR1, BR2 for area consistency ise_tl, jse_tl , & ! i,j of TL other edge relative to edge ise_bl, jse_bl , & ! i,j of BL other edge relative to edge ise_tr, jse_tr , & ! i,j of TR other edge relative to edge @@ -1908,11 +1908,11 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_bc = 0 jshift_bc = 0 - ! index shifts for TL1, BL2, TR1 and BR2 for triangle area consistency + ! index shifts for TL1, BL2, TR1 and BR2 for area consistency is_l = -1 js_l = 0 - is_r = +1 + is_r = 1 js_r = 0 ! index shifts for neighbor east edges @@ -1959,10 +1959,10 @@ subroutine locate_triangles (nx_block, ny_block, & ishift_bc = 0 jshift_bc = 0 - ! index shifts for TL1, BL2, TR1 and BR2 for triangle area consistency + ! index shifts for TL1, BL2, TR1 and BR2 for area consistency is_l = 0 - js_l = +1 + js_l = 1 is_r = 0 js_r = -1 From c56ca3ea41fe5ea58f587d93770b117d044aac24 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Mon, 11 Sep 2023 13:08:07 +0000 Subject: [PATCH 16/18] Added comment to explain latest change --- cicecore/cicedyn/dynamics/ice_transport_remap.F90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 00ed41ec7..6a1294a6e 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -2099,6 +2099,13 @@ subroutine locate_triangles (nx_block, ny_block, & !------------------------------------------------------------------- ! Locate triangles in TL cell (NW for north edge, NE for east edge) ! and BL cell (W for north edge, N for east edge). + ! + ! areafact_c or areafac_ce (areafact_c for the other edge) are used + ! (with shifted indices) to make sure that a flux area on one edge + ! is consistent with the analogous area on the other edge and to + ! ensure that areas add up when using l_fixed_area = T. See PR #849 + ! for details. + ! !------------------------------------------------------------------- if (yil > c0 .and. xdl < xcl .and. ydl >= c0) then From a571e82e02dd7ec6b7d27ab2f947f2ab05cdc13b Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 12 Sep 2023 12:46:54 +0000 Subject: [PATCH 17/18] Modification of bugcheck condition for l_fixed_area=T --- cicecore/cicedyn/dynamics/ice_transport_remap.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 6a1294a6e..f91f96a33 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -63,7 +63,7 @@ module ice_transport_remap ! if false, area flux is determined internally ! and is passed out - logical (kind=log_kind), parameter :: bugcheck = .true. + logical (kind=log_kind), parameter :: bugcheck = .false. !======================================================================= ! Here is some information about how the incremental remapping scheme @@ -3008,7 +3008,7 @@ subroutine locate_triangles (nx_block, ny_block, & do ij = 1, icellsd i = indxid(ij) j = indxjd(ij) - if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then + if ( abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j) .and. abs(edgearea(i,j)) > c0 ) then write(nu_diag,*) '' write(nu_diag,*) 'Areas do not add up: m, i, j, edge =', & my_task, i, j, trim(edge) From 261e9b671175e83de1d4c40b0f7bb4efc5e6776f Mon Sep 17 00:00:00 2001 From: apcraig Date: Thu, 14 Sep 2023 07:04:34 -0600 Subject: [PATCH 18/18] update areafac_c, areafac_ce in halo in dynamics --- .../cicedyn/dynamics/ice_transport_remap.F90 | 18 ++++++----- cicecore/cicedyn/infrastructure/ice_grid.F90 | 30 +++++++++++++------ 2 files changed, 31 insertions(+), 17 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index f91f96a33..b397b94b7 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -1927,17 +1927,18 @@ subroutine locate_triangles (nx_block, ny_block, & jse_br = 0 ! area scale factor + ! earea, narea valid on halo - do j = jlo-1, jhi - do i = ilo, ihi + do j = 1, ny_block + do i = 1, nx_block areafac_c(i,j) = narea(i,j) enddo enddo ! area scale factor for other edge (east) - do j = jlo-1, jhi+1 - do i = ilo-1, ihi + do j = 1, ny_block + do i = 1, nx_block areafac_ce(i,j) = earea(i,j) enddo enddo @@ -1978,17 +1979,18 @@ subroutine locate_triangles (nx_block, ny_block, & jse_br = -1 ! area scale factors + ! earea, narea valid on halo - do j = jlo, jhi - do i = ilo-1, ihi + do j = 1, ny_block + do i = 1, nx_block areafac_c(i,j) = earea(i,j) enddo enddo ! area scale factor for other edge (north) - do j = jlo-1, jhi - do i = ilo-1, ihi+1 + do j = 1, ny_block + do i = 1, nx_block areafac_ce(i,j) = narea(i,j) enddo enddo diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index 770ee9ed9..160e3cc64 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -82,14 +82,14 @@ module ice_grid dyE , & ! height of E-cell through the middle (m) HTE , & ! length of eastern edge of T-cell (m) HTN , & ! length of northern edge of T-cell (m) - tarea , & ! area of T-cell (m^2) - uarea , & ! area of U-cell (m^2) - narea , & ! area of N-cell (m^2) - earea , & ! area of E-cell (m^2) - tarear , & ! 1/tarea - uarear , & ! 1/uarea - narear , & ! 1/narea - earear , & ! 1/earea + tarea , & ! area of T-cell (m^2), valid in halo + uarea , & ! area of U-cell (m^2), valid in halo + narea , & ! area of N-cell (m^2), valid in halo + earea , & ! area of E-cell (m^2), valid in halo + tarear , & ! 1/tarea, valid in halo + uarear , & ! 1/uarea, valid in halo + narear , & ! 1/narea, valid in halo + earear , & ! 1/earea, valid in halo tarean , & ! area of NH T-cells tareas , & ! area of SH T-cells ULON , & ! longitude of velocity pts, NE corner of T pts (radians) @@ -101,7 +101,7 @@ module ice_grid ELON , & ! longitude of center of east face of T pts (radians) ELAT , & ! latitude of center of east face of T pts (radians) ANGLE , & ! for conversions between POP grid and lat/lon - ANGLET , & ! ANGLE converted to T-cells + ANGLET , & ! ANGLE converted to T-cells, valid in halo bathymetry , & ! ocean depth, for grounding keels and bergs (m) ocn_gridcell_frac ! only relevant for lat-lon grids ! gridcell value of [1 - (land fraction)] (T-cell) @@ -635,12 +635,24 @@ subroutine init_grid2 call ice_HaloUpdate (uarea, halo_info, & field_loc_NEcorner, field_type_scalar, & fillValue=c1, tripoleOnly=.true.) + call ice_HaloUpdate (narea, halo_info, & + field_loc_Nface, field_type_scalar, & + fillValue=c1, tripoleOnly=.true.) + call ice_HaloUpdate (earea, halo_info, & + field_loc_Eface, field_type_scalar, & + fillValue=c1, tripoleOnly=.true.) call ice_HaloUpdate (tarear, halo_info, & field_loc_center, field_type_scalar, & fillValue=c1, tripoleOnly=.true.) call ice_HaloUpdate (uarear, halo_info, & field_loc_NEcorner, field_type_scalar, & fillValue=c1, tripoleOnly=.true.) + call ice_HaloUpdate (narear, halo_info, & + field_loc_Nface, field_type_scalar, & + fillValue=c1, tripoleOnly=.true.) + call ice_HaloUpdate (earear, halo_info, & + field_loc_Eface, field_type_scalar, & + fillValue=c1, tripoleOnly=.true.) call ice_timer_stop(timer_bound)