Added a GME sink term to the MEKE budget. Cleaned up MOM_hor_visc.F90.

src/parameterizations/lateral/MOM_MEKE.F90

@@ -30,6 +30,7 @@ module MOM_MEKE
   ! Parameters
   real :: MEKE_FrCoeff  !< Efficiency of conversion of ME into MEKE (non-dim)
   real :: MEKE_GMcoeff  !< Efficiency of conversion of PE into MEKE (non-dim)
+  real :: MEKE_GMECoeff !< Efficiency of conversion of MEKE into ME by GME (non-dim)
   real :: MEKE_damping  !< Local depth-independent MEKE dissipation rate in s-1.
   real :: MEKE_Cd_scale !< The ratio of the bottom eddy velocity to the column mean
                         !! eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1
@@ -73,7 +74,7 @@ module MOM_MEKE
   !>@{ Diagnostic handles
   integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1
   integer :: id_Ub = -1, id_Ut = -1
-  integer :: id_GM_src = -1, id_mom_src = -1, id_decay = -1
+  integer :: id_GM_src = -1, id_mom_src = -1, id_GME_snk = -1, id_decay = -1
   integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1, id_Au = -1
   integer :: id_Le = -1, id_gamma_b = -1, id_gamma_t = -1
   integer :: id_Lrhines = -1, id_Leady = -1
@@ -113,6 +114,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
     MEKE_decay, &   ! The MEKE decay timescale, in s-1.
     MEKE_GM_src, &  ! The MEKE source from thickness mixing, in m2 s-3.
     MEKE_mom_src, & ! The MEKE source from momentum, in m2 s-3.
+    MEKE_GME_snk, & ! The MEKE sink from GME backscatter, in m2 s-3.
     drag_rate_visc, &
     drag_rate, &    ! The MEKE spindown timescale due to bottom drag, in s-1.
     LmixScale, &    ! Square of eddy mixing length, in m2.
@@ -165,6 +167,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
 
     if (CS%debug) then
       if (associated(MEKE%mom_src)) call hchksum(MEKE%mom_src, 'MEKE mom_src',G%HI)
+      if (associated(MEKE%GME_snk)) call hchksum(MEKE%GME_snk, 'MEKE GME_snk',G%HI)
       if (associated(MEKE%GM_src)) call hchksum(MEKE%GM_src, 'MEKE GM_src',G%HI)
       if (associated(MEKE%MEKE)) call hchksum(MEKE%MEKE, 'MEKE MEKE',G%HI)
       call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI)
@@ -292,6 +295,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
       enddo ; enddo
     endif
 
+    if (associated(MEKE%GME_snk)) then
+!$OMP do
+      do j=js,je ; do i=is,ie
+        src(i,j) = src(i,j) - CS%MEKE_GMECoeff*I_mass(i,j)*MEKE%GME_snk(i,j)
+      enddo ; enddo
+    endif
+
     if (associated(MEKE%GM_src)) then
 !$OMP do
       if (CS%GM_src_alt) then
@@ -578,6 +588,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, CS, hu, hv)
     if (CS%id_decay>0) call post_data(CS%id_decay, MEKE_decay, CS%diag)
     if (CS%id_GM_src>0) call post_data(CS%id_GM_src, MEKE%GM_src, CS%diag)
     if (CS%id_mom_src>0) call post_data(CS%id_mom_src, MEKE%mom_src, CS%diag)
+    if (CS%id_GME_snk>0) call post_data(CS%id_GME_snk, MEKE%GME_snk, CS%diag)
     if (CS%id_Le>0) call post_data(CS%id_Le, LmixScale, CS%diag)
     if (CS%id_gamma_b>0) then
       do j=js,je ; do i=is,ie
@@ -892,6 +903,10 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
                  "The efficiency of the conversion of mean energy into \n"//&
                  "MEKE.  If MEKE_FRCOEFF is negative, this conversion \n"//&
                  "is not used or calculated.", units="nondim", default=-1.0)
+  call get_param(param_file, mdl, "MEKE_GMECOEFF", CS%MEKE_GMECoeff, &
+                 "The efficiency of the conversion of MEKE into mean energy \n"//&
+                 "by GME.  If MEKE_GMECOEFF is negative, this conversion \n"//&
+                 "is not used or calculated.", units="nondim", default=-1.0)
   call get_param(param_file, mdl, "MEKE_BGSRC", CS%MEKE_BGsrc, &
                  "A background energy source for MEKE.", units="W kg-1", &
                  default=0.0)
