Correct the computation of FrictWork

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

@@ -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)

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)")

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

@@ -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"//&