Correct the computation of FrictWork in MOM_hor_visc

Add FRICTWORK_BUG parameter (default is true). If FRICTWORK_BUG is false, use the new way to compute FrictWork using the thickness flux. The previous version (realized by setting FRICTWORK_BUG as true) uses the velocity to compute FrictWork, which overestimates the total energy dissipation rate compared with the domain integral of KE_horvisc computed in MOM_diagnostics.

src/core/MOM_dynamics_split_RK2.F90

@@ -846,7 +846,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)
@@ -1530,7 +1530,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)

src/core/MOM_dynamics_split_RK2b.F90

@@ -553,7 +553,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)")
@@ -843,8 +843,8 @@ 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, &
- OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%ADp)
+ 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)")
 

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)
 

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)

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -94,6 +94,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.
@@ -247,7 +249,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.
@@ -257,6 +259,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]
@@ -615,7 +621,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
  !$OMP parallel do default(none) &
  !$OMP shared( &
- !$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, &
+ !$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, uh, vh, &
  !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &
  !$OMP is_vort, ie_vort, js_vort, je_vort, &
  !$OMP is_Kh, ie_Kh, js_Kh, je_Kh, &
@@ -1787,10 +1793,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) * &
@@ -1805,12 +1812,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)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
+ - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) &
+ - (str_xx(i,j)*CS%dx2h(i,j) * ( &
+ (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
+ - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) &
+ + (0.25*(((str_xy(I,J)*( &
+ (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) &
+ - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) &
+ + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) &
+ - (vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)))) )) &
+ +(str_xy(I-1,J-1)*( &
+ (CS%dx2q(I-1,J-1)*((uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) &
+ - (uh(I-1,j-1,k)*G%IareaCu(I-1,j-1)/(h_u(I-1,j-1)+h_neglect)))) &
+ + (CS%dy2q(I-1,J-1)*((vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) &
+ - (vh(i-1,J-1,k)*G%IareaCv(i-1,J-1)/(h_v(i-1,J-1)+h_neglect)))) )) ) &
+ +((str_xy(I-1,J)*( &
+ (CS%dx2q(I-1,J)*((uh(I-1,j+1,k)*G%IareaCu(I-1,j+1)/(h_u(I-1,j+1)+h_neglect)) &
+ - (uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)))) &
+ + (CS%dy2q(I-1,J)*((vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
+ - (vh(i-1,J,k)*G%IareaCv(i-1,J)/(h_v(i-1,J)+h_neglect)))) )) &
+ +(str_xy(I,J-1)*( &
+ (CS%dx2q(I,J-1)*((uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
+ - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) &
+ + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) &
+ - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) )
+ 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) * &
@@ -1825,7 +1864,38 @@ 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)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
+ - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) &
+ - (str_xx_GME(i,j)*CS%dx2h(i,j) * ( &
+ (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
+ - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) &
+ + (0.25*(((str_xy_GME(I,J)*( &
+ (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) &
+ - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) &
+ + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) &
+ - (vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)))) )) &
+ +(str_xy_GME(I-1,J-1)*( &
+ (CS%dx2q(I-1,J-1)*((uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) &
+ - (uh(I-1,j-1,k)*G%IareaCu(I-1,j-1)/(h_u(I-1,j-1)+h_neglect)))) &
+ + (CS%dy2q(I-1,J-1)*((vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) &
+ - (vh(i-1,J-1,k)*G%IareaCv(i-1,J-1)/(h_v(i-1,J-1)+h_neglect)))) )) ) &
+ +((str_xy_GME(I-1,J)*( &
+ (CS%dx2q(I-1,J)*((uh(I-1,j+1,k)*G%IareaCu(I-1,j+1)/(h_u(I-1,j+1)+h_neglect)) &
+ - (uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)))) &
+ + (CS%dy2q(I-1,J)*((vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
+ - (vh(i-1,J,k)*G%IareaCv(i-1,J)/(h_v(i-1,J)+h_neglect)))) )) &
+ +(str_xy_GME(I,J-1)*( &
+ (CS%dx2q(I,J-1)*((uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
+ - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) &
+ + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) &
+ - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) )
+
+ 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
@@ -2298,6 +2368,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"//&