@@ -1065,6 +1080,9 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS)
   CS%id_mom_src = register_diag_field('ocean_model', 'MEKE_mom_src',diag%axesT1, Time, &
      'MEKE energy available from momentum', 'W m-2')
   if (.not. associated(MEKE%mom_src)) CS%id_mom_src = -1
+  CS%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, Time, &
+     'MEKE energy lost to GME backscatter', 'W m-2')
+  if (.not. associated(MEKE%GME_snk)) CS%id_GME_snk = -1
   CS%id_Le = register_diag_field('ocean_model', 'MEKE_Le', diag%axesT1, Time, &
      'Eddy mixing length used in the MEKE derived eddy diffusivity', 'm')
   CS%id_Lrhines = register_diag_field('ocean_model', 'MEKE_Lrhines', diag%axesT1, Time, &
@@ -1097,7 +1115,7 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS)
   type(MOM_restart_CS),  pointer       :: restart_CS !< Restart control structure for MOM_MEKE.
 ! Local variables
   type(vardesc) :: vd
-  real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_KHCoeff, MEKE_viscCoeff
+  real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_GMECoeff, MEKE_KHCoeff, MEKE_viscCoeff
   logical :: useMEKE
   integer :: isd, ied, jsd, jed
 
@@ -1107,6 +1125,7 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS)
 ! Read these parameters to determine what should be in the restarts
   MEKE_GMcoeff =-1.; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff)
   MEKE_FrCoeff =-1.; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff)
+  MEKE_GMEcoeff =-1.; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff)
   MEKE_KhCoeff =1.; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff)
   MEKE_viscCoeff =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF",MEKE_viscCoeff)
 
@@ -1132,6 +1151,9 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS)
   if (MEKE_FrCoeff>=0.) then
     allocate(MEKE%mom_src(isd:ied,jsd:jed)) ; MEKE%mom_src(:,:) = 0.0
   endif
+  if (MEKE_GMECoeff>=0.) then
+    allocate(MEKE%GME_snk(isd:ied,jsd:jed)) ; MEKE%GME_snk(:,:) = 0.0
+  endif
   if (MEKE_KhCoeff>=0.) then
     allocate(MEKE%Kh(isd:ied,jsd:jed)) ; MEKE%Kh(:,:) = 0.0
     vd = var_desc("MEKE_Kh", "m2 s-1",hor_grid='h',z_grid='1', &
@@ -1166,6 +1188,7 @@ subroutine MEKE_end(MEKE, CS)
   if (associated(MEKE%MEKE)) deallocate(MEKE%MEKE)
   if (associated(MEKE%GM_src)) deallocate(MEKE%GM_src)
   if (associated(MEKE%mom_src)) deallocate(MEKE%mom_src)
+  if (associated(MEKE%GME_snk)) deallocate(MEKE%GME_snk)
   if (associated(MEKE%Kh)) deallocate(MEKE%Kh)
   if (associated(MEKE%Ku)) deallocate(MEKE%Ku)
   if (associated(MEKE%Au)) deallocate(MEKE%Au)

src/parameterizations/lateral/MOM_MEKE_types.F90

@@ -11,6 +11,7 @@ module MOM_MEKE_types
     MEKE => NULL(), &   !< Vertically averaged eddy kinetic energy, in m2 s-2.
     GM_src => NULL(), & !< MEKE source due to thickness mixing (GM), in W m-2.
     mom_src => NULL(),& !< MEKE source from lateral friction in the momentum equations, in W m-2.
+    GME_snk => NULL(),& !< MEKE sink from GME backscatter in the momentum equations, in W m-2.
     Kh => NULL(), &     !< The MEKE-derived lateral mixing coefficient in m2 s-1.
     Rd_dx_h => NULL()   !< The deformation radius compared with the grid spacing, nondim.
                         !! Rd_dx_h is copied from VarMix_CS.

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -76,7 +76,7 @@ module MOM_hor_visc
   logical :: dynamic_aniso   !< If true, the anisotropic viscosity is recomputed as a function
                              !! of state. This is set depending on ANISOTROPIC_MODE.
   logical :: use_GME         !< If true, use GME backscatter scheme.
-  real    :: GME_dt
+
   real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Kh_bg_xx
                       !< The background Laplacian viscosity at h points, in units
                       !! of m2 s-1. The actual viscosity may be the larger of this
@@ -326,11 +326,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
   real :: local_strain ! Local variable for interpolating computed strain rates (s-1).
   real :: epsilon
   real :: GME_coeff ! The GME (negative) viscosity coefficient (m2 s-1)
-  real :: GME_coeff_limiter
-  real :: FWfrac
+  real :: GME_coeff_limiter ! Maximum permitted value of the GME coefficient (m2 s-1)
+  real :: FWfrac ! Fraction of maximum theoretical energy transfer to use when scaling GME coefficient 
   real :: DY_dxBu, DX_dyBu
-  real :: H0
-  real :: dr, dv
+  real :: H0 ! Depth used to scale down GME coefficient in shallow areas (m)
   logical :: rescale_Kh, legacy_bound
   logical :: find_FrictWork
   logical :: apply_OBC = .false.
@@ -396,6 +395,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
     ! GME tapers off above this depth
     H0 = 1000.0
     FWfrac = 0.1
+    GME_coeff_limiter = 1e7
+
     ! initialize diag. array with zeros
     GME_coeff_h(:,:,:) = 0.0
     GME_coeff_q(:,:,:) = 0.0
@@ -1091,7 +1092,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
         if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then
           ! This is the maximum possible amount of energy that can be converted
           ! per unit time, according to theory (multiplied by h)
-          max_diss_rate(i,j,k) = 2.0*MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j))
+          max_diss_rate(i,j,k) = 2.0*MEKE%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j)) * (MIN(G%bathyT(i,j)/H0,1.0)**2)
 
           FrictWorkMax(i,j,k) = max_diss_rate(i,j,k) * h(i,j,k) * GV%H_to_kg_m2
 
