diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 8c016b11b0..ca94af2225 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -681,7 +681,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, & - OBC=CS%OBC, BT=CS%barotropic_CSp) + OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)") @@ -1149,7 +1149,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param .not. query_initialized(CS%diffv,"diffv",restart_CS)) then call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, & G, GV, US, CS%hor_visc_CSp, & - OBC=CS%OBC, BT=CS%barotropic_CSp) + OBC=CS%OBC, BT=CS%barotropic_CSp, & + TD=thickness_diffuse_CSp) else if ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & (US%m_to_L * US%s_to_T_restart**2 /= US%m_to_L_restart * US%s_to_T**2) ) then diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index ffd3b1ac63..d9b465ee0a 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -203,7 +203,7 @@ module MOM_hor_visc !! v[is-2:ie+2,js-2:je+2] !! h[is-1:ie+1,js-1:je+1] subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & - CS, OBC, BT) + CS, OBC, BT, TD) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & @@ -228,7 +228,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type type(barotropic_CS), optional, pointer :: BT !< Pointer to a structure containing !! barotropic velocities. - + type(thickness_diffuse_CS), optional, pointer :: TD !< Pointer to a structure containing + !! thickness diffusivities. ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & Del2u, & ! The u-compontent of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1] @@ -411,15 +412,13 @@ 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) - do j=js-1,je+1 ; do i=is-1,ie+1 + do j=js-1,je+2 ; do i=is-1,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 - call pass_vector(dudx_bt, dvdy_bt, G%Domain, stagger=BGRID_NE) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j) enddo ; enddo @@ -432,6 +431,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, - ubtav(I,j)*G%IdxCu(I,j)) enddo ; enddo + call pass_vector(dudx_bt, dvdy_bt, G%Domain, stagger=BGRID_NE) call pass_vector(dvdx_bt, dudy_bt, G%Domain, stagger=AGRID) if (CS%no_slip) then @@ -1063,48 +1063,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo if (CS%use_GME) then - if (CS%answers_2018) then - do j=js,je ; do i=is,ie - grad_vel_mag_h(i,j) = boundary_mask_h(i,j) * (dudx(i,j)**2 + dvdy(i,j)**2 + & - (0.25*((dvdx(I,J) + dvdx(I-1,J-1)) + (dvdx(I,J-1) + dvdx(I-1,J))) )**2 + & - (0.25*((dudy(I,J) + dudy(I-1,J-1)) + (dudy(I,J-1) + dudy(I-1,J))) )**2) - max_diss_rate_h(i,j,k) = 2.0 * MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j)) - enddo ; enddo - else ! This form is invariant to 90-degree rotations. - do j=js,je ; do i=is,ie - grad_vel_mag_h(i,j) = boundary_mask_h(i,j) * ((dudx(i,j)**2 + dvdy(i,j)**2) + & - ((0.25*((dvdx(I,J) + dvdx(I-1,J-1)) + (dvdx(I,J-1) + dvdx(I-1,J))) )**2 + & - (0.25*((dudy(I,J) + dudy(I-1,J-1)) + (dudy(I,J-1) + dudy(I-1,J))) )**2)) - max_diss_rate_h(i,j,k) = 2.0 * MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j)) - enddo ; enddo - endif - - if (CS%answers_2018) then - do J = G%JscB, G%JecB ; do I = G%IscB, G%IecB - grad_vel_mag_q(I,J) = boundary_mask_q(I,J) * (dudx(i,j)**2 + dvdy(i,j)**2 + & - (0.25*((dvdx(I,J)+dvdx(I-1,J-1)) + (dvdx(I,J-1)+dvdx(I-1,J))) )**2 + & - (0.25*((dudy(I,J)+dudy(I-1,J-1)) + (dudy(I,J-1)+dudy(I-1,J))) )**2) - - max_diss_rate_q(I,J,k) = 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i+1,j)+ & - MEKE%MEKE(i,j+1)+MEKE%MEKE(i+1,j+1)) * sqrt(grad_vel_mag_q(I,J)) - enddo ; enddo - else ! This form is rotationally invariant - do J = G%JscB, G%JecB ; do I = G%IscB, G%IecB - grad_vel_mag_q(I,J) = boundary_mask_q(I,J) * ((dudx(i,j)**2 + dvdy(i,j)**2) + & - ((0.25*((dvdx(I,J)+dvdx(I-1,J-1)) + (dvdx(I,J-1)+dvdx(I-1,J))) )**2 + & - (0.25*((dudy(I,J)+dudy(I-1,J-1)) + (dudy(I,J-1)+dudy(I-1,J))) )**2)) - - max_diss_rate_q(I,J,k) = 0.5*((MEKE%MEKE(i,j) + MEKE%MEKE(i+1,j+1)) + & - (MEKE%MEKE(i+1,j) + MEKE%MEKE(i,j+1))) * sqrt(grad_vel_mag_q(I,J)) - enddo ; enddo - endif + call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G) + call pass_vector(KH_u_GME, KH_v_GME, G%Domain) - do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if ((grad_vel_mag_bt_h(i,j)>0) .and. (max_diss_rate_h(i,j,k)>0)) then - GME_coeff = (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * CS%GME_efficiency*max_diss_rate_h(i,j,k) / & - grad_vel_mag_bt_h(i,j) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (grad_vel_mag_bt_h(i,j)>0) then + GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * & + (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))) else - GME_coeff = 1.0 + GME_coeff = 0.0 endif ! apply mask @@ -1117,10 +1084,9 @@ 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 - - if ((grad_vel_mag_bt_q(I,J)>0) .and. (max_diss_rate_q(I,J,k)>0)) then - GME_coeff = (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * CS%GME_efficiency*max_diss_rate_q(I,J,k) / & - grad_vel_mag_bt_q(I,J) + if (grad_vel_mag_bt_q(i,j)>0) then + GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * & + (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))) else GME_coeff = 0.0 endif @@ -1153,7 +1119,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (associated(MEKE%GME_snk)) then do j=js,je ; do i=is,ie - FrictWork_GME(i,j,k) = GME_coeff_h(i,j,k) * h(i,j,k) * GV%H_to_RZ * grad_vel_mag_bt_h(i,j) + FrictWork_GME(i,j,k) = GME_coeff_h(i,j,k) * h(i,j,k) * GV%H_to_kg_m2 * grad_vel_mag_bt_h(i,j) enddo ; enddo endif @@ -1395,8 +1361,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) ! valid parameters. logical :: split ! If true, use the split time stepping scheme. ! If false and USE_GME = True, issue a FATAL error. - logical :: use_MEKE ! If true, use the MEKE module for calculating eddy kinetic energy. - ! If false and USE_GME = True, issue a FATAL error. logical :: default_2018_answers character(len=64) :: inputdir, filename @@ -1667,14 +1631,6 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) if (.not. split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & "cannot be used with SPLIT=False.") - call get_param(param_file, mdl, "USE_MEKE", use_MEKE, & - "If true, turns on the MEKE scheme which calculates\n"// & - "a sub-grid mesoscale eddy kinetic energy budget.", & - default=.false.) - - if (.not. use_MEKE) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & - "cannot be used with USE_MEKE=False.") - call get_param(param_file, mdl, "GME_H0", CS%GME_h0, & "The strength of GME tapers quadratically to zero when the bathymetric "//& "depth is shallower than GME_H0.", units="m", scale=US%m_to_Z, &