mom-ocean · marshallward · Aug 31, 2021 · Dec 23, 2020 · Jan 3, 2021 · Jan 3, 2021
diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
@@ -2685,6 +2685,16 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
       ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k)
     enddo ; enddo ; enddo
   endif
+  if (present(ADp) .and. (associated(ADp%visc_rem_u))) then
+    do k=1,nz ; do j=js,je ; do I=is-1,ie
+      ADp%visc_rem_u(I,j,k) = visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+  endif
+  if (present(ADp) .and. (associated(ADp%visc_rem_u))) then
+    do k=1,nz ; do J=js-1,je ; do i=is,ie
+      ADp%visc_rem_v(i,J,k) = visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+  endif
 
   if (G%nonblocking_updates) then
     if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain)

diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
@@ -169,10 +169,12 @@ module MOM_dynamics_split_RK2
   integer :: id_h_PFu    = -1, id_h_PFv    = -1
   integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1
   integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1
+  integer :: id_PFu_visc_rem = -1, id_PFv_visc_rem = -1
   ! integer :: id_hf_CAu    = -1, id_hf_CAv    = -1
   integer :: id_h_CAu    = -1, id_h_CAv    = -1
   integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1
   integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1
+  integer :: id_CAu_visc_rem = -1, id_CAv_visc_rem = -1
 
   ! Split scheme only.
   integer :: id_uav        = -1, id_vav        = -1
@@ -181,6 +183,7 @@ module MOM_dynamics_split_RK2
   integer :: id_h_u_BT_accel    = -1, id_h_v_BT_accel    = -1
   integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1
   integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1
+  integer :: id_u_BT_accel_visc_rem    = -1, id_v_BT_accel_visc_rem    = -1
   !>@}
 
   type(diag_ctrl), pointer       :: diag !< A structure that is used to regulate the
@@ -360,6 +363,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   real, dimension(SZI_(G),SZJB_(G)) :: &
     intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [H L T-2 ~> m2 s-2].
 
+  ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
+  real, allocatable, dimension(:,:,:) :: &
+    PFu_visc_rem, PFv_visc_rem, & ! Pressure force accel. x visc_rem_[uv] [L T-2 ~> m s-2].
+    CAu_visc_rem, CAv_visc_rem, & ! Coriolis force accel. x visc_rem_[uv] [L T-2 ~> m s-2].
+    u_BT_accel_visc_rem, v_BT_accel_visc_rem ! barotropic correction accel. x visc_rem_[uv] [L T-2 ~> m s-2].
+
   real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
 
   logical :: dyn_p_surf
@@ -1108,6 +1117,61 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     deallocate(h_v_BT_accel)
   endif
 
