diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 4208bc1642..f104353c1f 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -75,6 +75,8 @@ module MOM_hor_visc !! Default is False to maintain answers with legacy experiments !! but should be changed to True for new experiments. logical :: anisotropic !< If true, allow anisotropic component to the viscosity. + logical :: add_LES_viscosity!< If true, adds the viscosity from Smagorinsky and Leith to + !! the background viscosity instead of taking the maximum. real :: Kh_aniso !< The anisotropic viscosity [L2 T-1 ~> m2 s-1]. logical :: dynamic_aniso !< If true, the anisotropic viscosity is recomputed as a function !! of state. This is set depending on ANISOTROPIC_MODE. @@ -679,8 +681,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif - call pass_var(vort_xy, G%Domain, position=CORNER, complete=.true.) - ! Vorticity gradient do J=js-2,Jeq+1 ; do i=is-1,Ieq+1 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) @@ -692,23 +692,17 @@ 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 - call pass_vector(vort_xy_dy, vort_xy_dx, G%Domain) - if (CS%modified_Leith) then ! Divergence do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 - div_xx(i,j) = 0.5*((G%dyCu(I,j) * u(I,j,k) * (h(i+1,j,k)+h(i,j,k)) - & - G%dyCu(I-1,j) * u(I-1,j,k) * (h(i-1,j,k)+h(i,j,k)) ) + & - (G%dxCv(i,J) * v(i,J,k) * (h(i,j,k)+h(i,j+1,k)) - & - G%dxCv(i,J-1)*v(i,J-1,k)*(h(i,j,k)+h(i,j-1,k))))*G%IareaT(i,j) / & - (h(i,j,k) + GV%H_subroundoff) + div_xx(i,j) = dudx(i,j) + dvdy(i,j) enddo ; enddo ! Divergence gradient - do j=Jsq,Jeq+1 ; do I=Isq-1,Ieq+1 + do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+1 div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j)) enddo ; enddo - do J=Jsq-1,Jeq+1 ; do i=Isq,Ieq+1 + do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j)) enddo ; enddo @@ -717,7 +711,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + & (0.5 * (div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2) enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq + !do J=js-1,Jeq ; do I=is-1,Ieq + do j=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1 grad_div_mag_q(I,J) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + & (0.5 * (div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2) enddo ; enddo @@ -727,13 +722,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1 div_xx_dx(I,j) = 0.0 enddo ; enddo - do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2 + do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2 div_xx_dy(i,J) = 0.0 enddo ; enddo do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 grad_div_mag_h(i,j) = 0.0 enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq + do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 grad_div_mag_q(I,J) = 0.0 enddo ; enddo @@ -802,8 +797,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Determine the Laplacian viscosity at h points, using the ! largest value from several parameterizations. Kh = CS%Kh_bg_xx(i,j) ! Static (pre-computed) background viscosity - if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag ) - if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3) + if (CS%add_LES_viscosity) then + if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag + if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3 + else + if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag ) + if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3) + endif ! All viscosity contributions above are subject to resolution scaling if (rescale_Kh) Kh = VarMix%Res_fn_h(i,j) * Kh if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j) @@ -827,7 +827,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%id_Kh_h>0) .or. find_FrictWork .or. CS%debug) Kh_h(i,j,k) = Kh if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j) -! if (CS%debug) sh_xx_3d(i,j,k) = sh_xx(i,j) str_xx(i,j) = -Kh * sh_xx(i,j) else ! not Laplacian @@ -965,8 +964,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Determine the Laplacian viscosity at q points, using the ! largest value from several parameterizations. Kh = CS%Kh_bg_xy(i,j) ! Static (pre-computed) background viscosity - if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag ) - if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3) + if (CS%add_LES_viscosity) then + if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag + if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3 + else + if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag ) + if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3) + endif ! All viscosity contributions above are subject to resolution scaling if (rescale_Kh) Kh = VarMix%Res_fn_q(i,j) * Kh if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(i,j) @@ -993,7 +997,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%id_Kh_q>0 .or. CS%debug) Kh_q(I,J,k) = Kh if (CS%id_vort_xy_q>0) vort_xy_q(I,J,k) = vort_xy(I,J) -! if (CS%debug) sh_xy_3d(I,J,k) = sh_xy(I,J) str_xy(I,J) = -Kh * sh_xy(I,J) else ! not Laplacian @@ -1294,8 +1297,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%Laplacian) then call hchksum(Kh_h, "Kh_h", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T) call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T) -! call Bchksum(sh_xy_3d, "shear_xy", G%HI, haloshift=0, scale=US%s_to_T) -! call hchksum(sh_xx_3d, "shear_xx", G%HI, haloshift=0, scale=US%s_to_T) endif if (CS%biharmonic) call hchksum(Ah_h, "Ah_h", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T) if (CS%biharmonic) call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T) @@ -1504,6 +1505,9 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, & "If true, allow anistropic viscosity in the Laplacian "//& "horizontal viscosity.", default=.false.) + call get_param(param_file, mdl, "ADD_LES_VISCOSITY", CS%add_LES_viscosity, & + "If true, adds the viscosity from Smagorinsky and Leith to the "//& + "background viscosity instead of taking the maximum.", default=.false.) endif if (CS%anisotropic .or. get_all) then call get_param(param_file, mdl, "KH_ANISO", CS%Kh_aniso, &