From e59e2d694d487f2d4cec2aa93d4e18e101b17716 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Sun, 15 Dec 2024 16:42:40 -0500 Subject: [PATCH 1/6] MEKE: Split src and src_GM compute loops The MEKE GM computation of `src` and `src_GM`, a diagnostic array, were placed in a single loop. The similar RHS of each expression made it unfavorable to use FMAs on the `src` update. Older production runs depending on this FMA were seeing answer changes. This patch restores the FMA loop update of `src` by separating `src` and `src_GM` into separate loops. --- src/parameterizations/lateral/MOM_MEKE.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 4bac09b7cc..886ce720d5 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -448,6 +448,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do j=js,je ; do i=is,ie src(i,j) = src(i,j) - CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j) + enddo ; enddo + + do j=js,je ; do i=is,ie src_GM(i,j) = -CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j) enddo ; enddo endif From 0a7ca20b6a97b25a4e6618e6d13d936b47564072 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Mon, 16 Dec 2024 00:48:22 -0500 Subject: [PATCH 2/6] MEKE: Conditionally compute biharmonic FrCoeff This patch makes several adjustments to MOM_MEKE.F90 and MOM_hor_visc.F90 to ensure that the Laplacian and biharmonic friction coefficients are computed separately, and only if their respective terms are enabled. This resolves some subtle bugs where the default biharmonic value of -1 was applied to the Laplacian case, even when the biharmonic MEKE friction was disabled. --- src/parameterizations/lateral/MOM_MEKE.F90 | 57 +++++++++++++++---- .../lateral/MOM_hor_visc.F90 | 26 ++++++--- 2 files changed, 64 insertions(+), 19 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 886ce720d5..bf02248876 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -227,6 +227,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h Kh_v, & ! The meridional diffusivity that is actually used [L2 T-1 ~> m2 s-1]. baroHv, & ! Depth integrated accumulated meridional mass flux [R Z L2 ~> kg]. drag_vel_v ! A piston velocity associated with bottom drag at v-points [H T-1 ~> m s-1 or kg m-2 s-1] + real :: bh_coeff ! Biharmonic part of efficiency conversion in total MEKE [nondim] real :: Kh_here ! The local horizontal viscosity [L2 T-1 ~> m2 s-1] real :: Inv_Kh_max ! The inverse of the local horizontal viscosity [T L-2 ~> s m-2] real :: K4_here ! The local horizontal biharmonic viscosity [L4 T-1 ~> m4 s-1] @@ -412,24 +413,53 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do j=js,je ; do i=is,ie src(i,j) = CS%MEKE_BGsrc - src_adv(i,j) = 0. - src_mom_K4(i,j) = 0. - src_btm_drag(i,j) = 0. - src_GM(i,j) = 0. - src_mom_lp(i,j) = 0. - src_mom_bh(i,j) = 0. enddo ; enddo + ! Initialize diagnostics + ! TODO: Remove these after the damping is only computed if the diagnostic + ! is registered. + if (CS%id_src_adv > 0) src_adv(is:ie, js:je) = 0. + if (CS%id_src_GM > 0) src_GM(is:ie, js:je) = 0. + if (CS%id_src_mom_lp > 0) src_mom_lp(is:ie, js:je) = 0. + if (CS%id_src_mom_bh > 0) src_mom_bh(is:ie, js:je) = 0. + if (CS%id_src_mom_K4 > 0) src_mom_K4(is:ie, js:je) = 0. + if (CS%id_src_btm_drag > 0) src_btm_drag(is:ie, js:je) = 0. + if (allocated(MEKE%mom_src)) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - src(i,j) = src(i,j) - CS%MEKE_FrCoeff*I_mass(i,j)*MEKE%mom_src(i,j) & - - (CS%MEKE_bhFrCoeff-CS%MEKE_FrCoeff)*I_mass(i,j)*MEKE%mom_src_bh(i,j) - src_mom_lp(i,j) = - CS%MEKE_FrCoeff*I_mass(i,j)*(MEKE%mom_src(i,j)-MEKE%mom_src_bh(i,j)) - src_mom_bh(i,j) = - CS%MEKE_bhFrCoeff*I_mass(i,j)*MEKE%mom_src_bh(i,j) + src(i,j) = src(i,j) - CS%MEKE_FrCoeff * I_mass(i,j) * MEKE%mom_src(i,j) enddo ; enddo endif + if (allocated(MEKE%mom_src_bh)) then + if (CS%MEKE_bhFrCoeff > 0. .and. CS%MEKE_FrCoeff > 0.) then + bh_coeff = CS%MEKE_bhFrCoeff - CS%MEKE_FrCoeff + else + bh_coeff = CS%MEKE_bhFrCoeff + endif + + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src(i,j) = src(i,j) - bh_coeff * I_mass(i,j) * MEKE%mom_src_bh(i,j) + enddo ; enddo + + if (CS%id_src_mom_lp > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_mom_lp(i,j) = -CS%MEKE_FrCoeff * I_mass(i,j) & + * (MEKE%mom_src(i,j) - MEKE%mom_src_bh(i,j)) + enddo ; enddo + endif + + if (CS%id_src_mom_bh > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_mom_bh(i,j) = -CS%MEKE_bhFrCoeff * I_mass(i,j) * MEKE%mom_src_bh(i,j) + enddo ; enddo + endif + endif + if (allocated(MEKE%GME_snk)) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie @@ -489,6 +519,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif + ! TODO: All of these diagnostics needs to be pulled out of the MEKE damping + ! loop and handled separately, and only if they are registered. + ! First stage of Strang splitting !$OMP parallel do default(shared) do j=js,je ; do i=is,ie @@ -1890,10 +1923,10 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) longname="Mesoscale Eddy Kinetic Energy", units="m2 s-2", conversion=US%L_T_to_m_s**2) if (MEKE_GMcoeff>=0.) allocate(MEKE%GM_src(isd:ied,jsd:jed), source=0.0) - if (MEKE_FrCoeff>=0. .or. MEKE_bhFrCoeff>=0. .or. MEKE_GMECoeff>=0.) then + if (MEKE_FrCoeff>=0. .or. MEKE_bhFrCoeff>=0. .or. MEKE_GMECoeff>=0.) & allocate(MEKE%mom_src(isd:ied,jsd:jed), source=0.0) + if (MEKE_bhFrCoeff >= 0.) & allocate(MEKE%mom_src_bh(isd:ied,jsd:jed), source=0.0) - endif if (MEKE_FrCoeff<0.) MEKE_FrCoeff = 0. if (MEKE_bhFrCoeff<0.) MEKE_bhFrCoeff = 0. if (MEKE_GMECoeff>=0.) allocate(MEKE%GME_snk(isd:ied,jsd:jed), source=0.0) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 92794c54e7..c2496e0d3a 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -2126,8 +2126,14 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (k==1) then do j=js,je ; do i=is,ie MEKE%mom_src(i,j) = 0. - MEKE%mom_src_bh(i,j) = 0. enddo ; enddo + + if (allocated(MEKE%mom_src_bh)) then + do j=js,je ; do i=is,ie + MEKE%mom_src_bh(i,j) = 0. + enddo ; enddo + endif + if (allocated(MEKE%GME_snk)) then do j=js,je ; do i=is,ie MEKE%GME_snk(i,j) = 0. @@ -2160,15 +2166,21 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + (FrictWork(i,j,k) - RoScl*FrictWork_bh(i,j,k)) - MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + & - (FrictWork_bh(i,j,k) - RoScl*FrictWork_bh(i,j,k)) + + if (allocated(MEKE%mom_src_bh)) & + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) & + + (FrictWork_bh(i,j,k) - RoScl * FrictWork_bh(i,j,k)) enddo ; enddo else + do j=js,je ; do i=is,ie + MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k) + enddo ; enddo - do j=js,je ; do i=is,ie - MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k) - MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + FrictWork_bh(i,j,k) - enddo ; enddo + if (allocated(MEKE%mom_src_bh)) then + do j=js,je ; do i=is,ie + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + FrictWork_bh(i,j,k) + enddo ; enddo + endif endif ! MEKE%backscatter_Ro_c if (CS%use_GME .and. allocated(MEKE%GME_snk)) then From 8a55ca4af4e1e24c54ceb9933784dc7ecfd7df94 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 17 Dec 2024 16:00:31 -0500 Subject: [PATCH 3/6] MEKE: Move flag test outside of FrictWork loop The if-test inside of the FrictWork loops are likely to impede performance. Even if the total work is reduced, they are likely to interrupt pipelines. When EY24_EBT_BS is disabled, they will clearly reduce performance. This patch moves those tests outside of the if-block and applies them separately. (Calculation would be slightly improved if the meaning of the flag were reversed, but I don't want to make additional changes.) --- .../lateral/MOM_hor_visc.F90 | 143 +++++++++--------- 1 file changed, 70 insertions(+), 73 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index c2496e0d3a..e4ea7f7ff9 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -480,12 +480,14 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 - visc_limit_h(:,:,:) = 0. - visc_limit_q(:,:,:) = 0. - visc_limit_h_flag(:,:,:) = 0. - visc_limit_q_flag(:,:,:) = 0. - visc_limit_h_frac(:,:,:) = 0. - visc_limit_q_frac(:,:,:) = 0. + if (CS%EY24_EBT_BS) then + visc_limit_h(:,:,:) = 0. + visc_limit_q(:,:,:) = 0. + visc_limit_h_flag(:,:,:) = 0. + visc_limit_q_flag(:,:,:) = 0. + visc_limit_h_frac(:,:,:) = 0. + visc_limit_q_frac(:,:,:) = 0. + endif m_leithy(:,:) = 0.0 ! Initialize @@ -1944,33 +1946,27 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (find_FrictWork) then if (CS%FrictWork_bug) then + ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion do j=js,je ; do i=is,ie - ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) - ! This is the old formulation that includes energy diffusion - if (visc_limit_h_flag(i,j,k) > 0) then - FrictWork(i,j,k) = 0 - else - FrictWork(i,j,k) = GV%H_to_RZ * ( & - ((str_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) & - - (str_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) & - + 0.25*(( (str_xy(I,J) * & - (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) & - + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) & - + (str_xy(I-1,J-1) * & - (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) & - + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) & - + ( (str_xy(I-1,J) * & - (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) & - + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) & - + (str_xy(I,J-1) * & - (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) & - + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) ) - endif + FrictWork(i,j,k) = GV%H_to_RZ * ( & + ((str_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) & + - (str_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) & + + 0.25*(( (str_xy(I,J) * & + (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) & + + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) & + + (str_xy(I-1,J-1) * & + (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) & + + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) & + + ( (str_xy(I-1,J) * & + (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) & + + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) & + + (str_xy(I,J-1) * & + (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) & + + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) ) enddo ; enddo - else ; do j=js,je ; do i=is,ie - if (visc_limit_h_flag(i,j,k) > 0) then - FrictWork(i,j,k) = 0 - else + else + do j=js,je ; do i=is,ie FrictWork(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & ((str_xx(i,j)*CS%dy2h(i,j) * ( & (uh(I,j,k)*G%dxCu(I,j)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) & @@ -1999,19 +1995,21 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) & - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) ) - endif - enddo ; enddo ; endif + enddo ; enddo + endif + + if (CS%EY24_EBT_BS) then + do j=js,je ; do i=is,ie + FrictWork(i,j,k) = (1. - visc_limit_h_flag(i,j,k)) * FrictWork(i,j,k) + enddo ; enddo + endif endif if (CS%id_FrictWork_bh>0 .or. CS%id_FrictWorkIntz_bh > 0 .or. allocated(MEKE%mom_src_bh)) then - if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie - ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) - ! This is the old formulation that includes energy diffusion - if (visc_limit_h_flag(i,j,k) > 0) then - FrictWork_bh(i,j,k) = 0 - else + if (CS%FrictWork_bug) then ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) ! This is the old formulation that includes energy diffusion !cyc + do j=js,je ; do i=is,ie FrictWork_bh(i,j,k) = GV%H_to_RZ * ( & (bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & - bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & @@ -2027,12 +2025,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, + bhstr_xy(I,J-1) * & ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) - endif - enddo ; enddo - else ; do j=js,je ; do i=is,ie - if (visc_limit_h_flag(i,j,k) > 0) then - FrictWork_bh(i,j,k) = 0 - else + enddo ; enddo + else + do j=js,je ; do i=is,ie ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) FrictWork_bh(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & ((bhstr_xx(i,j)*CS%dy2h(i,j) * ( & @@ -2061,11 +2056,15 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) & + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) & - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) ) - endif - enddo ; enddo ; endif - endif - + enddo ; enddo + endif + if (CS%EY24_EBT_BS) then + do j=js,je ; do i=is,ie + FrictWork_bh(i,j,k) = (1. - visc_limit_h_flag(i,j,k)) * FrictWork_bh(i,j,k) + enddo ; enddo + endif + endif if (CS%use_GME) then if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie @@ -2115,7 +2114,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) & + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) & - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) ) - enddo ; enddo ; endif endif @@ -2188,9 +2186,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, MEKE%GME_snk(i,j) = MEKE%GME_snk(i,j) + FrictWork_GME(i,j,k) enddo ; enddo endif - endif ! find_FrictWork and associated(mom_src) - enddo ! end of k loop ! Offer fields for diagnostic averaging. @@ -2219,16 +2215,15 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (CS%id_dudy_bt > 0) call post_data(CS%id_dudy_bt, dudy_bt, CS%diag) if (CS%id_dvdx_bt > 0) call post_data(CS%id_dvdx_bt, dvdx_bt, CS%diag) endif - if (CS%id_visc_limit_h>0) call post_data(CS%id_visc_limit_h, visc_limit_h, CS%diag) - if (CS%id_visc_limit_q>0) call post_data(CS%id_visc_limit_q, visc_limit_q, CS%diag) - if (CS%id_visc_limit_h_frac>0) call post_data(CS%id_visc_limit_h_frac, visc_limit_h_frac, CS%diag) - if (CS%id_visc_limit_q_frac>0) call post_data(CS%id_visc_limit_q_frac, visc_limit_q_frac, CS%diag) - if (CS%id_visc_limit_h_flag>0) call post_data(CS%id_visc_limit_h_flag, visc_limit_h_flag, CS%diag) - if (CS%id_visc_limit_q_flag>0) call post_data(CS%id_visc_limit_q_flag, visc_limit_q_flag, CS%diag) - if (CS%EY24_EBT_BS) then - if (CS%id_BS_coeff_h>0) call post_data(CS%id_BS_coeff_h, BS_coeff_h, CS%diag) - if (CS%id_BS_coeff_q>0) call post_data(CS%id_BS_coeff_q, BS_coeff_q, CS%diag) + if (CS%id_visc_limit_h>0) call post_data(CS%id_visc_limit_h, visc_limit_h, CS%diag) + if (CS%id_visc_limit_q>0) call post_data(CS%id_visc_limit_q, visc_limit_q, CS%diag) + if (CS%id_visc_limit_h_frac>0) call post_data(CS%id_visc_limit_h_frac, visc_limit_h_frac, CS%diag) + if (CS%id_visc_limit_q_frac>0) call post_data(CS%id_visc_limit_q_frac, visc_limit_q_frac, CS%diag) + if (CS%id_visc_limit_h_flag>0) call post_data(CS%id_visc_limit_h_flag, visc_limit_h_flag, CS%diag) + if (CS%id_visc_limit_q_flag>0) call post_data(CS%id_visc_limit_q_flag, visc_limit_q_flag, CS%diag) + if (CS%id_BS_coeff_h>0) call post_data(CS%id_BS_coeff_h, BS_coeff_h, CS%diag) + if (CS%id_BS_coeff_q>0) call post_data(CS%id_BS_coeff_q, BS_coeff_q, CS%diag) endif if (CS%debug) then @@ -3163,18 +3158,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) 'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T) CS%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, Time, & 'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim') - CS%id_visc_limit_h_flag = register_diag_field('ocean_model', 'visc_limit_h_flag', diag%axesTL, Time, & - 'Locations where the biharmonic viscosity reached the better_bound limiter at h points', 'nondim') - CS%id_visc_limit_q_flag = register_diag_field('ocean_model', 'visc_limit_q_flag', diag%axesBL, Time, & - 'Locations where the biharmonic viscosity reached the better_bound limiter at q points', 'nondim') - CS%id_visc_limit_h = register_diag_field('ocean_model', 'visc_limit_h', diag%axesTL, Time, & - 'Value of the biharmonic viscosity limiter at h points', 'nondim') - CS%id_visc_limit_q = register_diag_field('ocean_model', 'visc_limit_q', diag%axesBL, Time, & - 'Value of the biharmonic viscosity limiter at q points', 'nondim') - CS%id_visc_limit_h_frac = register_diag_field('ocean_model', 'visc_limit_h_frac', diag%axesTL, Time, & - 'Value of the biharmonic viscosity limiter at h points', 'nondim') - CS%id_visc_limit_q_frac = register_diag_field('ocean_model', 'visc_limit_q_frac', diag%axesBL, Time, & - 'Value of the biharmonic viscosity limiter at q points', 'nondim') + if (CS%EY24_EBT_BS) then + CS%id_visc_limit_h_flag = register_diag_field('ocean_model', 'visc_limit_h_flag', diag%axesTL, Time, & + 'Locations where the biharmonic viscosity reached the better_bound limiter at h points', 'nondim') + CS%id_visc_limit_q_flag = register_diag_field('ocean_model', 'visc_limit_q_flag', diag%axesBL, Time, & + 'Locations where the biharmonic viscosity reached the better_bound limiter at q points', 'nondim') + CS%id_visc_limit_h = register_diag_field('ocean_model', 'visc_limit_h', diag%axesTL, Time, & + 'Value of the biharmonic viscosity limiter at h points', 'nondim') + CS%id_visc_limit_q = register_diag_field('ocean_model', 'visc_limit_q', diag%axesBL, Time, & + 'Value of the biharmonic viscosity limiter at q points', 'nondim') + CS%id_visc_limit_h_frac = register_diag_field('ocean_model', 'visc_limit_h_frac', diag%axesTL, Time, & + 'Value of the biharmonic viscosity limiter at h points', 'nondim') + CS%id_visc_limit_q_frac = register_diag_field('ocean_model', 'visc_limit_q_frac', diag%axesBL, Time, & + 'Value of the biharmonic viscosity limiter at q points', 'nondim') + endif if (CS%id_grid_Re_Ah > 0) & ! Compute the smallest biharmonic viscosity capable of modifying the From 2788f40004b6bf755391198540984144c7938d08 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 18 Dec 2024 16:16:56 -0500 Subject: [PATCH 4/6] MEKE: Move damping diagnostics outside MEKE loop The damping MEKE loop also included updates to multiple diagnostics, even if they were not registered. This would presumably have a negative impact on performance. This patch moves each diagnostic into a separate loop. It also conditionally precomputes the damping and damp_rate parameters, which are now stored as 2d arrays rather than in-loop scalars. As before, the MEKE calculation is left unchanged in order to preserve bit reproducibility. --- src/parameterizations/lateral/MOM_MEKE.F90 | 235 ++++++++++++++++----- 1 file changed, 179 insertions(+), 56 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index bf02248876..4e874294ed 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -212,8 +212,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h barotrFac2, & ! Ratio of EKE_barotropic / EKE [nondim] bottomFac2, & ! Ratio of EKE_bottom / EKE [nondim] tmp, & ! Temporary variable for computation of diagnostic velocities [L T-1 ~> m s-1] - equilibrium_value ! The equilibrium value of MEKE to be calculated at each - ! time step [L2 T-2 ~> m2 s-2] + equilibrium_value, & ! The equilibrium value of MEKE to be calculated at + ! each time step [L2 T-2 ~> m2 s-2] + damp_rate, & ! The MEKE damping rate [T-1 ~> s-1] + damping ! The net damping of a field after sdt_damp [nondim] real, dimension(SZIB_(G),SZJ_(G)) :: & MEKE_uflux, & ! The zonal advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m2 s-3]. @@ -238,9 +240,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate) real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split). real :: damp_step ! Size of damping timestep relative to sdt [nondim] - real :: damp_rate ! The MEKE damping rate [T-1 ~> s-1]. - real :: damping ! The net damping of a field after sdt_damp [nondim] logical :: use_drag_rate ! Flag to indicate drag_rate is finite + logical :: any_damping_diags_s1 ! True if any damped diagnostics are enabled in first stage + logical :: any_damping_diags ! True if any damped diagnostics are enabled integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features ! needed for the machine learning inference, with different @@ -416,8 +418,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo ! Initialize diagnostics - ! TODO: Remove these after the damping is only computed if the diagnostic - ! is registered. if (CS%id_src_adv > 0) src_adv(is:ie, js:je) = 0. if (CS%id_src_GM > 0) src_GM(is:ie, js:je) = 0. if (CS%id_src_mom_lp > 0) src_mom_lp(is:ie, js:je) = 0. @@ -425,7 +425,22 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (CS%id_src_mom_K4 > 0) src_mom_K4(is:ie, js:je) = 0. if (CS%id_src_btm_drag > 0) src_btm_drag(is:ie, js:je) = 0. - if (allocated(MEKE%mom_src)) then + ! Identify any damped diagnostics in first stage of Strang splitting + any_damping_diags_s1 = any([ & + CS%id_src_GM > 0, & + CS%id_src_mom_lp > 0, & + CS%id_src_mom_bh > 0, & + CS%id_src_btm_drag > 0 & + ]) + + ! Identify any damped diagnostics + any_damping_diags = any([ & + any_damping_diags_s1, & + CS%id_src_adv > 0, & + CS%id_src_mom_K4 > 0 & + ]) + + if (CS%MEKE_FrCoeff > 0.) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie src(i,j) = src(i,j) - CS%MEKE_FrCoeff * I_mass(i,j) * MEKE%mom_src(i,j) @@ -519,36 +534,75 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif - ! TODO: All of these diagnostics needs to be pulled out of the MEKE damping - ! loop and handled separately, and only if they are registered. - ! First stage of Strang splitting + !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) damp_rate = 0. + damp_rate(i,j) = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + + if (MEKE%MEKE(i,j) < 0.) damp_rate(i,j) = 0. ! notice that the above line ensures a damping only if MEKE is positive, ! while leaving MEKE unchanged if it is negative + enddo ; enddo - damping = 1. / (1. + sdt_damp * damp_rate) + ! NOTE: MEKE%MEKE cannot use `damping` since we must preserve the existing + ! bit-reproducible solution. + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate(i,j)) + enddo ; enddo - ! NOTE: MEKE%MEKE should use `damping` but we must preserve the existing - ! expression for bit reproducibility - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate) - MEKE_decay(i,j) = damp_rate * G%mask2dT(i,j) + if (any_damping_diags_s1) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + damping(i,j) = 1. / (1. + sdt_damp * damp_rate(i,j)) + enddo ; enddo - src_GM(i,j) = src_GM(i,j) * damping - src_mom_lp(i,j) = src_mom_lp(i,j) * damping - src_mom_bh(i,j) = src_mom_bh(i,j) * damping + if (CS%id_decay > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE_decay(i,j) = damp_rate(i,j) * G%mask2dT(i,j) + enddo ; enddo + endif - src_btm_drag(i,j) = - MEKE_current(i,j) * ( & - damp_step * (damp_rate * damping) & - ) + if (CS%id_src_GM > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_GM(i,j) = src_GM(i,j) * damping(i,j) + enddo ; enddo + endif - ! Store the effective damping rate if sdt is split - if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) & - damp_rate_s1(i,j) = damp_rate * damping - enddo ; enddo + if (CS%id_src_mom_lp > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_mom_lp(i,j) = src_mom_lp(i,j) * damping(i,j) + enddo ; enddo + endif + + if (CS%id_src_mom_bh > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_mom_bh(i,j) = src_mom_bh(i,j) * damping(i,j) + enddo ; enddo + endif + + if (CS%id_src_btm_drag > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_btm_drag(i,j) = -MEKE_current(i,j) * ( & + damp_step * (damp_rate(i,j) * damping(i,j)) & + ) + enddo ; enddo + + ! Store the effective damping rate if sdt is split + if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + damp_rate_s1(i,j) = damp_rate(i,j) * damping(i,j) + enddo ; enddo + endif + endif + endif if (CS%kh_flux_enabled .or. CS%MEKE_K4 >= 0.0) then ! Update MEKE in the halos for lateral or bi-harmonic diffusion @@ -687,10 +741,16 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h MEKE%MEKE(i,j) = MEKE%MEKE(i,j) + (sdt*(G%IareaT(i,j)*I_mass(i,j))) * & ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) - src_adv(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * & - ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & - (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) enddo ; enddo + + if (CS%id_src_adv > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_adv(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * & + ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & + (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) + enddo ; enddo + endif endif ! MEKE_KH>0 ! Add on bi-harmonic tendency @@ -711,30 +771,80 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) enddo ; enddo endif + !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) damp_rate = 0. + damp_rate(i,j) = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + + if (MEKE%MEKE(i,j) < 0.) damp_rate(i,j) = 0. ! notice that the above line ensures a damping only if MEKE is positive, ! while leaving MEKE unchanged if it is negative + enddo ; enddo - damping = 1. / (1. + sdt_damp * damp_rate) + ! NOTE: MEKE%MEKE cannot use `damping` since we must preserve the + ! existing bit-reproducible solution. + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate(i,j)) + enddo ; enddo - ! NOTE: As above, MEKE%MEKE should use `damping` but we must preserve - ! the existing expression for bit reproducibility. - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*damp_rate) - MEKE_decay(i,j) = damp_rate*G%mask2dT(i,j) + if (any_damping_diags) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + damping(i,j) = 1. / (1. + sdt_damp * damp_rate(i,j)) + enddo ; enddo - src_GM(i,j) = src_GM(i,j) * damping - src_mom_lp(i,j) = src_mom_lp(i,j) * damping - src_mom_bh(i,j) = src_mom_bh(i,j) * damping - src_adv(i,j) = src_adv(i,j) * damping - src_mom_K4(i,j) = src_mom_K4(i,j) * damping + if (CS%id_decay > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE_decay(i,j) = damp_rate(i,j) * G%mask2dT(i,j) + enddo ; enddo + endif - src_btm_drag(i,j) = -MEKE_current(i,j) * ( & - damp_step * damping * (damp_rate + damp_rate_s1(i,j)) & - ) - enddo ; enddo + if (CS%id_src_GM > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_GM(i,j) = src_GM(i,j) * damping(i,j) + enddo ; enddo + endif + + if (CS%id_src_mom_lp > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_mom_lp(i,j) = src_mom_lp(i,j) * damping(i,j) + enddo ; enddo + endif + + if (CS%id_src_mom_bh > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_mom_bh(i,j) = src_mom_bh(i,j) * damping(i,j) + enddo ; enddo + endif + + if (CS%id_src_adv > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_adv(i,j) = src_adv(i,j) * damping(i,j) + enddo ; enddo + endif + + if (CS%id_src_mom_K4 > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_mom_K4(i,j) = src_mom_K4(i,j) * damping(i,j) + enddo ; enddo + endif + + if (CS%id_src_btm_drag > 0) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + src_btm_drag(i,j) = -MEKE_current(i,j) * (damp_step & + * ((damp_rate(i,j) + damp_rate_s1(i,j)) * damping(i,j)) & + ) + enddo ; enddo + endif + endif endif ! MEKE_KH>=0 if (CS%debug) then @@ -1522,20 +1632,33 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME if (.not. allocated(MEKE%MEKE)) CS%id_Ut = -1 CS%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, Time, & 'MEKE energy source', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) - !add diagnostics for the terms in the MEKE budget + CS%id_src_adv = register_diag_field('ocean_model', 'MEKE_src_adv', diag%axesT1, Time, & 'MEKE energy source from the horizontal advection of MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) - CS%id_src_mom_K4 = register_diag_field('ocean_model', 'MEKE_src_mom_K4', diag%axesT1, Time, & - 'MEKE energy source from the biharmonic of MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_btm_drag = register_diag_field('ocean_model', 'MEKE_src_btm_drag', diag%axesT1, Time, & 'MEKE energy source from the bottom drag acting on MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) - CS%id_src_GM = register_diag_field('ocean_model', 'MEKE_src_GM', diag%axesT1, Time, & - 'MEKE energy source from the thickness mixing (GM scheme)', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) - CS%id_src_mom_lp = register_diag_field('ocean_model', 'MEKE_src_mom_lp', diag%axesT1, Time, & - 'MEKE energy source from the Laplacian of resolved flows', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) - CS%id_src_mom_bh = register_diag_field('ocean_model', 'MEKE_src_mom_bh', diag%axesT1, Time, & - 'MEKE energy source from the biharmonic of resolved flows', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) - !end + + if (CS%MEKE_K4 >= 0.) & + CS%id_src_mom_K4 = register_diag_field('ocean_model', 'MEKE_src_mom_K4', & + diag%axesT1, Time, 'MEKE energy source from the biharmonic of MEKE', & + 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + + if (CS%MEKE_GMcoeff >= 0.) & + CS%id_src_GM = register_diag_field('ocean_model', 'MEKE_src_GM', & + diag%axesT1, Time, 'MEKE energy source from the thickness mixing (GM scheme)', & + 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + + if (CS%MEKE_FrCoeff >= 0.) & + CS%id_src_mom_lp = register_diag_field('ocean_model', 'MEKE_src_mom_lp', & + diag%axesT1, Time, 'MEKE energy source from the Laplacian of resolved flows', & + 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + + if (CS%MEKE_bhFrCoeff >= 0.) & + CS%id_src_mom_bh = register_diag_field('ocean_model', 'MEKE_src_mom_bh', & + diag%axesT1, Time, 'MEKE energy source from the biharmonic of resolved flows', & + 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, Time, & 'MEKE decay rate', 's-1', conversion=US%s_to_T) CS%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, Time, & From 5db4f28424796ab245dc175e2e445589569cecbc Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 19 Dec 2024 11:17:22 -0500 Subject: [PATCH 5/6] Diffusivity: Revert int_tide_CS to pointer The redefining of int_tide_CS control struct in set_diffusivity_init caused errors in debug-mode for Intel compilers. The issue appears to be an internal function that expects a pointer rather than the type. This patch reverts this back to a pointer. We can revisit this if there is a need to reduce reliance on pointers. --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index d15a364267..fc6d8380c5 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2201,7 +2201,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output. type(set_diffusivity_CS), pointer :: CS !< pointer set to point to the module control !! structure. - type(int_tide_CS), intent(in), target :: int_tide_CSp !< Internal tide control structure + type(int_tide_CS), pointer :: int_tide_CSp !< Internal tide control structure integer, intent(out) :: halo_TS !< The halo size of tracer points that must be !! valid for the calculations in set_diffusivity. logical, intent(out) :: double_diffuse !< This indicates whether some version From 79c922a60e41300ba995c87daf258c5b595f8aeb Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 2 Jan 2025 15:59:56 -0500 Subject: [PATCH 6/6] FMA rotation symmetric form of legacy FrictWork_bh This patch updates the expression for FrictWork_bh (biharmonic frictional work) when the FrictWork_bug flag is enabled. The new form is symmetric to rotations when FMA instructions are enabled. --- .../lateral/MOM_hor_visc.F90 | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index e4ea7f7ff9..a11b8a232e 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -2011,20 +2011,20 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, ! This is the old formulation that includes energy diffusion !cyc do j=js,je ; do i=is,ie FrictWork_bh(i,j,k) = GV%H_to_RZ * ( & - (bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & - - bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & - + 0.25*((bhstr_xy(I,J) * & - ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & - + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)) & - + bhstr_xy(I-1,J-1) * & - ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & - + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & - + (bhstr_xy(I-1,J) * & - ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & - + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & - + bhstr_xy(I,J-1) * & - ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & - + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) + ((bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) & + - (bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) & + + 0.25*(( (bhstr_xy(I,J) * & + (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) & + + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) & + + (bhstr_xy(I-1,J-1) * & + (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) & + + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) & + + ( (bhstr_xy(I-1,J) * & + (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) & + + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) & + + (bhstr_xy(I,J-1) * & + (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) & + + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) ) enddo ; enddo else do j=js,je ; do i=is,ie