diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 4bbd03a46a..14b0942009 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -48,7 +48,7 @@ module MOM_dynamics_split_RK2 use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS +use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS, hor_visc_vel_stencil use MOM_hor_visc, only : hor_visc_init, hor_visc_end use MOM_interface_heights, only : thickness_to_dz, find_col_avg_SpV use MOM_lateral_mixing_coeffs, only : VarMix_CS @@ -401,7 +401,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s logical :: showCallTree, sym integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz - integer :: cont_stencil, obc_stencil + integer :: cont_stencil, obc_stencil, vel_stencil is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -468,19 +468,20 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (associated(CS%OBC)) then if (CS%OBC%oblique_BCs_exist_globally) obc_stencil = 3 endif + vel_stencil = max(2, obc_stencil, hor_visc_vel_stencil(CS%hor_visc)) call cpu_clock_begin(id_clock_pass) call create_group_pass(CS%pass_eta, eta, G%Domain, halo=1) call create_group_pass(CS%pass_visc_rem, CS%visc_rem_u, CS%visc_rem_v, G%Domain, & To_All+SCALAR_PAIR, CGRID_NE, halo=max(1,cont_stencil)) call create_group_pass(CS%pass_uvp, up, vp, G%Domain, halo=max(1,cont_stencil)) call create_group_pass(CS%pass_hp_uv, hp, G%Domain, halo=2) - call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=max(2,obc_stencil)) - call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil)) + call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=vel_stencil) + call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=vel_stencil) call create_group_pass(CS%pass_uv, u, v, G%Domain, halo=max(2,cont_stencil)) call create_group_pass(CS%pass_h, h, G%Domain, halo=max(2,cont_stencil)) - call create_group_pass(CS%pass_av_uvh, u_av, v_av, G%Domain, halo=max(2,obc_stencil)) - call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil)) + call create_group_pass(CS%pass_av_uvh, u_av, v_av, G%Domain, halo=vel_stencil) + call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=vel_stencil) call cpu_clock_end(id_clock_pass) !--- end set up for group halo pass @@ -841,7 +842,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (CS%debug) then call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym) - call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) + call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=vel_stencil, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_MKS) ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US) call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index d59e6b3871..78c3971d87 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -31,7 +31,7 @@ module MOM_hor_visc #include -public horizontal_viscosity, hor_visc_init, hor_visc_end +public horizontal_viscosity, hor_visc_init, hor_visc_end, hor_visc_vel_stencil !> Control structure for horizontal viscosity type, public :: hor_visc_CS ; private @@ -1198,10 +1198,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, 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=js_Kh,je_Kh ; do i=is_Kh,ie_Kh + do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh AhSm = Shear_mag(i,j) * (CS%Biharm_const_xx(i,j) & - + CS%Biharm_const2_xx(i,j) * Shear_mag(i,j) & - ) + + CS%Biharm_const2_xx(i,j) * Shear_mag(i,j)) Ah(i,j) = max(Ah(i,j), AhSm) enddo ; enddo else @@ -1432,10 +1431,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Pass the velocity gradients and thickness to ZB2020 if (CS%use_ZB2020) then - call ZB2020_copy_gradient_and_thickness( & - sh_xx, sh_xy, vort_xy, & - hq, & - G, GV, CS%ZB2020, k) + call ZB2020_copy_gradient_and_thickness(sh_xx, sh_xy, vort_xy, hq, G, GV, CS%ZB2020, k) endif if (CS%Laplacian) then @@ -1575,8 +1571,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%bound_Coriolis) then do J=js-1,Jeq ; do I=is-1,Ieq AhSm = Shear_mag(I,J) * (CS%Biharm_const_xy(I,J) & - + CS%Biharm_const2_xy(I,J) * Shear_mag(I,J) & - ) + + CS%Biharm_const2_xy(I,J) * Shear_mag(I,J)) Ah(I,J) = max(Ah(I,J), AhSm) enddo ; enddo else @@ -1605,8 +1600,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! *Add* the MEKE contribution do J=js-1,Jeq ; do I=is-1,Ieq Ah(I,J) = Ah(I,J) + 0.25 * ( & - (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) & - ) + (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) ) enddo ; enddo endif @@ -1897,11 +1891,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%debug) then 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 hchksum(Kh_h, "Kh_h", G%HI, haloshift=1, scale=US%L_to_m**2*US%s_to_T) + call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**2*US%s_to_T) + endif + if (CS%biharmonic) then + call hchksum(Ah_h, "Ah_h", G%HI, haloshift=1, scale=US%L_to_m**4*US%s_to_T) + call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**4*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) endif if (CS%id_FrictWorkIntz > 0) then @@ -2403,14 +2399,31 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) 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 - ALLOC_(CS%Re_Ah_const_xy(IsdB:IedB,JsdB:JedB)); CS%Re_Ah_const_xy(:,:) = 0.0 + ALLOC_(CS%Re_Ah_const_xx(isd:ied,jsd:jed)) ; CS%Re_Ah_const_xx(:,:) = 0.0 + ALLOC_(CS%Re_Ah_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Re_Ah_const_xy(:,:) = 0.0 endif endif do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 CS%dx2q(I,J) = G%dxBu(I,J)*G%dxBu(I,J) ; CS%dy2q(I,J) = G%dyBu(I,J)*G%dyBu(I,J) - CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) enddo ; enddo + + if (((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) .and. & + ((G%isc-G%isd < 3) .or. (G%isc-G%isd < 3))) call MOM_error(FATAL, & + "The minimum halo size is 3 when a Leith viscosity is being used.") + if (CS%use_Leithy) then + do J=js-3,Jeq+2 ; do I=is-3,Ieq+2 + CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) + enddo ; enddo + elseif ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) + enddo ; enddo + else + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) + enddo ; enddo + endif + do j=js-2,Jeq+2 ; do i=is-2,Ieq+2 CS%dx2h(i,j) = G%dxT(i,j)*G%dxT(i,j) ; CS%dy2h(i,j) = G%dyT(i,j)*G%dyT(i,j) CS%DX_dyT(i,j) = G%dxT(i,j)*G%IdyT(i,j) ; CS%DY_dxT(i,j) = G%dyT(i,j)*G%IdxT(i,j) @@ -2541,12 +2554,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) endif endif if (CS%Leith_Ah) then - CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3) + 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 + 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 @@ -2571,12 +2584,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) endif endif if ((CS%Leith_Ah) .or. (CS%use_Leithy))then - CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3) + 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)) if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xy(i,j) = grid_sp_q3 / CS%Re_Ah if (Ah_time_scale > 0.) CS%Ah_bg_xy(i,j) = & - MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale) + MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then CS%Ah_Max_xy(I,J) = Ah_Limit * (grid_sp_q2 * grid_sp_q2) CS%Ah_bg_xy(I,J) = MIN(CS%Ah_bg_xy(I,J), CS%Ah_Max_xy(I,J)) @@ -2822,6 +2835,18 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) end subroutine hor_visc_init +!> hor_visc_vel_stencil returns the horizontal viscosity input velocity stencil size +function hor_visc_vel_stencil(CS) result(stencil) + type(hor_visc_CS), intent(in) :: CS !< Control structure for horizontal viscosity + integer :: stencil !< The horizontal viscosity velocity stencil size with the current settings. + + stencil = 2 + + if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then + stencil = 3 + endif +end function hor_visc_vel_stencil + !> Calculates factors in the anisotropic orientation tensor to be align with the grid. !! With n1=1 and n2=0, this recovers the approach of Large et al, 2001. subroutine align_aniso_tensor_to_grid(CS, n1, n2)