+  if (CS%id_PFu_visc_rem > 0) then
+    allocate(PFu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
+    PFu_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      PFu_visc_rem(I,j,k) = CS%PFu(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_PFu_visc_rem, PFu_visc_rem, CS%diag)
+    deallocate(PFu_visc_rem)
+  endif
+  if (CS%id_PFv_visc_rem > 0) then
+    allocate(PFv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
+    PFv_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      PFv_visc_rem(i,J,k) = CS%PFv(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_PFv_visc_rem, PFv_visc_rem, CS%diag)
+    deallocate(PFv_visc_rem)
+  endif
+  if (CS%id_CAu_visc_rem > 0) then
+    allocate(CAu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
+    CAu_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      CAu_visc_rem(I,j,k) = CS%CAu(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_CAu_visc_rem, CAu_visc_rem, CS%diag)
+    deallocate(CAu_visc_rem)
+  endif
+  if (CS%id_CAv_visc_rem > 0) then
+    allocate(CAv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
+    CAv_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      CAv_visc_rem(i,J,k) = CS%CAv(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_CAv_visc_rem, CAv_visc_rem, CS%diag)
+    deallocate(CAv_visc_rem)
+  endif
+  if (CS%id_u_BT_accel_visc_rem > 0) then
+    allocate(u_BT_accel_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
+    u_BT_accel_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      u_BT_accel_visc_rem(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_u_BT_accel_visc_rem, u_BT_accel_visc_rem, CS%diag)
+    deallocate(u_BT_accel_visc_rem)
+  endif
+  if (CS%id_v_BT_accel_visc_rem > 0) then
+    allocate(v_BT_accel_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
+    v_BT_accel_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      v_BT_accel_visc_rem(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_v_BT_accel_visc_rem, v_BT_accel_visc_rem, CS%diag)
+    deallocate(v_BT_accel_visc_rem)
+  endif
+
   if (CS%debug) then
     call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
     call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
@@ -1612,6 +1676,33 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
+  CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, &
+      'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if (CS%id_PFu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
+  CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, &
+      'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+
+  CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, &
+      'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if (CS%id_CAu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
+  CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, &
+      'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+
+  CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, &
+      'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if (CS%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
+  CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, &
+      'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
   id_clock_continuity = cpu_clock_id('(Ocean continuity equation)',      grain=CLOCK_MODULE)
   id_clock_pres       = cpu_clock_id('(Ocean pressure force)',           grain=CLOCK_MODULE)

diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
@@ -195,6 +195,9 @@ module MOM_variables
   real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points
   real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points
 
+  real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
+  real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
+
 end type accel_diag_ptrs
 
 !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances.

diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
@@ -2341,7 +2341,6 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
     call safe_alloc_ptr(ADp%gradKEu,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%gradKEv,isd,ied,JsdB,JedB,nz)
   endif
-
   if (associated(CS%KE_visc)) then
     call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
@@ -2395,7 +2394,7 @@ subroutine MOM_diagnostics_end(CS, ADp, CDp)
   if (associated(CS%vhGM_Rlay))  deallocate(CS%vhGM_Rlay)
 
   if (associated(ADp%gradKEu))    deallocate(ADp%gradKEu)
-  if (associated(ADp%gradKEu))    deallocate(ADp%gradKEu)
+  if (associated(ADp%gradKEv))    deallocate(ADp%gradKEv)
   if (associated(ADp%du_dt_visc)) deallocate(ADp%du_dt_visc)
   if (associated(ADp%dv_dt_visc)) deallocate(ADp%dv_dt_visc)
   if (associated(ADp%du_dt_str))  deallocate(ADp%du_dt_str)

diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -133,8 +133,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     MEKE_decay, &   ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1].
     drag_rate_visc, & ! Near-bottom velocity contribution to bottom dratg [L T-1 ~> m s-1]
     drag_rate, &    ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
-    drag_rate_J15, &  ! The MEKE spindown timescale due to bottom drag with the Jansen 2015 scheme.
-                    ! Unfortunately, as written the units seem inconsistent. [T-1 ~> s-1].
     del2MEKE, &     ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2].
     del4MEKE, &     ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2].
     LmixScale, &    ! Eddy mixing length [L ~> m].
@@ -179,8 +177,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
   if (.not.associated(MEKE)) call MOM_error(FATAL, &
          "MOM_MEKE: MEKE must be initialized before it is used.")
 
-  if ((CS%MEKE_damping > 0.0) .or. (CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) &
-      .or. CS%visc_drag) then
+  if ((CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) .or. CS%visc_drag) then
     use_drag_rate = .true.
   else
     use_drag_rate = .false.
@@ -236,12 +233,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       enddo
     endif
 
-    if (CS%MEKE_Cd_scale == 0.0 .and. .not. CS%visc_drag) then
-      !$OMP parallel do default(shared) private(ldamping)
-      do j=js,je ; do i=is,ie
-        drag_rate(i,j) = 0. ; drag_rate_J15(i,j) = 0.
-      enddo ; enddo
-    endif
 
     ! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow
     if (CS%visc_drag) then
@@ -370,6 +361,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
         drag_rate(i,j) = (US%L_to_Z*Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
                  cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
       enddo ; enddo
+    else
+      !$OMP parallel do default(shared)
+      do j=js,je ; do i=is,ie
+        drag_rate(i,j) = 0.
+      enddo ; enddo
     endif
 
     ! First stage of Strang splitting
@@ -538,23 +534,19 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
             drag_rate(i,j) = (US%L_to_Z*Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
                    cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
           enddo ; enddo
-          !$OMP parallel do default(shared)
-          do j=js,je ; do i=is,ie
-            ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
-            if (MEKE%MEKE(i,j) < 0.) ldamping = 0.
-            ! notice that the above line ensures a damping only if MEKE is positive,
-            ! while leaving MEKE unchanged if it is negative
-            MEKE%MEKE(i,j) =  MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
-            MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
-          enddo ; enddo
         endif
+        !$OMP parallel do default(shared)
+        do j=js,je ; do i=is,ie
+          ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
+          if (MEKE%MEKE(i,j) < 0.) ldamping = 0.
+          ! notice that the above line ensures a damping only if MEKE is positive,
+          ! while leaving MEKE unchanged if it is negative
+          MEKE%MEKE(i,j) =  MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
+          MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
+        enddo ; enddo
       endif
     endif ! MEKE_KH>=0
 
- !   do j=js,je ; do i=is,ie
- !     MEKE%MEKE(i,j) =  MAX(MEKE%MEKE(i,j),0.0)
- !   enddo ; enddo
-
     call cpu_clock_begin(CS%id_clock_pass)
     call do_group_pass(CS%pass_MEKE, G%Domain)
     call cpu_clock_end(CS%id_clock_pass)

diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -190,6 +190,7 @@ module MOM_hor_visc
   integer :: id_h_diffu  = -1, id_h_diffv      = -1
   integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1
   integer :: id_intz_diffu_2d = -1, id_intz_diffv_2d = -1
+  integer :: id_diffu_visc_rem = -1, id_diffv_visc_rem = -1
   integer :: id_Ah_h      = -1, id_Ah_q          = -1
   integer :: id_Kh_h      = -1, id_Kh_q          = -1
   integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1
@@ -282,6 +283,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
   real, allocatable, dimension(:,:,:) :: h_diffu ! h x diffu [H L T-2 ~> m2 s-2]
   real, allocatable, dimension(:,:,:) :: h_diffv ! h x diffv [H L T-2 ~> m2 s-2]
+  real, allocatable, dimension(:,:,:) :: diffu_visc_rem ! diffu x visc_rem_u [L T-2 ~> m s-2]
+  real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [L T-2 ~> m s-2]
 
   real, dimension(SZIB_(G),SZJB_(G)) :: &
     dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
@@ -1726,6 +1729,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     deallocate(h_diffv)
   endif
 
+  if (present(ADp) .and. (CS%id_diffu_visc_rem > 0)) then
+    allocate(diffu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
+    diffu_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      diffu_visc_rem(I,j,k) = diffu(I,j,k) * ADp%visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_diffu_visc_rem, diffu_visc_rem, CS%diag)
+    deallocate(diffu_visc_rem)
+  endif
+  if (present(ADp) .and. (CS%id_diffv_visc_rem > 0)) then
+    allocate(diffv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
+    diffv_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      diffv_visc_rem(i,J,k) = diffv(i,J,k) * ADp%visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_diffv_visc_rem, diffv_visc_rem, CS%diag)
+    deallocate(diffv_visc_rem)
+  endif
+
 end subroutine horizontal_viscosity
 
 !> Allocates space for and calculates static variables used by horizontal_viscosity().
@@ -2470,6 +2492,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
     call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
   endif
 
+  CS%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, Time, &
+      'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if ((CS%id_diffu_visc_rem > 0) .and. (present(ADp))) then
+    call safe_alloc_ptr(ADp%visc_rem_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
+  endif
+
+  CS%id_diffv_visc_rem = register_diag_field('ocean_model', 'diffv_visc_rem', diag%axesCvL, Time, &
+      'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if ((CS%id_diffv_visc_rem > 0) .and. (present(ADp))) then
+    call safe_alloc_ptr(ADp%visc_rem_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
+  endif
+
   if (CS%biharmonic) then
     CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time,    &
         'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, &