@@ -1111,64 +1112,24 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
 
 
     if (CS%use_GME) then
+
+      if (.not. (associated(MEKE))) call MOM_error(FATAL, &
+        "MEKE must be enabled for GME to be used.")
 
-      do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      if (.not. (associated(MEKE%mom_src))) call MOM_error(FATAL, &
+        "MEKE%mom_src must be enabled for GME to be used.")
 
-!        GME_coeff = (MIN(G%bathyT(i,j)/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)) * &
-!                    (0.5*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)),0.0) * &
-!                    ( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I-1,j,k)) )**2 + &
-!                      (0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i,J-1,k)) )**2 ) / &
-!                    ( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
-!                      (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
-!                      (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + &
-!                     epsilon))
-
-!        GME_coeff = 2.0 * (MIN(G%bathyT(i,j)/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)) * &
-!                    sqrt(0.5*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)),0.0)) / &
-!                   sqrt( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
-!                      (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
-!                      (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + &
-!                     epsilon)
-
-      if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then
-!        GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * target_FrictWork_GME(i,j,k) / ( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
-!                      (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
-!                      (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + &
-!                      epsilon)          
-!        GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * max_diss_rate(i,j,k) * (G%areaT(i,j)/0.01) 
-
-!           GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * MEKE%MEKE(i,j) / MAX(0.25*(VarMix%SN_u(I,j)+VarMIX%SN_u(I-1,j)+VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)),epsilon)
-        GME_coeff = 1e10; dr = 1e6
-
-!        if (G%mask2dT(i,j)) 
-!          do while (dr >= max_diss_rate(i,j,k))
-!            GME_coeff = GME_coeff / 2.0
-            dr = GME_coeff * grad_vel_mag_h(i,j)
+      do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
 
-            if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_h(i,j)>0) ) then
-              GME_coeff = FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_h(i,j)      
-!              GME_coeff = max_diss_rate_bt(i,j) / grad_vel_mag_bt_h(i,j)
-            else
-              GME_coeff = 0.0
-            endif
-
-!            GME_coeff = (MIN(G%bathyT(i,j)*G%Zd_to_m/H0,1.0)**2)*GME_coeff
-
-!           GME_coeff = 10.0
-!          enddo     
-!        else
-!          GME_coeff = 0.0
-!        endif
-
-
-      endif ; endif
+        if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_h(i,j)>0) ) then
+          GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_h(i,j)      
+        else
+          GME_coeff = 0.0
+        endif
 
         ! apply mask
         GME_coeff = GME_coeff * (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))
 
-        GME_coeff_limiter = 1e7
         GME_coeff = MIN(GME_coeff,GME_coeff_limiter)
 
         if ((CS%id_GME_coeff_h>0) .or. find_FrictWork) GME_coeff_h(i,j,k) = GME_coeff
@@ -1180,74 +1141,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
 
       do J=js-1,Jeq ; do I=is-1,Ieq
 
