diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index b315916ec5..08afa29c00 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -841,7 +841,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f ! diffu = horizontal viscosity terms (u_av) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & + call horizontal_viscosity(u_av, v_av, h_av, uh, vh, CS%diffu, CS%diffv, & MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt, & OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & ADp=CS%ADp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v) @@ -1511,7 +1511,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. & .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, & + call horizontal_viscosity(u, v, h, uh, vh, CS%diffu, CS%diffv, MEKE, VarMix, G, GV, US, CS%hor_visc, & tv, dt, OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v) call set_initialized(CS%diffu, "diffu", restart_CS) diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index e55e2e3f96..2220d0e32a 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -552,7 +552,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc ! diffu = horizontal viscosity terms (u_av) 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, & + call horizontal_viscosity(u_av, v_av, h_av, uh, vh, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, & tv, dt, OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%AD_pred) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with predictor horizontal_viscosity (step_MOM_dyn_split_RK2b)") @@ -838,7 +838,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc ! diffu = horizontal viscosity terms (u_av) 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, tv, dt, & + call horizontal_viscosity(u_av, v_av, h_av, uh, vh, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt, & OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%ADp) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2b)") diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 579ddead2d..72c7dbe6cd 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -263,7 +263,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! diffu = horizontal viscosity terms (u,h) call enable_averages(dt, Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt) + call horizontal_viscosity(u, v, h, uh, vh, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 65b3bdf50e..cc37f1c2bc 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -275,7 +275,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! diffu = horizontal viscosity terms (u,h) call enable_averages(dt,Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u_in, v_in, h_in, CS%diffu, CS%diffv, MEKE, VarMix, & + call horizontal_viscosity(u_in, v_in, h_in, uh, vh, CS%diffu, CS%diffv, MEKE, VarMix, & G, GV, US, CS%hor_visc, tv, dt) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 2eef171bf5..9b9ea878c9 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -88,6 +88,8 @@ module MOM_hor_visc !! in setting the corner-point viscosities when USE_KH_BG_2D=True. real :: Kh_bg_min !< The minimum value allowed for Laplacian horizontal !! viscosity [L2 T-1 ~> m2 s-1]. The default is 0.0. + logical :: FrictWork_bug !< If true, retain an answer-changing bug + !! in calculating FrictWork, which cancels the h in thickness flux and the h at velocity point. logical :: use_land_mask !< Use the land mask for the computation of thicknesses !! at velocity locations. This eliminates the dependence on !! arbitrary values over land or outside of the domain. @@ -241,7 +243,7 @@ module MOM_hor_visc !! u(is-2:ie+2,js-2:je+2) !! v(is-2:ie+2,js-2:je+2) !! h(is-1:ie+1,js-1:je+1) or up to h(is-2:ie+2,js-2:je+2) with some Leith options. -subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & +subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, GV, US, & CS, tv, dt, OBC, BT, TD, ADp, hu_cont, hv_cont) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -251,6 +253,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: uh !< The zonal volume transport [H L2 T-1 ~> m3 s-1]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in) :: vh !< The meridional volume transport [H L2 T-1 ~> m3 s-1]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(out) :: diffu !< Zonal acceleration due to convergence of !! along-coordinate stress tensor [L T-2 ~> m s-2] @@ -1781,10 +1787,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo endif - if (find_FrictWork) then ; do j=js,je ; do i=is,ie + if (find_FrictWork) then + if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie ! 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 * ( & + 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) * & @@ -1799,12 +1806,44 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, + 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 + enddo ; enddo + else ; do j=js,je ; do i=is,ie + FrictWork(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & + (str_xx(i,j)*CS%dy2h(i,j) * ( & + uh(I,j,k)*G%dxCu(I,j)/(h_u(I,j)+h_neglect)*G%IdyCu(I,j)*G%IareaCu(I,j) & + -uh(I-1,j,k)*G%dxCu(I-1,j)/(h_u(I-1,j)+h_neglect)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)) & + -str_xx(i,j)*CS%dx2h(i,j) * ( & + vh(i,J,k)*G%dyCv(i,J)/(h_v(i,J)+h_neglect)*G%IdxCv(i,J)*G%IareaCv(i,J) & + -vh(i,J-1,k)*G%dyCv(i,J-1)/(h_v(i,J-1)+h_neglect)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1))) & + +0.25*((str_xy(I,J)*( & + CS%dx2q(I,J)*(uh(I,j+1,k)/(h_u(I,j+1)+h_neglect)*G%IareaCu(I,j+1) & + -uh(I,j,k)/(h_u(I,j)+h_neglect)*G%IareaCu(I,j)) & + +CS%dy2q(I,J)*(vh(i+1,J,k)/(h_v(i+1,J)+h_neglect)*G%IareaCv(i+1,J) & + -vh(i,J,k)/(h_v(i,J)+h_neglect)*G%IareaCv(i,J)) ) & + +str_xy(I-1,J-1)*( & + CS%dx2q(I-1,J-1)*(uh(I-1,j,k)/(h_u(I-1,j)+h_neglect)*G%IareaCu(I-1,j) & + -uh(I-1,j-1,k)/(h_u(I-1,j-1)+h_neglect)*G%IareaCu(I-1,j-1)) & + +CS%dy2q(I-1,J-1)*(vh(i,J-1,k)/(h_v(i,J-1)+h_neglect)*G%IareaCv(i,J-1) & + -vh(i-1,J-1,k)/(h_v(i-1,J-1)+h_neglect)*G%IareaCv(i-1,J-1)) )) & + +(str_xy(I-1,J)*( & + CS%dx2q(I-1,J)*(uh(I-1,j+1,k)/(h_u(I-1,j+1)+h_neglect)*G%IareaCu(I-1,j+1) & + -uh(I-1,j,k)/(h_u(I-1,j)+h_neglect)*G%IareaCu(I-1,j)) & + +CS%dy2q(I-1,J)*(vh(i,J,k)/(h_v(i,J)+h_neglect)*G%IareaCv(i,J) & + -vh(i-1,J,k)/(h_v(i-1,J)+h_neglect)*G%IareaCv(i-1,J)) ) & + +str_xy(I,J-1)*( & + CS%dx2q(I,J-1)*(uh(I,j,k)/(h_u(I,j)+h_neglect)*G%IareaCu(I,j) & + -uh(I,j-1,k)/(h_u(I,j-1)+h_neglect)*G%IareaCu(I,j-1)) & + +CS%dy2q(I,J-1)*(vh(i+1,J-1,k)/(h_v(i+1,J-1)+h_neglect)*G%IareaCv(i+1,J-1) & + -vh(i,J-1,k)/(h_v(i,J-1)+h_neglect)*G%IareaCv(i,J-1)) )) ) ) + enddo ; enddo ; endif + endif - if (CS%use_GME) then ; do j=js,je ; do i=is,ie + + if (CS%use_GME) then + if (CS%FrictWork_bug) 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 * ( & + 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) * & @@ -1819,7 +1858,37 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, + 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 + enddo ; enddo + else ; do j=js,je ; do i=is,ie + FrictWork_GME(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & + (str_xx_GME(i,j)*CS%dy2h(i,j) * ( & + uh(I,j,k)*G%dxCu(I,j)/(h_u(I,j)+h_neglect)*G%IdyCu(I,j)*G%IareaCu(I,j) & + -uh(I-1,j,k)*G%dxCu(I-1,j)/(h_u(I-1,j)+h_neglect)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)) & + -str_xx_GME(i,j)*CS%dx2h(i,j) * ( & + vh(i,J,k)*G%dyCv(i,J)/(h_v(i,J)+h_neglect)*G%IdxCv(i,J)*G%IareaCv(i,J) & + -vh(i,J-1,k)*G%dyCv(i,J-1)/(h_v(i,J-1)+h_neglect)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1))) & + +0.25*((str_xy_GME(I,J)*( & + CS%dx2q(I,J)*(uh(I,j+1,k)/(h_u(I,j+1)+h_neglect)*G%IareaCu(I,j+1) & + -uh(I,j,k)/(h_u(I,j)+h_neglect)*G%IareaCu(I,j)) & + +CS%dy2q(I,J)*(vh(i+1,J,k)/(h_v(i+1,J)+h_neglect)*G%IareaCv(i+1,J) & + -vh(i,J,k)/(h_v(i,J)+h_neglect)*G%IareaCv(i,J)) ) & + +str_xy_GME(I-1,J-1)*( & + CS%dx2q(I-1,J-1)*(uh(I-1,j,k)/(h_u(I-1,j)+h_neglect)*G%IareaCu(I-1,j) & + -uh(I-1,j-1,k)/(h_u(I-1,j-1)+h_neglect)*G%IareaCu(I-1,j-1)) & + +CS%dy2q(I-1,J-1)*(vh(i,J-1,k)/(h_v(i,J-1)+h_neglect)*G%IareaCv(i,J-1) & + -vh(i-1,J-1,k)/(h_v(i-1,J-1)+h_neglect)*G%IareaCv(i-1,J-1)) )) & + +(str_xy_GME(I-1,J)*( & + CS%dx2q(I-1,J)*(uh(I-1,j+1,k)/(h_u(I-1,j+1)+h_neglect)*G%IareaCu(I-1,j+1) & + -uh(I-1,j,k)/(h_u(I-1,j)+h_neglect)*G%IareaCu(I-1,j)) & + +CS%dy2q(I-1,J)*(vh(i,J,k)/(h_v(i,J)+h_neglect)*G%IareaCv(i,J) & + -vh(i-1,J,k)/(h_v(i-1,J)+h_neglect)*G%IareaCv(i-1,J)) ) & + +str_xy_GME(I,J-1)*( & + CS%dx2q(I,J-1)*(uh(I,j,k)/(h_u(I,j)+h_neglect)*G%IareaCu(I,j) & + -uh(I,j-1,k)/(h_u(I,j-1)+h_neglect)*G%IareaCu(I,j-1)) & + +CS%dy2q(I,J-1)*(vh(i+1,J-1,k)/(h_v(i+1,J-1)+h_neglect)*G%IareaCv(i+1,J-1) & + -vh(i,J-1,k)/(h_v(i,J-1)+h_neglect)*G%IareaCv(i,J-1)) )) ) ) + enddo ; enddo ; endif + endif ! Make a similar calculation as for FrictWork above but accumulating into ! the vertically integrated MEKE source term, and adjusting for any @@ -2282,6 +2351,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "If true, retain an answer-changing horizontal indexing bug in setting "//& "the corner-point viscosities when USE_KH_BG_2D=True. This is "//& "not recommended.", default=.false., do_not_log=.not.CS%use_Kh_bg_2d) + call get_param(param_file, mdl, "FRICTWORK_BUG", CS%FrictWork_bug, & + "If true, retain an answer-changing bug in calculating "//& + "the FrictWork, which cancels the h in thickness flux and the h at velocity point. This is"//& + "not recommended.", default=.true.) call get_param(param_file, mdl, "USE_GME", CS%use_GME, & "If true, use the GM+E backscatter scheme in association \n"//&