diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 78f48975e5..8e26014996 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -273,20 +273,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [R L2 T-3 ~> W m-2] grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] - Del2vort_h, & ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1] grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1] dudx, dvdy, & ! components in the horizontal tension [T-1 ~> s-1] - grad_vel_mag_h, & ! Magnitude of the velocity gradient tensor squared at h-points [T-2 ~> s-2] - grad_vel_mag_bt_h, & ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2] - grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2] GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim] - htot, & ! The total thickness of all layers [Z ~> m] - boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] - - real, allocatable, dimension(:,:,:) :: h_diffu ! h x diffu [H L T-2 ~> m2 s-2] - real, allocatable, dimension(:,:,:) :: h_diffv ! h x diffv [H L T-2 ~> m2 s-2] - real, allocatable, dimension(:,:,:) :: diffu_visc_rem ! diffu x visc_rem_u [L T-2 ~> m s-2] - real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [L T-2 ~> m s-2] + htot ! The total thickness of all layers [Z ~> m] + real :: Del2vort_h ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1] + real :: grad_vel_mag_bt_h ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2] + real :: boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1] @@ -304,12 +297,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1] Del2vort_q, & ! Laplacian of vorticity at q-points [L-2 T-1 ~> m-2 s-1] grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1] - grad_vel_mag_q, & ! Magnitude of the velocity gradient tensor squared at q-points [T-2 ~> s-2] hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] ! This form guarantees that hq/hu < 4. - grad_vel_mag_bt_q, & ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2] - GME_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim] - boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim] + GME_effic_q ! The filtered efficiency of the GME terms at q points [nondim] + real :: grad_vel_mag_bt_q ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2] + real :: boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim] real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: & Ah_q, & ! biharmonic viscosity at corner points [L4 T-1 ~> m4 s-1] @@ -347,8 +339,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! points; these are first interpolated to u or v velocity ! points where masks are applied [H ~> m or kg m-2]. real :: h_arith_q ! The arithmetic mean total thickness at q points [Z ~> m] - real :: h_harm_q ! The harmonic mean total thickness at q points [Z ~> m] - real :: I_hq ! The inverse of the arithmetic mean total thickness at q points [Z-1 ~> m-1] real :: I_GME_h0 ! The inverse of GME tapering scale [Z-1 ~> m-1] real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2] real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6] @@ -391,13 +381,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Fields evaluated on active layers, used for constructing 3D stress fields ! NOTE: The position of these declarations can impact performance, due to the ! very large number of stack arrays in this function. Move with caution! + ! NOTE: Several of these are declared with the memory extent of q-points, but the + ! same arrays are also used at h-points to reduce the memory footprint of this + ! module, so they should never be used in halo point or checksum calls. real, dimension(SZIB_(G),SZJB_(G)) :: & Ah, & ! biharmonic viscosity (h or q) [L4 T-1 ~> m4 s-1] - Kh, & ! Laplacian viscosity [L2 T-1 ~> m2 s-1] - Shear_mag, & ! magnitude of the shear [T-1 ~> s-1] - vert_vort_mag, & ! magnitude of the vertical vorticity gradient [L-1 T-1 ~> m-1 s-1] + Kh, & ! Laplacian viscosity (h or q) [L2 T-1 ~> m2 s-1] + Shear_mag, & ! magnitude of the shear (h or q) [T-1 ~> s-1] + vert_vort_mag, & ! magnitude of the vertical vorticity gradient (h or q) [L-1 T-1 ~> m-1 s-1] hrat_min, & ! h_min divided by the thickness at the stress point (h or q) [nondim] - visc_bound_rem ! fraction of overall viscous bounds that remain to be applied [nondim] + visc_bound_rem ! fraction of overall viscous bounds that remain to be applied (h or q) [nondim] real, dimension(SZIB_(G),SZJ_(G)) :: & hf_diffu_2d, & ! Depth sum of hf_diffu, hf_diffv [L T-2 ~> m s-2] @@ -456,14 +449,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%use_GME) then - do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 - boundary_mask_h(i,j) = (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1)) - enddo ; enddo - - do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1 - boundary_mask_q(I,J) = (G%mask2dCv(i,J) * G%mask2dCv(i+1,J) * G%mask2dCu(I,j) * G%mask2dCu(I,j+1)) - enddo ; enddo - ! initialize diag. array with zeros GME_coeff_h(:,:,:) = 0.0 GME_coeff_q(:,:,:) = 0.0 @@ -474,13 +459,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call barotropic_get_tav(BT, ubtav, vbtav, G, US) call pass_vector(ubtav, vbtav, G%Domain) + call pass_var(h, G%domain, halo=2) ! Calculate the barotropic horizontal tension - do j=Jsq-2,Jeq+2 ; do i=Isq-2,Ieq+2 + do J=js-2,je+2 ; do I=is-2,ie+2 dudx_bt(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * ubtav(I,j) - & G%IdyCu(I-1,j) * ubtav(I-1,j)) dvdy_bt(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * vbtav(i,J) - & G%IdxCv(i,J-1) * vbtav(i,J-1)) + enddo ; enddo + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j) enddo ; enddo @@ -492,12 +480,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, - ubtav(I,j)*G%IdxCu(I,j)) enddo ; enddo - ! post barotropic tension and strain - if (CS%id_dudx_bt > 0) call post_data(CS%id_dudx_bt, dudx_bt, CS%diag) - if (CS%id_dvdy_bt > 0) call post_data(CS%id_dvdy_bt, dvdy_bt, CS%diag) - 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) - if (CS%no_slip) then do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 sh_xy_bt(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_bt(I,J) + dudy_bt(I,J) ) @@ -508,20 +490,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif - do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 - grad_vel_mag_bt_h(i,j) = G%mask2dT(I,J) * boundary_mask_h(i,j) * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + & - (0.25*((dvdx_bt(I,J)+dvdx_bt(I-1,J-1))+(dvdx_bt(I,J-1)+dvdx_bt(I-1,J))))**2 + & - (0.25*((dudy_bt(I,J)+dudy_bt(I-1,J-1))+(dudy_bt(I,J-1)+dudy_bt(I-1,J))))**2) - enddo ; enddo - - do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1 - grad_vel_mag_bt_q(I,J) = G%mask2dBu(I,J) * boundary_mask_q(I,J) * (dvdx_bt(I,J)**2 + dudy_bt(I,J)**2 + & - (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1))+(dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + & - (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2) - enddo ; enddo - - call pass_var(h, G%domain, halo=2) - do j=js-2,je+2 ; do i=is-2,ie+2 htot(i,j) = 0.0 enddo ; enddo @@ -531,21 +499,29 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, I_GME_h0 = 1.0 / CS%GME_h0 do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 - if (grad_vel_mag_bt_h(i,j)>0) then - GME_effic_h(i,j) = CS%GME_efficiency * G%mask2dT(I,J) * & - (MIN(htot(i,j) * I_GME_h0, 1.0)**2) + boundary_mask_h = (G%mask2dCu(I,j) * G%mask2dCu(I-1,j)) * (G%mask2dCv(i,J) * G%mask2dCv(i,J-1)) + grad_vel_mag_bt_h = G%mask2dT(I,J) * boundary_mask_h * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + & + (0.25*((dvdx_bt(I,J)+dvdx_bt(I-1,J-1)) + (dvdx_bt(I,J-1)+dvdx_bt(I-1,J))))**2 + & + (0.25*((dudy_bt(I,J)+dudy_bt(I-1,J-1)) + (dudy_bt(I,J-1)+dudy_bt(I-1,J))))**2) + ! Probably the following test could be simplified to + ! if (boundary_mask_h * G%mask2dT(I,J) > 0.0) then + if (grad_vel_mag_bt_h > 0.0) then + GME_effic_h(i,j) = CS%GME_efficiency * G%mask2dT(I,J) * (MIN(htot(i,j) * I_GME_h0, 1.0)**2) else GME_effic_h(i,j) = 0.0 endif enddo ; enddo - do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1 - if (grad_vel_mag_bt_q(I,J)>0) then + do J=js-2,je+1 ; do I=is-2,ie+1 + boundary_mask_q = (G%mask2dCv(i,J) * G%mask2dCv(i+1,J)) * (G%mask2dCu(I,j) * G%mask2dCu(I,j+1)) + grad_vel_mag_bt_q = G%mask2dBu(I,J) * boundary_mask_q * (dvdx_bt(I,J)**2 + dudy_bt(I,J)**2 + & + (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1)) + (dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + & + (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1)) + (dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2) + ! Probably the following test could be simplified to + ! if (boundary_mask_q * G%mask2dBu(I,J) > 0.0) then + if (grad_vel_mag_bt_q > 0.0) then h_arith_q = 0.25 * ((htot(i,j) + htot(i+1,j+1)) + (htot(i+1,j) + htot(i,j+1))) - I_hq = 1.0 / h_arith_q - h_harm_q = 0.25 * h_arith_q * ((htot(i,j)*I_hq + htot(i+1,j+1)*I_hq) + & - (htot(i+1,j)*I_hq + htot(i,j+1)*I_hq)) - GME_effic_q(I,J) = CS%GME_efficiency * G%mask2dBu(I,J) * (MIN(h_harm_q * I_GME_h0, 1.0)**2) + GME_effic_q(I,J) = CS%GME_efficiency * G%mask2dBu(I,J) * (MIN(h_arith_q * I_GME_h0, 1.0)**2) else GME_effic_q(I,J) = 0.0 endif @@ -565,7 +541,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, & !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, & !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & - !$OMP use_MEKE_Ku, use_MEKE_Au, boundary_mask_h, boundary_mask_q, & + !$OMP use_MEKE_Ku, use_MEKE_Au, & !$OMP backscat_subround, GME_coeff_limiter, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI6, H0_GME, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & @@ -580,8 +556,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP vort_xy, vort_xy_dx, vort_xy_dy, div_xx, div_xx_dx, div_xx_dy, & !$OMP grad_div_mag_h, grad_div_mag_q, grad_vort_mag_h, grad_vort_mag_q, & !$OMP grad_vort, grad_vort_qg, grad_vort_mag_h_2d, grad_vort_mag_q_2d, & - !$OMP grad_vel_mag_h, grad_vel_mag_q, sh_xx_sq, sh_xy_sq, & - !$OMP grad_vel_mag_bt_h, grad_vel_mag_bt_q, grad_d2vel_mag_h, & + !$OMP sh_xx_sq, sh_xy_sq, & !$OMP meke_res_fn, Shear_mag, Shear_mag_bc, vert_vort_mag, h_min, hrat_min, visc_bound_rem, & !$OMP grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & @@ -831,9 +806,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Del2vort_q(I,J) = DY_dxBu * (vort_xy_dx(i+1,J) * G%IdyCv(i+1,J) - vort_xy_dx(i,J) * G%IdyCv(i,J)) + & DX_dyBu * (vort_xy_dy(I,j+1) * G%IdyCu(I,j+1) - vort_xy_dy(I,j) * G%IdyCu(I,j)) enddo ; enddo - do J=Jsq,Jeq+1 ; do I=Isq,Ieq+1 - Del2vort_h(i,j) = 0.25*(Del2vort_q(I,J) + Del2vort_q(I-1,J) + Del2vort_q(I,J-1) + Del2vort_q(I-1,J-1)) - enddo ; enddo if (CS%modified_Leith) then @@ -1017,8 +989,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (Kh(i,j) >= hrat_min(i,j) * CS%Kh_Max_xx(i,j)) then visc_bound_rem(i,j) = 0.0 Kh(i,j) = hrat_min(i,j) * CS%Kh_Max_xx(i,j) - else - ! ### NOTE: The denominator could be zero here - AJA ### + else ! if (Kh(i,j) > 0.0) then !### Change this to avoid a zero denominator. visc_bound_rem(i,j) = 1.0 - Kh(i,j) / (hrat_min(i,j) * CS%Kh_Max_xx(i,j)) endif enddo ; enddo @@ -1094,7 +1065,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%Leith_Ah) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h(i,j)) * inv_PI6 + Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + & + (Del2vort_q(I-1,J) + Del2vort_q(I,J-1))) + AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h) * inv_PI6 Ah(i,j) = max(Ah(i,j), AhLth) enddo ; enddo endif @@ -1215,7 +1188,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%better_bound_Kh) then do J=js-1,Jeq ; do I=is-1,Ieq - visc_bound_rem(i,j) = 1.0 + visc_bound_rem(I,J) = 1.0 enddo ; enddo endif endif @@ -1264,17 +1237,17 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Static (pre-computed) background viscosity do J=js-1,Jeq ; do I=is-1,Ieq - Kh(i,j) = CS%Kh_bg_xy(i,j) + Kh(I,J) = CS%Kh_bg_xy(I,J) enddo ; enddo if (CS%Smagorinsky_Kh) then if (CS%add_LES_viscosity) then do J=js-1,Jeq ; do I=is-1,Ieq - Kh(I,J) = Kh(I,J) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j) + Kh(I,J) = Kh(I,J) + CS%Laplac2_const_xy(I,J) * Shear_mag(I,J) enddo ; enddo else do J=js-1,Jeq ; do I=is-1,Ieq - Kh(I,J) = max(Kh(I,J), CS%Laplac2_const_xy(I,J) * Shear_mag(i,j) ) + Kh(I,J) = max(Kh(I,J), CS%Laplac2_const_xy(I,J) * Shear_mag(I,J) ) enddo ; enddo endif endif @@ -1282,7 +1255,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%Leith_Kh) then if (CS%add_LES_viscosity) then do J=js-1,Jeq ; do I=is-1,Ieq - Kh(I,J) = Kh(I,J) + CS%Laplac3_const_xx(i,j) * vert_vort_mag(I,J) * inv_PI3 ! Is this right? -AJA + Kh(I,J) = Kh(I,J) + CS%Laplac3_const_xy(I,J) * vert_vort_mag(I,J) * inv_PI3 ! Is this right? -AJA enddo ; enddo else do J=js-1,Jeq ; do I=is-1,Ieq @@ -1314,17 +1287,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, (MEKE%Ku(i+1,j) + MEKE%Ku(i,j+1)) ) * meke_res_fn endif - ! Older method of bounding for stability if (CS%anisotropic) & ! *Add* the shear component of anisotropic viscosity Kh(I,J) = Kh(I,J) + CS%Kh_aniso * CS%n1n2_q(I,J)**2 ! Newer method of bounding for stability if (CS%better_bound_Kh) then - if (Kh(i,j) >= hrat_min(I,J) * CS%Kh_Max_xy(I,J)) then - visc_bound_rem(i,j) = 0.0 - Kh(i,j) = hrat_min(I,J) * CS%Kh_Max_xy(I,J) - elseif (hrat_min(I,J)*CS%Kh_Max_xy(I,J)>0.) then + if (Kh(I,J) >= hrat_min(I,J) * CS%Kh_Max_xy(I,J)) then + visc_bound_rem(I,J) = 0.0 + Kh(I,J) = hrat_min(I,J) * CS%Kh_Max_xy(I,J) + elseif (hrat_min(I,J)*CS%Kh_Max_xy(I,J)>0.) then !### Change to elseif (Kh(I,J) > 0.0) then visc_bound_rem(I,J) = 1.0 - Kh(I,J) / (hrat_min(I,J) * CS%Kh_Max_xy(I,J)) endif endif @@ -1340,7 +1312,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq - str_xy(I,J) = -Kh(i,j) * sh_xy(I,J) + str_xy(I,J) = -Kh(I,J) * sh_xy(I,J) enddo ; enddo else do J=js-1,Jeq ; do I=is-1,Ieq @@ -1353,7 +1325,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Horizontal-tension averaged to q-points local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) ) ! *Add* the tension contribution to the xy-component of stress - str_xy(I,J) = str_xy(I,J) - CS%Kh_aniso * CS%n1n2_q(i,j) * CS%n1n1_m_n2n2_q(i,j) * local_strain + str_xy(I,J) = str_xy(I,J) - CS%Kh_aniso * CS%n1n2_q(I,J) * CS%n1n1_m_n2n2_q(I,J) * local_strain enddo ; enddo endif @@ -1361,7 +1333,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Determine the biharmonic viscosity at q points, using the ! largest value from several parameterizations. do J=js-1,Jeq ; do I=is-1,Ieq - Ah(i,j) = CS%Ah_bg_xy(I,J) + Ah(I,J) = CS%Ah_bg_xy(I,J) enddo ; enddo if (CS%Smagorinsky_Ah .or. CS%Leith_Ah) then @@ -1371,7 +1343,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, AhSm = Shear_mag(I,J) * (CS%Biharm_const_xy(I,J) & + CS%Biharm_const2_xy(I,J) * Shear_mag(I,J) & ) - Ah(i,j) = max(Ah(I,J), AhSm) + Ah(I,J) = max(Ah(I,J), AhSm) enddo ; enddo else do J=js-1,Jeq ; do I=is-1,Ieq @@ -1407,18 +1379,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%Re_Ah > 0.0) then do J=js-1,Jeq ; do I=is-1,Ieq KE = 0.125 * ((u(I,j,k) + u(I,j+1,k))**2 + (v(i,J,k) + v(i+1,J,k))**2) - Ah(I,J) = sqrt(KE) * CS%Re_Ah_const_xy(i,j) + Ah(I,J) = sqrt(KE) * CS%Re_Ah_const_xy(I,J) enddo ; enddo endif if (CS%better_bound_Ah) then if (CS%better_bound_Kh) then do J=js-1,Jeq ; do I=is-1,Ieq - Ah(I,J) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(I,J) * CS%Ah_Max_xy(I,J)) + Ah(I,J) = min(Ah(I,J), visc_bound_rem(I,J) * hrat_min(I,J) * CS%Ah_Max_xy(I,J)) enddo ; enddo else do J=js-1,Jeq ; do I=is-1,Ieq - Ah(I,J) = min(Ah(i,j), hrat_min(I,J) * CS%Ah_Max_xy(I,J)) + Ah(I,J) = min(Ah(I,J), hrat_min(I,J) * CS%Ah_Max_xy(I,J)) enddo ; enddo endif endif @@ -1441,6 +1413,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif if (CS%use_GME) then + ! The wider halo here is to permit one pass of smoothing without a halo update. do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 GME_coeff = GME_effic_h(i,j) * 0.25 * & ((KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) @@ -1450,7 +1423,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xx_GME(i,j) = GME_coeff * sh_xx_bt(i,j) enddo ; enddo - do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1 + ! The wider halo here is to permit one pass of smoothing without a halo update. + do J=js-2,je+1 ; do I=is-2,ie+1 GME_coeff = GME_effic_q(I,J) * 0.25 * & ((KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) GME_coeff = MIN(GME_coeff, CS%GME_limiter) @@ -1463,7 +1437,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call smooth_GME(CS, G, GME_flux_h=str_xx_GME) call smooth_GME(CS, G, GME_flux_q=str_xy_GME) - do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = (str_xx(i,j) + str_xx_GME(i,j)) * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo @@ -1496,10 +1470,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent. do j=js,je ; do I=Isq,Ieq - diffu(I,j,k) = ((G%IdyCu(I,j)*(CS%dy2h(i,j) *str_xx(i,j) - & - CS%dy2h(i+1,j)*str_xx(i+1,j)) + & - G%IdxCu(I,j)*(CS%dx2q(I,J-1)*str_xy(I,J-1) - & - CS%dx2q(I,J) *str_xy(I,J))) * & + diffu(I,j,k) = ((G%IdyCu(I,j)*(CS%dy2h(i,j)*str_xx(i,j) - CS%dy2h(i+1,j)*str_xx(i+1,j)) + & + G%IdxCu(I,j)*(CS%dx2q(I,J-1)*str_xy(I,J-1) - CS%dx2q(I,J)*str_xy(I,J))) * & G%IareaCu(I,j)) / (h_u(I,j) + h_neglect) enddo ; enddo @@ -1518,10 +1490,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Evaluate 1/h y.Div(h Grad u) or the biharmonic equivalent. do J=Jsq,Jeq ; do i=is,ie - diffv(i,J,k) = ((G%IdyCv(i,J)*(CS%dy2q(I-1,J)*str_xy(I-1,J) - & - CS%dy2q(I,J) *str_xy(I,J)) - & - G%IdxCv(i,J)*(CS%dx2h(i,j) *str_xx(i,j) - & - CS%dx2h(i,j+1)*str_xx(i,j+1))) * & + diffv(i,J,k) = ((G%IdyCv(i,J)*(CS%dy2q(I-1,J)*str_xy(I-1,J) - CS%dy2q(I,J)*str_xy(I,J)) - & + G%IdxCv(i,J)*(CS%dx2h(i,j)*str_xx(i,j) - CS%dx2h(i,j+1)*str_xx(i,j+1))) * & G%IareaCv(i,J)) / (h_v(i,J) + h_neglect) enddo ; enddo @@ -1542,41 +1512,40 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! 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 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) )) ) ) + (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 ; endif - if (CS%use_GME) then - do j=js,je ; do i=is,ie - ! Diagnose str_xx_GME*d_x u - str_yy_GME*d_y v + str_xy_GME*(d_y u + d_x v) - ! This is the old formulation that includes energy diffusion - FrictWork_GME(i,j,k) = GV%H_to_RZ * ( & + if (CS%use_GME) then ; do j=js,je ; do i=is,ie + ! Diagnose str_xx_GME*d_x u - str_yy_GME*d_y v + str_xy_GME*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion + FrictWork_GME(i,j,k) = GV%H_to_RZ * ( & (str_xx_GME(i,j)*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & - -str_xx_GME(i,j)*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & - +0.25*((str_xy_GME(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_GME(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_GME(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_GME(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) )) ) ) + - str_xx_GME(i,j)*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & + + 0.25*((str_xy_GME(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_GME(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_GME(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_GME(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 ; endif ! Make a similar calculation as for FrictWork above but accumulating into @@ -1621,18 +1590,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + GV%H_to_RZ * ( & ((str_xx(i,j)-RoScl*bhstr_xx(i,j))*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & -(str_xx(i,j)-RoScl*bhstr_xx(i,j))*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & - +0.25*(((str_xy(I,J)-RoScl*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) ) & - +(str_xy(I-1,J-1)-RoScl*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) )) & - +((str_xy(I-1,J)-RoScl*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) ) & - +(str_xy(I,J-1)-RoScl*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) )) ) ) + + 0.25*(((str_xy(I,J)-RoScl*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) ) & + + (str_xy(I-1,J-1)-RoScl*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)) ) & + + ((str_xy(I-1,J)-RoScl*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)) & + + (str_xy(I,J-1)-RoScl*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 endif ! MEKE%backscatter_Ro_c @@ -1656,19 +1625,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%id_diffu>0) call post_data(CS%id_diffu, diffu, CS%diag) if (CS%id_diffv>0) call post_data(CS%id_diffv, diffv, CS%diag) if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag) - if (CS%id_FrictWork_GME>0) call post_data(CS%id_FrictWork_GME, FrictWork_GME, CS%diag) if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag) if (CS%id_grid_Re_Ah>0) call post_data(CS%id_grid_Re_Ah, grid_Re_Ah, CS%diag) if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag) if (CS%id_vort_xy_q>0) call post_data(CS%id_vort_xy_q, vort_xy_q, CS%diag) - if (CS%id_sh_xx_h>0) call post_data(CS%id_sh_xx_h, sh_xx_h, CS%diag) - if (CS%id_sh_xy_q>0) call post_data(CS%id_sh_xy_q, sh_xy_q, CS%diag) + if (CS%id_sh_xx_h>0) call post_data(CS%id_sh_xx_h, sh_xx_h, CS%diag) + if (CS%id_sh_xy_q>0) call post_data(CS%id_sh_xy_q, sh_xy_q, CS%diag) if (CS%id_Ah_q>0) call post_data(CS%id_Ah_q, Ah_q, CS%diag) if (CS%id_Kh_h>0) call post_data(CS%id_Kh_h, Kh_h, CS%diag) if (CS%id_grid_Re_Kh>0) call post_data(CS%id_grid_Re_Kh, grid_Re_Kh, CS%diag) if (CS%id_Kh_q>0) call post_data(CS%id_Kh_q, Kh_q, CS%diag) - if (CS%id_GME_coeff_h > 0) call post_data(CS%id_GME_coeff_h, GME_coeff_h, CS%diag) - if (CS%id_GME_coeff_q > 0) call post_data(CS%id_GME_coeff_q, GME_coeff_q, CS%diag) + if (CS%use_GME) then ! post barotropic tension and strain + if (CS%id_GME_coeff_h > 0) call post_data(CS%id_GME_coeff_h, GME_coeff_h, CS%diag) + if (CS%id_GME_coeff_q > 0) call post_data(CS%id_GME_coeff_q, GME_coeff_q, CS%diag) + if (CS%id_FrictWork_GME>0) call post_data(CS%id_FrictWork_GME, FrictWork_GME, CS%diag) + if (CS%id_dudx_bt > 0) call post_data(CS%id_dudx_bt, dudx_bt, CS%diag) + if (CS%id_dvdy_bt > 0) call post_data(CS%id_dvdy_bt, dvdy_bt, CS%diag) + 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%debug) then if (CS%Laplacian) then @@ -2003,7 +1978,8 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) call get_param(param_file, mdl, "USE_KH_BG_2D", CS%use_Kh_bg_2d, & "If true, read a file containing 2-d background harmonic "//& "viscosities. The final viscosity is the maximum of the other "//& - "terms and this background value.", default=.false.) ! ###do_not_log=.not.CS%Laplacian? + "terms and this background value.", default=.false., do_not_log=.not.CS%Laplacian) + if (.not.CS%Laplacian) CS%use_Kh_bg_2d = .false. call get_param(param_file, mdl, "KH_BG_2D_BUG", CS%Kh_bg_2d_bug, & "If true, retain an answer-changing horizontal indexing bug in setting "//& "the corner-point viscosities when USE_KH_BG_2D=True. This is"//& @@ -2514,25 +2490,25 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) CS%min_grid_Kh = spacing(1.) * min_grid_sp_h2 * Idt endif if (CS%use_GME) then - CS%id_dudx_bt = register_diag_field('ocean_model', 'dudx_bt', diag%axesT1, Time, & + CS%id_dudx_bt = register_diag_field('ocean_model', 'dudx_bt', diag%axesT1, Time, & 'Zonal component of the barotropic shearing strain at h points', 's-1', & conversion=US%s_to_T) - CS%id_dudy_bt = register_diag_field('ocean_model', 'dudy_bt', diag%axesB1, Time, & + CS%id_dudy_bt = register_diag_field('ocean_model', 'dudy_bt', diag%axesB1, Time, & 'Zonal component of the barotropic shearing strain at q points', 's-1', & conversion=US%s_to_T) - CS%id_dvdy_bt = register_diag_field('ocean_model', 'dvdy_bt', diag%axesT1, Time, & + CS%id_dvdy_bt = register_diag_field('ocean_model', 'dvdy_bt', diag%axesT1, Time, & 'Meridional component of the barotropic shearing strain at h points', 's-1', & conversion=US%s_to_T) - CS%id_dvdx_bt = register_diag_field('ocean_model', 'dvdx_bt', diag%axesB1, Time, & + CS%id_dvdx_bt = register_diag_field('ocean_model', 'dvdx_bt', diag%axesB1, Time, & 'Meridional component of the barotropic shearing strain at q points', 's-1', & conversion=US%s_to_T) - CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, & + CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, & 'GME coefficient at h Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) - CS%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, Time, & + CS%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, Time, & 'GME coefficient at q Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) - CS%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,Time,& - 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', & - 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + CS%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,Time,& + 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) endif CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,& 'Integral work done by lateral friction terms. If GME is turned on, this '//& @@ -2577,57 +2553,62 @@ subroutine smooth_GME(CS, G, GME_flux_h, GME_flux_q) real, dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original real :: wc, ww, we, wn, ws ! averaging weights for smoothing - integer :: i, j, k, s + integer :: i, j, k, s, halosz + integer :: xh, xq ! The number of valid extra halo points for h and q points. integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + xh = 0 ; xq = 0 + do s=1,CS%num_smooth_gme if (present(GME_flux_h)) then - ! Update halos if needed - if (s >= 2) call pass_var(GME_flux_h, G%Domain) + if (xh < 0) then + ! Update halos if needed, but avoid doing so more often than is needed. + halosz = min(G%isc-G%isd, G%jsc-G%jsd, 2+CS%num_smooth_gme-s) + call pass_var(GME_flux_h, G%Domain, halo=halosz) + xh = halosz - 2 + endif GME_flux_h_original(:,:) = GME_flux_h(:,:) ! apply smoothing on GME - do j = Jsq, Jeq+1 - do i = Isq, Ieq+1 - ! skip land points - if (G%mask2dT(i,j)==0.) cycle - ! compute weights - ww = 0.125 * G%mask2dT(i-1,j) - we = 0.125 * G%mask2dT(i+1,j) - ws = 0.125 * G%mask2dT(i,j-1) - wn = 0.125 * G%mask2dT(i,j+1) - wc = 1.0 - ((ww+we)+(wn+ws)) - GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) & - + ((ww * GME_flux_h_original(i-1,j) & - + we * GME_flux_h_original(i+1,j)) & - + (ws * GME_flux_h_original(i,j-1) & - + wn * GME_flux_h_original(i,j+1))) - enddo - enddo + do j=Jsq-xh,Jeq+1+xh ; do i=Isq-xh,Ieq+1+xh + ! skip land points + if (G%mask2dT(i,j)==0.) cycle + ! compute weights + ww = 0.125 * G%mask2dT(i-1,j) + we = 0.125 * G%mask2dT(i+1,j) + ws = 0.125 * G%mask2dT(i,j-1) + wn = 0.125 * G%mask2dT(i,j+1) + wc = 1.0 - ((ww+we)+(wn+ws)) + GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) & + + ((ww * GME_flux_h_original(i-1,j) + we * GME_flux_h_original(i+1,j)) & + + (ws * GME_flux_h_original(i,j-1) + wn * GME_flux_h_original(i,j+1))) + enddo ; enddo + xh = xh - 1 endif if (present(GME_flux_q)) then - ! Update halos if needed - if (s >= 2) call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true.) + if (xq < 0) then + ! Update halos if needed, but avoid doing so more often than is needed. + halosz = min(G%isc-G%isd, G%jsc-G%jsd, 2+CS%num_smooth_gme-s) + call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true., halo=halosz) + xq = halosz - 2 + endif GME_flux_q_original(:,:) = GME_flux_q(:,:) ! apply smoothing on GME - do J = Jsq-1, Jeq - do I = Isq-1, Ieq - ! skip land points - if (G%mask2dBu(I,J)==0.) cycle - ! compute weights - ww = 0.125 * G%mask2dBu(I-1,J) - we = 0.125 * G%mask2dBu(I+1,J) - ws = 0.125 * G%mask2dBu(I,J-1) - wn = 0.125 * G%mask2dBu(I,J+1) - wc = 1.0 - ((ww+we)+(wn+ws)) - GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) & - + ((ww * GME_flux_q_original(I-1,J) & - + we * GME_flux_q_original(I+1,J)) & - + (ws * GME_flux_q_original(I,J-1) & - + wn * GME_flux_q_original(I,J+1))) - enddo - enddo + do J=js-1-xq,je+xq ; do I=is-1-xq,ie+xq + ! skip land points + if (G%mask2dBu(I,J)==0.) cycle + ! compute weights + ww = 0.125 * G%mask2dBu(I-1,J) + we = 0.125 * G%mask2dBu(I+1,J) + ws = 0.125 * G%mask2dBu(I,J-1) + wn = 0.125 * G%mask2dBu(I,J+1) + wc = 1.0 - ((ww+we)+(wn+ws)) + GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) & + + ((ww * GME_flux_q_original(I-1,J) + we * GME_flux_q_original(I+1,J)) & + + (ws * GME_flux_q_original(I,J-1) + wn * GME_flux_q_original(I,J+1))) + enddo ; enddo + xq = xq - 1 endif enddo ! s-loop end subroutine smooth_GME