From 53c293e448ac9ea314c13d79abbae210760863d5 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Tue, 18 Jul 2023 10:35:45 -0600 Subject: [PATCH 1/4] Add Leith+E This commit adds the 2D Leith+E closure, which uses a modified 2D Leith biharmonic viscosity paired with a harmonic backscatter. ('Modified' here is not used in the same sense as 'modified 2D Leith'; it just means that the biharmonic coefficient is modified to account for enstrophy backscatter.) Variables are often named 'leithy' to refer to Leith+E. The parameterization is controlled by three main entries in user_nl_mom: 1. USE_LEITHY = True 2. LEITH_CK = 1.0 3. LEITH_BI_CONST = 8.0 To use Leith+E you should have LAPLACIAN=True and BIHARMONIC=True. (It doesn't hurt to be explicit and also set LEITH_AH=False, along with any other viscous closures, but this is not required. If USE_LEITHY=True it will not use any of the other schemes. It does use the background value of the biharmonic coefficient as a minimum, but ignores the background harmonic value.) LEITH_CK is the fraction of energy dissipated by the biharmonic term that gets backscattered by the harmonic term (it's a target; the backscatter rate is not exact.) Recommended values between 0 and 1. LEITH_BI_CONST is Upsilon^6 where Upsilon is the ratio between the grid scale and the dissipation scale for enstrophy. Values should be greater than or equal to 1; 8 is a good place to start. The code is sensitive to the background value of Ah; specifically, if Ah is too large, the code is unstable. This is because the backscatter coefficient is proportional to Ah, and if Ah is large then you get large backscatter. If your code is unstable, consider reducing, e.g., `AH_VEL_SCALE`. --- .../lateral/MOM_hor_visc.F90 | 378 ++++++++++++++++-- 1 file changed, 344 insertions(+), 34 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index e6dd131a99..2af86488eb 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -68,6 +68,11 @@ module MOM_hor_visc logical :: use_beta_in_Leith !< If true, includes the beta term in the Leith viscosity logical :: Leith_Ah !< If true, use a biharmonic form of 2D Leith !! nonlinear eddy viscosity. AH is the background. + logical :: use_Leithy !< If true, use a biharmonic form of 2D Leith + !! nonlinear eddy viscosity with harmonic backscatter. + !! Ah is the background. Leithy = Leith+E + real :: c_K !< Fraction of energy dissipated by the biharmonic term + !! that gets backscattered in the Leith+E scheme. [nondim] logical :: use_QG_Leith_visc !< If true, use QG Leith nonlinear eddy viscosity. !! KH is the background value. logical :: bound_Coriolis !< If true & SMAGORINSKY_AH is used, the biharmonic @@ -145,10 +150,12 @@ module MOM_hor_visc n1n1_m_n2n2_q !< Factor n1**2-n2**2 in the anisotropic direction tensor at q-points [nondim] real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & - dx2h, & !< Pre-calculated dx^2 at h points [L2 ~> m2] - dy2h, & !< Pre-calculated dy^2 at h points [L2 ~> m2] - dx_dyT, & !< Pre-calculated dx/dy at h points [nondim] - dy_dxT !< Pre-calculated dy/dx at h points [nondim] + dx2h, & !< Pre-calculated dx^2 at h points [L2 ~> m2] + dy2h, & !< Pre-calculated dy^2 at h points [L2 ~> m2] + dx_dyT, & !< Pre-calculated dx/dy at h points [nondim] + dy_dxT, & !< Pre-calculated dy/dx at h points [nondim] + m_const_leithy, & !< Pre-calculated .5*sqrt(c_K)*max{dx,dy} [L ~> m] + m_leithy_max !< Pre-calculated 4./max(dx,dy)^2 at h points [L-2 ~> m-2] real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & dx2q, & !< Pre-calculated dx^2 at q points [L2 ~> m2] dy2q, & !< Pre-calculated dy^2 at q points [L2 ~> m2] @@ -257,18 +264,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Del2u, & ! The u-component of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1] h_u, & ! Thickness interpolated to u points [H ~> m or kg m-2]. vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] + vort_xy_dy_smooth, & ! y-derivative of smoothed vertical vorticity [L-1 T-1 ~> m-1 s-1] div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1] - ubtav ! zonal barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1] + ubtav, & ! zonal barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1] + u_smooth ! Zonal velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G)) :: & Del2v, & ! The v-component of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1] h_v, & ! Thickness interpolated to v points [H ~> m or kg m-2]. vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] + vort_xy_dx_smooth, & ! x-derivative of smoothed vertical vorticity [L-1 T-1 ~> m-1 s-1] div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1] - vbtav ! meridional barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1] + vbtav, & ! meridional barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1] + v_smooth ! Meridional velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G)) :: & dudx_bt, dvdy_bt, & ! components in the barotropic horizontal tension [T-1 ~> s-1] div_xx, & ! Estimate of horizontal divergence at h-points [T-1 ~> s-1] sh_xx, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] + sh_xx_smooth, & ! horizontal tension from smoothed velocity including metric terms [T-1 ~> s-1] sh_xx_bt, & ! barotropic horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] str_xx,& ! str_xx is the diagonal term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2], but ! at some points in the code it is not yet layer integrated, so is in [L2 T-2 ~> m2 s-2]. @@ -279,23 +291,28 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 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] + dudx_smooth, dvdy_smooth, & ! components in the horizontal tension from smoothed velocity [T-1 ~> s-1] GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim] - htot ! The total thickness of all layers [Z ~> m] + htot, & ! The total thickness of all layers [Z ~> m] + m_leithy ! Kh=m_leithy*Ah in Leith+E parameterization [L-2 ~> m-2] 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] + dvdx_smooth, dudy_smooth, & ! components in the shearing strain from smoothed velocity [T-1 ~> s-1] dDel2vdx, dDel2udy, & ! Components in the biharmonic equivalent of the shearing strain [L-2 T-1 ~> m-2 s-1] dvdx_bt, dudy_bt, & ! components in the barotropic shearing strain [T-1 ~> s-1] sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) including metric terms [T-1 ~> s-1] + sh_xy_smooth, & ! horizontal shearing strain from smoothed velocity including metric terms [T-1 ~> s-1] sh_xy_bt, & ! barotropic horizontal shearing strain (du/dy + dv/dx) inc. metric terms [T-1 ~> s-1] str_xy, & ! str_xy is the cross term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2], but ! at some points in the code it is not yet layer integrated, so is in [L2 T-2 ~> m2 s-2]. str_xy_GME, & ! smoothed cross term in the stress tensor from GME [L2 T-2 ~> m2 s-2] bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1] + vort_xy_smooth, & ! Vertical vorticity including metric terms, smoothed [T-1 ~> s-1] grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1] 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] @@ -331,6 +348,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, GME_coeff_h ! GME coefficient at h-points [L2 T-1 ~> m2 s-1] real :: AhSm ! Smagorinsky biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: AhLth ! 2D Leith biharmonic viscosity [L4 T-1 ~> m4 s-1] + real :: AhLthy ! 2D Leith+E biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: Shear_mag_bc ! Shear_mag value in backscatter [T-1 ~> s-1] real :: sh_xx_sq ! Square of tension (sh_xx) [T-2 ~> s-2] real :: sh_xy_sq ! Square of shearing strain (sh_xy) [T-2 ~> s-2] @@ -382,6 +400,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, 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] + vert_vort_mag_smooth, & ! magnitude of gradient of smoothed vertical vorticity (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 (h or q) [nondim] @@ -394,6 +413,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 + m_leithy(:,:) = 0. ! Initialize + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally apply_OBC = .true. @@ -546,7 +567,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & !$OMP dDel2vdx, dDel2udy, Del2vort_q, Del2vort_h, KE, & - !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff & + !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff, & + !$OMP dudx_smooth, dudy_smooth, dvdx_smooth, dvdy_smooth, & + !$OMP vort_xy_smooth, vort_xy_dx_smooth, vort_xy_dy_smooth, & + !$OMP sh_xx_smooth, sh_xy_smooth, u_smooth, v_smooth, & + !$OMP vert_vort_mag_smooth, m_leithy, AhLthy & !$OMP ) do k=1,nz @@ -575,6 +600,30 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j)) enddo ; enddo + if (CS%use_Leithy) then + ! Smooth the velocity. Right now it happens twice. In the future + ! one might make the number of smoothing cycles a user-specified parameter + u_smooth(:,:) = u(:,:,k) + v_smooth(:,:) = v(:,:,k) + call smooth_x9(CS, G, field_u=u_smooth,field_v=v_smooth) ! one call applies the filter twice + ! Calculate horizontal tension from smoothed velocity + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 + dudx_smooth(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u_smooth(I,j) - & + G%IdyCu(I-1,j) * u_smooth(I-1,j)) + dvdy_smooth(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v_smooth(i,J) - & + G%IdxCv(i,J-1) * v_smooth(i,J-1)) + sh_xx_smooth(i,j) = dudx_smooth(i,j) - dvdy_smooth(i,j) + enddo ; enddo + + ! Components for the shearing strain from smoothed velocity + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + dvdx_smooth(I,J) = CS%DY_dxBu(I,J) * & + (v_smooth(i+1,J)*G%IdyCv(i+1,J) - v_smooth(i,J)*G%IdyCv(i,J)) + dudy_smooth(I,J) = CS%DX_dyBu(I,J) * & + (u_smooth(I,j+1)*G%IdxCu(I,j+1) - u_smooth(I,j)*G%IdxCu(I,j)) + enddo ; enddo + end if ! use Leith+E + if (CS%id_normstress > 0) then do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 NoSt(i,j,k) = sh_xx(i,j) @@ -728,6 +777,20 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + if (CS%use_Leithy) then + ! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask). + ! dudy_smooth and dvdx_smooth do not (yet) include modifications at OBCs from above. + if (CS%no_slip) then + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + sh_xy_smooth(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_smooth(I,J) + dudy_smooth(I,J) ) + enddo ; enddo + else + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + sh_xy_smooth(I,J) = G%mask2dBu(I,J) * ( dvdx_smooth(I,J) + dudy_smooth(I,J) ) + enddo ; enddo + endif + endif ! use Leith+E + ! Evaluate Del2u = x.Div(Grad u) and Del2v = y.Div( Grad u) if (CS%biharmonic) then do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1 @@ -765,12 +828,24 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + if (CS%use_Leithy) then + if (CS%no_slip) then + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + vort_xy_smooth(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_smooth(I,J) - dudy_smooth(I,J) ) + enddo ; enddo + else + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + vort_xy_smooth(I,J) = G%mask2dBu(I,J) * ( dvdx_smooth(I,J) - dudy_smooth(I,J) ) + enddo ; enddo + endif + endif + ! Divergence do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 div_xx(i,j) = dudx(i,j) + dvdy(i,j) enddo ; enddo - if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then + if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then ! Vorticity gradient do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 @@ -783,6 +858,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1)) enddo ; enddo + if (CS%use_Leithy) then + ! Gradient of smoothed vorticity + do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 + DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) + vort_xy_dx_smooth(i,J) = DY_dxBu * & + (vort_xy_smooth(I,J) * G%IdyCu(I,j) - vort_xy_smooth(I-1,J) * G%IdyCu(I-1,j)) + enddo ; enddo + + do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+1 + DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) + vort_xy_dy_smooth(I,j) = DX_dyBu * & + (vort_xy_smooth(I,J) * G%IdxCv(i,J) - vort_xy_smooth(I,J-1) * G%IdxCv(i,J-1)) + enddo ; enddo + endif ! If Leithy + ! Laplacian of vorticity do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) @@ -865,6 +955,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, (0.5*(vort_xy_dy(I,j) + vort_xy_dy(I,j+1)))**2 ) enddo ; enddo + if (CS%use_Leithy) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + vert_vort_mag_smooth(i,j) = SQRT((0.5*(vort_xy_dx_smooth(i,J) + & + vort_xy_dx_smooth(i,J-1)))**2 + & + (0.5*(vort_xy_dy_smooth(I,j) + & + vort_xy_dy_smooth(I-1,j)))**2 ) + enddo ; enddo + endif ! Leithy + endif ! CS%Leith_Kh if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then @@ -890,6 +989,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif if (CS%Laplacian) then + ! Determine the Laplacian viscosity at h points, using the + ! largest value from several parameterizations. Also get + ! the Laplacian component of str_xx. + if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -904,9 +1007,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif - ! Determine the Laplacian viscosity at h points, using the - ! largest value from several parameterizations. - ! Static (pre-computed) background viscosity do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Kh(i,j) = CS%Kh_bg_xx(i,j) @@ -980,6 +1080,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + ! In Leith+E parameterization Kh is computed after Ah in the biharmonic loop. + ! The harmonic component of str_xx is added in the biharmonic loop. + if (CS%use_Leithy) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = 0. + enddo ; enddo + end if + if (CS%id_Kh_h>0 .or. CS%debug) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Kh_h(i,j,k) = Kh(i,j) @@ -1013,7 +1121,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = 0.0 enddo ; enddo - endif + endif ! Get Kh at h points and get Laplacian component of str_xx if (CS%anisotropic) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1026,12 +1134,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%biharmonic) then ! Determine the biharmonic viscosity at h points, using the - ! largest value from several parameterizations. + ! largest value from several parameterizations. Also get the + ! biharmonic component of str_xx. do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Ah(i,j) = CS%Ah_bg_xx(i,j) enddo ; enddo - if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah)) then + if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then if (CS%Smagorinsky_Ah) then if (CS%bound_Coriolis) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1057,12 +1166,49 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + if (CS%use_Leithy) then + ! Get m_leithy + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + 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) * inv_PI6 * abs(Del2vort_h) + if (AhLth .le. Ah(i,j)) then + m_leithy(i,j) = 0.0 + else + if ((CS%m_const_leithy(i,j)*vert_vort_mag(i,j)) .lt. abs(vort_xy_smooth(i,j))) then + m_leithy(i,j) = CS%c_K * (vert_vort_mag(i,j) / vort_xy_smooth(i,j))**2 + else + m_leithy(i,j) = CS%m_leithy_max(i,j) + endif + endif + enddo ; enddo + ! Smooth m_leithy + call smooth_x9(CS, G, field_h=m_leithy, zero_land=.true.) + ! Get Ah + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + & + (Del2vort_q(I-1,J) + Del2vort_q(I,J-1))) + AhLthy = CS%Biharm6_const_xx(i,j) * inv_PI6 * & + sqrt(max(0.,Del2vort_h**2 - m_leithy(i,j)*vert_vort_mag_smooth(i,j)**2)) + Ah(i,j) = max(Ah(i,j), AhLthy) + enddo ; enddo + ! Smooth Ah before applying upper bound + ! square, then smooth, then square root + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah_h(i,j,k) = Ah(i,j)**2 + enddo ; enddo + call smooth_x9(CS, G, field_h=Ah_h(:,:,k)) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = sqrt(Ah_h(i,j,k)) + enddo ; enddo + endif + if (CS%bound_Ah .and. .not. CS%better_bound_Ah) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xx(i,j)) enddo ; enddo endif - endif ! Smagorinsky_Ah or Leith_Ah + endif ! Smagorinsky_Ah or Leith_Ah or Leith+E if (use_MEKE_Au) then ! *Add* the MEKE contribution @@ -1096,6 +1242,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif + if (CS%use_Leithy) then + ! Compute Leith+E Kh after bounds have been applied to Ah + ! and after it has been smoothed. Kh = -m_leithy * Ah + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = -m_leithy(i,j) * Ah(i,j) + Kh_h(i,j,k) = Kh(i,j) + enddo ; enddo + endif + if (CS%id_grid_Re_Ah>0) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 KE = 0.125 * ((u(I,j,k) + u(I-1,j,k))**2 + (v(i,J,k) + v(i,J-1,k))**2) @@ -1111,10 +1266,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xx(i,j) = str_xx(i,j) + d_str + if (CS%use_Leithy) str_xx(i,j) = str_xx(i,j) - Kh(i,j) * sh_xx_smooth(i,j) + ! Keep a copy of the biharmonic contribution for backscatter parameterization bhstr_xx(i,j) = d_str * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo - endif + endif ! Get biharmonic coefficient at h points and biharmonic part of str_xx if (CS%biharmonic) then ! Gradient of Laplacian, for use in bi-harmonic term @@ -1203,6 +1360,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif if (CS%Laplacian) then + ! Determine the Laplacian viscosity at q points, using the + ! largest value from several parameterizations. Also get the + ! Laplacian component of str_xy. + if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then do J=js-1,Jeq ; do I=is-1,Ieq @@ -1217,9 +1378,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif - ! Determine the Laplacian viscosity at q points, using the - ! largest value from several parameterizations. - ! Static (pre-computed) background viscosity do J=js-1,Jeq ; do I=is-1,Ieq Kh(I,J) = CS%Kh_bg_xy(I,J) @@ -1286,6 +1444,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif + ! Leith+E doesn't recompute Kh at q points, it just interpolates it from h to q points + if (CS%use_Leithy) then + Kh(I,J) = Kh_h(i+1,j+1,k) + end if + if (CS%id_Kh_q>0 .or. CS%debug) & Kh_q(I,J,k) = Kh(I,J) @@ -1296,14 +1459,20 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, sh_xy_q(I,J,k) = sh_xy(I,J) enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq - str_xy(I,J) = -Kh(I,J) * sh_xy(I,J) - enddo ; enddo + if ( .not. CS%use_Leithy) then + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = -Kh(I,J) * sh_xy(I,J) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = -Kh(I,J) * sh_xy_smooth(I,J) + enddo ; enddo + endif else do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = 0. enddo ; enddo - endif + endif ! get harmonic coefficient Kh at q points and harmonic part of str_xy if (CS%anisotropic) then do J=js-1,Jeq ; do I=is-1,Ieq @@ -1316,7 +1485,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%biharmonic) then ! Determine the biharmonic viscosity at q points, using the - ! largest value from several parameterizations. + ! largest value from several parameterizations. Also get the + ! biharmonic component of str_xy. do J=js-1,Jeq ; do I=is-1,Ieq Ah(I,J) = CS%Ah_bg_xy(I,J) enddo ; enddo @@ -1380,6 +1550,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif + ! Leith+E doesn't recompute Ah at q points, it just interpolates it from h to q points + if (CS%use_Leithy) then + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(I,J) = Ah_h(i+1,j+1,k) + enddo ; enddo + end if + if (CS%id_Ah_q>0 .or. CS%debug) then do J=js-1,Jeq ; do I=is-1,Ieq Ah_q(I,J,k) = Ah(I,J) @@ -1395,7 +1572,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Keep a copy of the biharmonic contribution for backscatter parameterization bhstr_xy(I,J) = d_str * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) enddo ; enddo - endif + endif ! Get Ah at q points and biharmonic part of str_xy if (CS%use_GME) then ! The wider halo here is to permit one pass of smoothing without a halo update. @@ -1907,6 +2084,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "If true, use a biharmonic Leith nonlinear eddy "//& "viscosity.", default=.false., do_not_log=.not.CS%biharmonic) if (.not.CS%biharmonic) CS%Leith_Ah = .false. + call get_param(param_file, mdl, "USE_LEITHY", CS%use_Leithy, & + "If true, use a biharmonic Leith nonlinear eddy "//& + "viscosity together with a harmonic backscatter.", & + default=.false., do_not_log=.not.(CS%biharmonic .and. CS%Laplacian)) call get_param(param_file, mdl, "BOUND_AH", CS%bound_Ah, & "If true, the biharmonic coefficient is locally limited "//& "to be stable.", default=.true., do_not_log=.not.CS%biharmonic) @@ -1965,12 +2146,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "Coriolis acceleration. The default is set by MAXVEL.", & units="m s-1", default=maxvel*US%L_T_to_m_s, scale=US%m_s_to_L_T, & do_not_log=.not.(CS%Smagorinsky_Ah .and. CS%bound_Coriolis)) - call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, & "The nondimensional biharmonic Leith constant, "//& "typical values are thus far undetermined.", units="nondim", default=0.0, & - fail_if_missing=CS%Leith_Ah, do_not_log=.not.CS%Leith_Ah) - + fail_if_missing=(CS%Leith_Ah .or. CS%use_Leithy), & + do_not_log=.not.(CS%Leith_Ah .or. CS%use_Leithy)) call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", CS%use_land_mask, & "If true, use the land mask for the computation of thicknesses "//& "at velocity locations. This eliminates the dependence on arbitrary "//& @@ -2002,6 +2182,16 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "with the Gent and McWilliams parameterization.", default=.false.) call get_param(param_file, mdl, "SPLIT", split, & "Use the split time stepping if true.", default=.true., do_not_log=.true.) + if (CS%use_Leithy) then + if (.not.(CS%biharmonic .and. CS%Laplacian)) then + call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//& + "LAPLACIAN and BIHARMONIC must both be True when USE_LEITHY=True.") + endif + call get_param(param_file, mdl, "LEITHY_CK", CS%c_K, & + "Fraction of biharmonic dissipation that gets backscattered, "//& + "in Leith+E.", units="nondim", default=1.0) + endif + if (CS%use_GME .and. .not.split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & "cannot be used with SPLIT=False.") @@ -2120,9 +2310,13 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%Biharm_const2_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_const2_xy(:,:) = 0.0 endif endif - if (CS%Leith_Ah) then - ALLOC_(CS%biharm6_const_xx(isd:ied,jsd:jed)) ; CS%biharm6_const_xx(:,:) = 0.0 - ALLOC_(CS%biharm6_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm6_const_xy(:,:) = 0.0 + if ((CS%Leith_Ah) .or. (CS%use_Leithy)) then + ALLOC_(CS%biharm6_const_xx(isd:ied,jsd:jed)) ; CS%biharm6_const_xx(:,:) = 0.0 + ALLOC_(CS%biharm6_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm6_const_xy(:,:) = 0.0 + endif + if (CS%use_Leithy) then + ALLOC_(CS%m_const_leithy(isd:ied,jsd:jed)) ; CS%m_const_leithy(:,:) = 0.0 + ALLOC_(CS%m_leithy_max(isd:ied,jsd:jed)) ; CS%m_leithy_max(:,:) = 0.0 endif if (CS%Re_Ah > 0.0) then ALLOC_(CS%Re_Ah_const_xx(isd:ied,jsd:jed)); CS%Re_Ah_const_xx(:,:) = 0.0 @@ -2265,6 +2459,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) if (CS%Leith_Ah) then CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3) endif + if (CS%use_Leithy) then + CS%biharm6_const_xx(i,j) = Leith_bi_const * max(G%dxT(i,j),G%dyT(i,j))**6 + CS%m_const_leithy(i,j) = 0.5 * sqrt(CS%c_K) * max(G%dxT(i,j),G%dyT(i,j)) + CS%m_leithy_max(i,j) = 4. / max(G%dxT(i,j),G%dyT(i,j))**2 + endif CS%Ah_bg_xx(i,j) = MAX(Ah, Ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2)) if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xx(i,j) = grid_sp_h3 / CS%Re_Ah if (Ah_time_scale > 0.) CS%Ah_bg_xx(i,j) = & @@ -2287,7 +2486,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) (abs(G%CoriolisBu(I,J)) * BoundCorConst) endif endif - if (CS%Leith_Ah) then + if ((CS%Leith_Ah) .or. (CS%use_Leithy))then CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3) endif CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2)) @@ -2629,6 +2828,113 @@ subroutine smooth_GME(CS, G, GME_flux_h, GME_flux_q) enddo ! s-loop end subroutine smooth_GME +!> Apply a 9-point smoothing filter twice to reduce horizontal two-grid-point noise +!! Note that this subroutine does not conserve mass or angular momentum, so don't use it +!! in situations where you need conservation. Also can't apply it to Ah and Kh in the +!! horizontal_viscosity subroutine because they are not supposed to be halo-updated. +!! But you _can_ apply them to Kh_h and Ah_h. +subroutine smooth_x9(CS, G, field_h, field_u, field_v, field_q, zero_land) + type(hor_visc_CS), intent(in) :: CS !< Control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid + real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: field_h !< field to be smoothed + !! at h points + real, dimension(SZIB_(G),SZJ_(G)), optional, intent(inout) :: field_u !< field to be smoothed + !! at u points + real, dimension(SZI_(G),SZJB_(G)), optional, intent(inout) :: field_v !< field to be smoothed + !! at v points + real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: field_q !< field to be smoothed + !! at q points + logical, optional, intent(in) :: zero_land !< An optional argument + !! indicating whether to set values + !! on land to zero (.true.) or + !! whether to ignore land values + !! (.false. or not present) + ! local variables. It would be good to make the _original variables allocatable. + real, dimension(SZI_(G),SZJ_(G)) :: field_h_original + real, dimension(SZIB_(G),SZJ_(G)) :: field_u_original + real, dimension(SZI_(G),SZJB_(G)) :: field_v_original + real, dimension(SZIB_(G),SZJB_(G)) :: field_q_original + real, dimension(3,3) :: weights, local_weights ! averaging weights for smoothing, nondimensional + logical :: zero_land_val ! actual value of zero_land optional argument + integer :: i, j, s + 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 + + weights = reshape([1., 2., 1., 2., 4., 2., 1., 2., 1.],shape(weights))/16. + + if (present(zero_land)) then + zero_land_val = zero_land + else + zero_land_val = .false. + endif + + if (present(field_h)) then + call pass_var(field_h, G%Domain, halo=2) ! Halo size 2 ensures that you can smooth twice + do s=1,0,-1 + field_h_original(:,:) = field_h(:,:) + ! apply smoothing on field_h + do j=js-s,je+s ; do i=is-s,ie+s + ! skip land points + if (G%mask2dT(i,j)==0.) cycle + ! compute local weights + local_weights = weights*G%mask2dT(i-1:i+1,j-1:j+1) + if (zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16) + field_h(i,j) = sum(local_weights*field_h_original(i-1:i+1,j-1:j+1)) + enddo ; enddo + enddo + call pass_var(field_h, G%Domain) + endif + + if (present(field_u)) then + call pass_vector(field_u, field_v, G%Domain, halo=2) + do s=1,0,-1 + field_u_original(:,:) = field_u(:,:) + ! apply smoothing on field_u + do j=js-s,je+s ; do I=Isq-s,Ieq+s + ! skip land points + if (G%mask2dCu(I,j)==0.) cycle + ! compute local weights + local_weights = weights*G%mask2dCu(I-1:I+1,j-1:j+1) + if (zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16) + field_u(I,j) = sum(local_weights*field_u_original(I-1:I+1,j-1:j+1)) + enddo ; enddo + + field_v_original(:,:) = field_v(:,:) + ! apply smoothing on field_v + do J=Jsq-s,Jeq+s ; do i=is-s,ie+s + ! skip land points + if (G%mask2dCv(i,J)==0.) cycle + ! compute local weights + local_weights = weights*G%mask2dCv(i-1:i+1,J-1:J+1) + if (zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16) + field_v(i,J) = sum(local_weights*field_v_original(i-1:i+1,J-1:J+1)) + enddo ; enddo + enddo + call pass_vector(field_u, field_v, G%Domain) + endif + + if (present(field_q)) then + call pass_var(field_q, G%Domain, halo=2, position=CORNER) + do s=1,0,-1 + field_q_original(:,:) = field_q(:,:) + ! apply smoothing on field_q + do J=Jsq-s,Jeq+s ; do I=Isq-s,Ieq+s + ! skip land points + if (G%mask2dBu(I,J)==0.) cycle + ! compute local weights + local_weights = weights*G%mask2dBu(I-1:I+1,J-1:J+1) + if (zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16) + field_q(I,J) = sum(local_weights*field_q_original(I-1:I+1,J-1:J+1)) + enddo ; enddo + enddo + call pass_var(field_q, G%Domain, position=CORNER) + endif + +end subroutine smooth_x9 + !> Deallocates any variables allocated in hor_visc_init. subroutine hor_visc_end(CS) type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure @@ -2661,9 +2967,13 @@ subroutine hor_visc_end(CS) if (CS%Smagorinsky_Ah) then DEALLOC_(CS%Biharm_const_xx) ; DEALLOC_(CS%Biharm_const_xy) endif - if (CS%Leith_Ah) then + if ((CS%Leith_Ah) .or. (CS%use_Leithy)) then DEALLOC_(CS%Biharm6_const_xx) ; DEALLOC_(CS%Biharm6_const_xy) endif + if (CS%use_Leithy) then + DEALLOC_(CS%m_const_leithy) + DEALLOC_(CS%m_leithy_max) + endif if (CS%Re_Ah > 0.0) then DEALLOC_(CS%Re_Ah_const_xx) ; DEALLOC_(CS%Re_Ah_const_xy) endif From bd13ba6024213bdd83789d2fb4c2c6844ff4217a Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Tue, 18 Jul 2023 18:38:49 -0600 Subject: [PATCH 2/4] Background Ah This commit updates the code so that it uses the background Ah as a minimum. Previously, if `SMAGORINSKY_AH = True`, Leith+E would use the Smag value of Ah as the minimum, which is incorrect. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 2af86488eb..40a12c7598 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -1172,7 +1172,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, 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) * inv_PI6 * abs(Del2vort_h) - if (AhLth .le. Ah(i,j)) then + if (AhLth .le. CS%Ah_bg_xx(i,j)) then m_leithy(i,j) = 0.0 else if ((CS%m_const_leithy(i,j)*vert_vort_mag(i,j)) .lt. abs(vort_xy_smooth(i,j))) then @@ -1190,7 +1190,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, (Del2vort_q(I-1,J) + Del2vort_q(I,J-1))) AhLthy = CS%Biharm6_const_xx(i,j) * inv_PI6 * & sqrt(max(0.,Del2vort_h**2 - m_leithy(i,j)*vert_vort_mag_smooth(i,j)**2)) - Ah(i,j) = max(Ah(i,j), AhLthy) + Ah(i,j) = max(CS%Ah_bg_xx(i,j), AhLthy) enddo ; enddo ! Smooth Ah before applying upper bound ! square, then smooth, then square root From 754d5430e91f80b5e2da410735398493b7e50149 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Thu, 24 Aug 2023 07:48:05 -0600 Subject: [PATCH 3/4] Improve logging Removed `do_not_log` condition on `USE_LEITHY` --- src/parameterizations/lateral/MOM_hor_visc.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 40a12c7598..de92a006b7 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -2087,7 +2087,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) call get_param(param_file, mdl, "USE_LEITHY", CS%use_Leithy, & "If true, use a biharmonic Leith nonlinear eddy "//& "viscosity together with a harmonic backscatter.", & - default=.false., do_not_log=.not.(CS%biharmonic .and. CS%Laplacian)) + default=.false.) call get_param(param_file, mdl, "BOUND_AH", CS%bound_Ah, & "If true, the biharmonic coefficient is locally limited "//& "to be stable.", default=.true., do_not_log=.not.CS%biharmonic) From 4ebeabede9546582e6de80804ac6648faf24cb86 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Mon, 28 Aug 2023 12:47:04 -0600 Subject: [PATCH 4/4] Fix Leithy Logic Added one line to fix the fact that the code would only work as intended if either (i) writing out Ah_h, or (ii) in debug mode. Also swapped .le. and .lt. for <= and <. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index de92a006b7..40435359dc 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -1172,10 +1172,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, 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) * inv_PI6 * abs(Del2vort_h) - if (AhLth .le. CS%Ah_bg_xx(i,j)) then + if (AhLth <= CS%Ah_bg_xx(i,j)) then m_leithy(i,j) = 0.0 else - if ((CS%m_const_leithy(i,j)*vert_vort_mag(i,j)) .lt. abs(vort_xy_smooth(i,j))) then + if ((CS%m_const_leithy(i,j)*vert_vort_mag(i,j)) < abs(vort_xy_smooth(i,j))) then m_leithy(i,j) = CS%c_K * (vert_vort_mag(i,j) / vort_xy_smooth(i,j))**2 else m_leithy(i,j) = CS%m_leithy_max(i,j) @@ -1199,7 +1199,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo call smooth_x9(CS, G, field_h=Ah_h(:,:,k)) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - Ah(i,j) = sqrt(Ah_h(i,j,k)) + Ah_h(i,j,k) = sqrt(Ah_h(i,j,k)) + Ah(i,j) = Ah_h(i,j,k) enddo ; enddo endif