-!        GME_coeff = (MIN(G%bathyT(i,j)/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)) * &
-!                     ( 0.25*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I,j+1,k) + &
-!                            VarMix%N2_v(i,J,k)+VarMix%N2_v(i+1,J,k)),0.0) * &
-!                    ( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I,j+1,k)) )**2 + &
-!                      (0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i+1,J,k)) )**2 ) / &
-!                      ( dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
-!                        (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
-!                       (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + &
-!                         epsilon))
-
-!        GME_coeff = 2.0 * (MIN(G%bathyT(i,j)/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)) * &
-!                     sqrt( 0.25*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I,j+1,k) + &
-!                            VarMix%N2_v(i,J,k)+VarMix%N2_v(i+1,J,k)),0.0)) / &
-!                      sqrt(dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
-!                        (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
-!                       (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + &
-!                         epsilon)
-
-      if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then
-!        GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * target_FrictWork_GME(i,j,k) / (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
-!                        (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
-!                        (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + &
-!                        epsilon)
-!        GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * target_diss_rate_GME(i,j,k) * G%areaT(i,j)  / &
-!                       (0.1**2)
-!         GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * target_diss_rate_GME(i,j,k) * sqrt(G%areaT(i,j))/0.0001
-!
-!           GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * MEKE%MEKE(i,j) / MAX(0.25*(VarMix%SN_u(I,j)+VarMIX%SN_u(I-1,j)+VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)),epsilon)        
-!
-
-!      if (G%mask2dT(i,j))
-        GME_coeff = 1e10; dr = 1e6
-
-            if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_q(i,j)>0) ) then
-              GME_coeff = FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_q(I,J)
-!              GME_coeff = max_diss_rate_bt(i,j) / grad_vel_mag_bt_q(I,J)
-            else
-              GME_coeff = 0.0
-            endif
-
-!            dr = GME_coeff * bt_grad_mag_q(I,J)
-!            if ((max_diss_rate(i,j,k) > 0) .and. (dr>0) ) then
-!              dv = dr / (FWfrac*max_diss_rate(i,j,k))
-!              GME_coeff = GME_coeff / dv
-!            else
-!              GME_coeff = 0.0
-!            endif
+        if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_q(i,j)>0) ) then
+          GME_coeff = (MIN(G%bathyT(i,j)/H0,1.0)**2) * FWfrac*max_diss_rate(i,j,k) / grad_vel_mag_bt_q(I,J)
+        else 
+          GME_coeff = 0.0
+        endif
 
-!            GME_coeff = (MIN(G%bathyT(i,j)*G%Zd_to_m/H0,1.0)**2)*GME_coeff
-!        do while (dr >=  max_diss_rate(i,j,k)) 
-!             GME_coeff = GME_coeff / 2.0
-!             dr = GME_coeff * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
-!             (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
-!             (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2)
-!        enddo
-!      else
-!        GME_coeff = 0.0
-!      endif
-!            GME_coeff = 10.0
-
-      endif ; endif
-
         ! apply mask
         GME_coeff = GME_coeff * (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))
 
-        GME_coeff_limiter = 1e7
         GME_coeff = MIN(GME_coeff,GME_coeff_limiter)
 
         if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff
@@ -1272,11 +1174,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
         endif
       enddo ; enddo
 
-    if (find_FrictWork) 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_kg_m2 * grad_vel_mag_bt_h(i,j)
-      enddo ; enddo
-    endif
+      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_kg_m2 * grad_vel_mag_bt_h(i,j)
+        enddo ; enddo
+      endif
 
     else ! use_GME
       do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
@@ -1401,11 +1303,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
          ! MEKE%mom_src now is sign definite because it only uses the dissipation 
          MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + MAX(-FrictWorkMax(i,j,k),FrictWork_diss(i,j,k))
         enddo ; enddo
+
         if (CS%use_GME) then
-          do j=js,je ; do i=is,ie
-            ! MEKE%mom_src now is sign definite because it only uses the dissipation 
-            MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork_GME(i,j,k)
-          enddo ; enddo
+          if (associated(MEKE%GME_snk)) then
+            do j=js,je ; do i=is,ie
+              ! MEKE%mom_src now is sign definite because it only uses the dissipation 
+              MEKE%GME_snk(i,j) = MEKE%GME_snk(i,j) + FrictWork_GME(i,j,k)
+            enddo ; enddo
+          endif
         endif
       endif
     endif ; endif
@@ -1864,8 +1769,6 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS)
       CS%reduction_xy(I,J) = G%dx_Cv(i+1,J) / G%dxCv(i+1,J)
   enddo ; enddo
 
-  CS%GME_dt = dt
-
   if (CS%Laplacian) then
    ! The 0.3 below was 0.4 in MOM1.10.  The change in hq requires
    ! this to be less than 1/3, rather than 1/2 as before.