From 3acd777490645758a0844649dd3ad26a532adb1e Mon Sep 17 00:00:00 2001
From: Alistair Adcroft <Alistair.Adcroft@noaa.gov>
Date: Wed, 23 Dec 2020 16:01:44 -0500
Subject: [PATCH 01/24] Draft. Compiles but FPE in advection

---
 src/tracer/MOM_tracer_flow_control.F90 |  24 +-
 src/tracer/nw2_tracers.F90             | 319 +++++++++++++++++++++++++
 2 files changed, 338 insertions(+), 5 deletions(-)
 create mode 100644 src/tracer/nw2_tracers.F90

diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90
index 4c7c27c7e6..c1885644cc 100644
--- a/src/tracer/MOM_tracer_flow_control.F90
+++ b/src/tracer/MOM_tracer_flow_control.F90
@@ -62,6 +62,8 @@ module MOM_tracer_flow_control
 use boundary_impulse_tracer, only : boundary_impulse_tracer_column_physics, boundary_impulse_tracer_surface_state
 use boundary_impulse_tracer, only : boundary_impulse_stock, boundary_impulse_tracer_end
 use boundary_impulse_tracer, only : boundary_impulse_tracer_CS
+use nw2_tracers, only : nw2_tracers_CS, register_nw2_tracers, nw2_tracer_column_physics
+use nw2_tracers, only : initialize_nw2_tracers, nw2_tracers_end
 
 implicit none ; private
 
@@ -84,6 +86,7 @@ module MOM_tracer_flow_control
   logical :: use_pseudo_salt_tracer = .false.      !< If true, use the psuedo_salt tracer  package
   logical :: use_boundary_impulse_tracer = .false. !< If true, use the boundary impulse tracer package
   logical :: use_dyed_obc_tracer = .false.         !< If true, use the dyed OBC tracer package
+  logical :: use_nw2_tracers = .false.             !< If true, use the ideal age tracer package
   !>@{ Pointers to the control strucures for the tracer packages
   type(USER_tracer_example_CS), pointer :: USER_tracer_example_CSp => NULL()
   type(DOME_tracer_CS), pointer :: DOME_tracer_CSp => NULL()
@@ -98,6 +101,7 @@ module MOM_tracer_flow_control
   type(pseudo_salt_tracer_CS), pointer :: pseudo_salt_tracer_CSp => NULL()
   type(boundary_impulse_tracer_CS), pointer :: boundary_impulse_tracer_CSp => NULL()
   type(dyed_obc_tracer_CS), pointer :: dyed_obc_tracer_CSp => NULL()
+  type(nw2_tracers_CS), pointer :: nw2_tracers_CSp => NULL()
   !>@}
 end type tracer_flow_control_CS
 
@@ -206,6 +210,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
   call get_param(param_file, mdl, "USE_DYED_OBC_TRACER", CS%use_dyed_obc_tracer, &
                  "If true, use the dyed_obc_tracer tracer package.", &
                  default=.false.)
+  call get_param(param_file, mdl, "USE_NW2_TRACERS", CS%use_nw2_tracers, &
+                 "If true, use the NeverWorld2 tracers.", &
+                 default=.false.)
 
 !    Add other user-provided calls to register tracers for restarting here. Each
 !  tracer package registration call returns a logical false if it cannot be run
@@ -249,7 +256,8 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS)
   if (CS%use_dyed_obc_tracer) CS%use_dyed_obc_tracer = &
     register_dyed_obc_tracer(HI, GV, param_file, CS%dyed_obc_tracer_CSp, &
                              tr_Reg, restart_CS)
-
+  if (CS%use_nw2_tracers) CS%use_ideal_age = &
+    register_nw2_tracers(HI, GV, param_file,  CS%nw2_tracers_CSp, tr_Reg, restart_CS)
 
 end subroutine call_tracer_register
 
@@ -328,6 +336,8 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag
                                 sponge_CSp, tv)
   if (CS%use_dyed_obc_tracer) &
     call initialize_dyed_obc_tracer(restart, day, G, GV, h, diag, OBC, CS%dyed_obc_tracer_CSp)
+  if (CS%use_nw2_tracers) &
+    call initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS%nw2_tracers_CSp)
 
 end subroutine tracer_flow_control_init
 
@@ -485,8 +495,11 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV,
                                       G, GV, US, CS%dyed_obc_tracer_CSp, &
                                       evap_CFL_limit=evap_CFL_limit, &
                                       minimum_forcing_depth=minimum_forcing_depth)
-
-
+    if (CS%use_nw2_tracers) &
+      call nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+                                     G, GV, US, tv, CS%nw2_tracers_CSp, &
+                                     evap_CFL_limit=evap_CFL_limit, &
+                                     minimum_forcing_depth=minimum_forcing_depth)
   else ! Apply tracer surface fluxes using ea on the first layer
     if (CS%use_USER_tracer_example) &
       call tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
@@ -531,10 +544,10 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV,
     if (CS%use_dyed_obc_tracer) &
       call dyed_obc_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
                                       G, GV, US, CS%dyed_obc_tracer_CSp)
-
+    if (CS%use_nw2_tracers) call nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+                                                           G, GV, US, tv, CS%nw2_tracers_CSp)
   endif
 
-
 end subroutine call_tracer_column_fns
 
 !> This subroutine calls all registered tracer packages to enable them to
@@ -774,6 +787,7 @@ subroutine tracer_flow_control_end(CS)
   if (CS%use_pseudo_salt_tracer) call pseudo_salt_tracer_end(CS%pseudo_salt_tracer_CSp)
   if (CS%use_boundary_impulse_tracer) call boundary_impulse_tracer_end(CS%boundary_impulse_tracer_CSp)
   if (CS%use_dyed_obc_tracer) call dyed_obc_tracer_end(CS%dyed_obc_tracer_CSp)
+  if (CS%use_nw2_tracers) call nw2_tracers_end(CS%nw2_tracers_CSp)
 
   if (associated(CS)) deallocate(CS)
 end subroutine tracer_flow_control_end
diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
new file mode 100644
index 0000000000..e71387f9c7
--- /dev/null
+++ b/src/tracer/nw2_tracers.F90
@@ -0,0 +1,319 @@
+!> Ideal tracers designed to help diagnose a tracer diffusivity tensor in NeverWorld2
+module nw2_tracers
+
+! This file is part of MOM6. See LICENSE.md for the license.
+
+use MOM_diag_mediator, only : diag_ctrl
+use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
+use MOM_forcing_type, only : forcing
+use MOM_grid, only : ocean_grid_type
+use MOM_hor_index, only : hor_index_type
+use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc
+use MOM_restart, only : query_initialized, MOM_restart_CS
+use MOM_time_manager, only : time_type, time_type_to_real
+use MOM_tracer_registry, only : register_tracer, tracer_registry_type
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
+use MOM_unit_scaling, only : unit_scale_type
+use MOM_variables, only : surface, thermo_var_ptrs
+use MOM_verticalGrid, only : verticalGrid_type
+
+implicit none ; private
+
+#include <MOM_memory.h>
+
+public register_nw2_tracers
+public initialize_nw2_tracers
+public nw2_tracer_column_physics
+public nw2_tracers_end
+
+integer, parameter :: NTR_MAX = 20 !< the maximum number of tracers in this module.
+
+!> The control structure for the nw2_tracers package
+type, public :: nw2_tracers_CS ; private
+  integer :: ntr = NTR_MAX  !< The number of tracers that are actually used.
+  type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
+  type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry
+  real, pointer :: tr(:,:,:,:) => NULL()   !< The array of tracers used in this package, in g m-3?
+  real, dimension(NTR_MAX) :: restore_rate !< The exponential growth rate for restoration value [year-1].
+  type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
+                                             !! regulate the timing of diagnostic output.
+  type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart controls structure
+end type nw2_tracers_CS
+
+contains
+
+!> Register the NW2 tracer fields to be used with MOM.
+logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS)
+  type(hor_index_type),       intent(in) :: HI   !< A horizontal index type structure
+  type(verticalGrid_type),    intent(in) :: GV   !< The ocean's vertical grid structure
+  type(param_file_type),      intent(in) :: param_file !< A structure to parse for run-time parameters
+  type(nw2_tracers_CS),       pointer    :: CS !< The control structure returned by a previous
+                                               !! call to register_nw2_tracer.
+  type(tracer_registry_type), pointer    :: tr_Reg !< A pointer that is set to point to the control
+                                                  !! structure for the tracer advection and
+                                                  !! diffusion module
+  type(MOM_restart_CS),       pointer    :: restart_CS !< A pointer to the restart control structure
+
+! This include declares and sets the variable "version".
+#include "version_variable.h"
+  character(len=40)  :: mdl = "nw2_example" ! This module's name.
+  character(len=200) :: inputdir ! The directory where the input files are.
+  character(len=8)  :: var_name ! The variable's name.
+  real, pointer :: tr_ptr(:,:,:) => NULL()
+  logical :: do_nw2
+  integer :: isd, ied, jsd, jed, nz, m
+  type(vardesc) :: tr_desc ! Descriptions and metadata for the tracers
+  isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
+
+  if (associated(CS)) then
+    call MOM_error(FATAL, "register_nw2_tracer called with an "// &
+                          "associated control structure.")
+    return
+  endif
+  allocate(CS)
+
+  ! Read all relevant parameters and write them to the model log.
+  call log_version(param_file, mdl, version, "")
+
+  allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%tr(:,:,:,:) = 0.0
+
+  do m=1,CS%ntr
+    write(var_name(1:8),'(a6,i2.2)') 'tracer',m
+    tr_desc = var_desc(var_name, "1", "Ideal Tracer", caller=mdl)
+    ! This is needed to force the compiler not to do a copy in the registration
+    ! calls.  Curses on the designers and implementers of Fortran90.
+    tr_ptr => CS%tr(:,:,:,m)
+    call query_vardesc(tr_desc, name=var_name, caller="register_nw2_tracers")
+    ! Register the tracer for horizontal advection, diffusion, and restarts.
+    call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, &
+                         registry_diags=.true., restart_CS=restart_CS, mandatory=.true.)
+    if (m<=10) then
+      CS%restore_rate(m) = 1.0 / 6.0 ! 6 years
+    else
+      CS%restore_rate(m) = 1.0 ! 1 year
+    endif
+  enddo
+
+  CS%tr_Reg => tr_Reg
+  CS%restart_CSp => restart_CS
+  register_nw2_tracers = .true.
+end function register_nw2_tracers
+
+!> Sets the NW2 traces to their initial values and sets up the tracer output
+subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS)
+  logical,                            intent(in) :: restart !< .true. if the fields have already
+                                                         !! been read from a restart file.
+  type(time_type),            target, intent(in) :: day  !< Time of the start of the run.
+  type(ocean_grid_type),              intent(in) :: G    !< The ocean's grid structure
+  type(verticalGrid_type),            intent(in) :: GV   !< The ocean's vertical grid structure
+  type(unit_scale_type),              intent(in) :: US   !< A dimensional unit scaling type
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
+                                      intent(in) :: h    !< Layer thicknesses [H ~> m or kg m-2]
+  type(thermo_var_ptrs),              intent(in) :: tv   !< A structure pointing to various
+                                                         !! thermodynamic variables
+  type(diag_ctrl),            target, intent(in) :: diag !< A structure that is used to regulate
+                                                         !! diagnostic output.
+  type(nw2_tracers_CS),               pointer    :: CS !< The control structure returned by a previous
+                                                       !! call to register_nw2_tracer.
+  ! Local variables
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights
+  real :: rscl ! z* scaling factor
+  integer :: i, j, k, m
+
+  if (.not.associated(CS)) return
+
+  CS%Time => day
+  CS%diag => diag
+
+  ! Calculate z* interface positions
+  if (GV%Boussinesq) then
+    ! First calculate interface positions in z-space (m)
+    do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      eta(i,j,GV%ke+1) = - G%mask2dT(i,j) * G%bathyT(i,j)
+    enddo ; enddo
+    do k=GV%ke,1,-1 ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h(i,j,k) * GV%H_to_Z
+    enddo ; enddo ; enddo
+    ! Re-calculate for interface positions in z*-space (m)
+    do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      if (G%bathyT(i,j)>0.) then
+        rscl = G%bathyT(i,j) / ( eta(i,j,1) + G%bathyT(i,j) )
+        do K=GV%ke, 1, -1
+          eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h(i,j,k) * GV%H_to_Z * rscl
+        enddo
+      endif
+    enddo ; enddo
+  else
+    call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode")
+  endif
+
+  if (.not.restart) then
+    do m=1,CS%ntr
+      do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
+     !  if (G%mask2dT(i,j) < 0.5) then
+     !    CS%tr(i,j,k,m) = 0.
+     !  else
+          CS%tr(i,j,k,m) = nw2_tracer_dist(m, G, GV, eta, i, j, k)
+     !  endif
+      enddo ; enddo ; enddo
+    enddo ! Tracer loop
+  endif ! restart
+
+end subroutine initialize_nw2_tracers
+
+!> Applies diapycnal diffusion, aging and regeneration at the surface to the NW2 tracers
+subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, tv, CS, &
+              evap_CFL_limit, minimum_forcing_depth)
+  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure
+  type(verticalGrid_type), intent(in) :: GV   !< The ocean's vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
+                           intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
+                           intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
+                           intent(in) :: ea   !< an array to which the amount of fluid entrained
+                                              !! from the layer above during this call will be
+                                              !! added [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
+                           intent(in) :: eb   !< an array to which the amount of fluid entrained
+                                              !! from the layer below during this call will be
+                                              !! added [H ~> m or kg m-2].
+  type(forcing),           intent(in) :: fluxes !< A structure containing pointers to thermodynamic
+                                              !! and tracer forcing fields.  Unused fields have NULL ptrs.
+  real,                    intent(in) :: dt   !< The amount of time covered by this call [T ~> s]
+  type(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
+  type(thermo_var_ptrs),   intent(in) :: tv   !< A structure pointing to various thermodynamic variables
+  type(nw2_tracers_CS),    pointer    :: CS   !< The control structure returned by a previous
+                                              !! call to register_nw2_tracer.
+  real,          optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can
+                                              !! be fluxed out of the top layer in a timestep [nondim]
+  real,          optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
+                                              !! fluxes can be applied [H ~> m or kg m-2]
+!   This subroutine applies diapycnal diffusion and any other column
+! tracer physics or chemistry to the tracers from this file.
+! This is a simple example of a set of advected passive tracers.
+
+! The arguments to this subroutine are redundant in that
+!     h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
+  ! Local variables
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights
+  integer :: i, j, k, m
+  real :: dt_x_rate ! dt * restoring rate
+  real :: rscl ! z* scaling factor
+  real :: target_value ! tracer value
+
+! if (.not.associated(CS)) return
+
+  if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+    do m=1,CS%ntr
+      do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
+        h_work(i,j,k) = h_old(i,j,k)
+      enddo ; enddo ; enddo
+      call applyTracerBoundaryFluxesInOut(G, GV, CS%tr(:,:,:,m), dt, fluxes, h_work, &
+                                          evap_CFL_limit, minimum_forcing_depth)
+      call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+    enddo
+  else
+    do m=1,CS%ntr
+      call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+    enddo
+  endif
+
+  ! Calculate z* interface positions
+  if (GV%Boussinesq) then
+    ! First calculate interface positions in z-space (m)
+    do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      eta(i,j,GV%ke+1) = - G%mask2dT(i,j) * G%bathyT(i,j)
+    enddo ; enddo
+    do k=GV%ke,1,-1 ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h_new(i,j,k) * GV%H_to_Z
+    enddo ; enddo ; enddo
+    ! Re-calculate for interface positions in z*-space (m)
+    do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      if (G%bathyT(i,j)>0.) then
+        rscl = G%bathyT(i,j) / ( eta(i,j,1) + G%bathyT(i,j) )
+        do K=GV%ke, 1, -1
+          eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h_new(i,j,k) * GV%H_to_Z * rscl
+        enddo
+      endif
+    enddo ; enddo
+  else
+    call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode")
+  endif
+
+  do m=1,CS%ntr
+    dt_x_rate = dt * ( CS%restore_rate(m)*((365.0*86400.0)*US%s_to_T) )
+!$OMP parallel do default(private) shared(CS,G,dt,dt_x_rate)
+    do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
+      target_value = nw2_tracer_dist(m, G, GV, eta, i, j, k)
+      CS%tr(i,j,k,m) = CS%tr(i,j,k,m) + G%mask2dT(i,j) * dt_x_rate * ( target_value - CS%tr(i,j,k,m) )
+    enddo ; enddo ; enddo
+  enddo
+
+end subroutine nw2_tracer_column_physics
+
+!> The target value of a NeverWorld2 tracer label m at non-dimensional
+!! position x=lon/Lx, y=lat/Ly, z=eta/H
+real function nw2_tracer_dist(m, G, GV, eta, i, j, k)
+  integer, intent(in) :: m !< Indicates the NW2 tracer
+  type(ocean_grid_type),   intent(in) :: G   !< The ocean's grid structure
+  type(verticalGrid_type), intent(in) :: GV  !< The ocean's vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),0:SZK_(G)), &
+                           intent(in) :: eta !< Interface position [m]
+  integer, intent(in) :: i !< Cell index i
+  integer, intent(in) :: j !< Cell index j
+  integer, intent(in) :: k !< Layer index k
+  ! Local variables
+  real :: pi ! 3.1415...
+  real :: x, y, z ! non-dimensional positions
+  pi = 2.*acos(0.)
+  x = ( G%geolonT(i,j) - G%west_lon ) / G%len_lon ! 0 ... 1
+  y = -G%geolatT(i,j) / G%south_lat ! -1 ... 1
+  z = - 0.5 * ( eta(i,j,K-1) + eta(i,j,K) ) / GV%max_depth ! 0 ... 1
+  select case ( mod(m-1,10) )
+  case (0) ! y/L
+    nw2_tracer_dist = y
+  case (1) ! -z/L
+    nw2_tracer_dist = -z
+  case (2) ! cos(2 pi x/L)
+    nw2_tracer_dist = cos( 2.0 * pi * x )
+  case (3) ! sin(2 pi x/L)
+    nw2_tracer_dist = sin( 2.0 * pi * x )
+  case (4) ! sin(4 pi x/L)
+    nw2_tracer_dist = sin( 2.0 * pi * mod( 2.0*x, 1.0 ) )
+  case (5) ! sin(pi y/L)
+    nw2_tracer_dist = sin( pi * y )
+  case (6) ! cos(2 pi y/L)
+    nw2_tracer_dist = cos( 2.0 * pi * y )
+  case (7) ! sin(2 pi y/L)
+    nw2_tracer_dist = sin( 2.0 * pi * y )
+  case (8) ! cos(pi z/H)
+    nw2_tracer_dist = cos( pi * z )
+  case (9) ! sin(pi z/H)
+    nw2_tracer_dist = sin( pi * z )
+  case default
+    stop 'This should not happen. Died in nw2_tracer_dist()!'
+  end select
+  nw2_tracer_dist = nw2_tracer_dist * G%mask2dT(i,j)
+end function nw2_tracer_dist
+
+!> Deallocate any memory associated with this tracer package
+subroutine nw2_tracers_end(CS)
+  type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous
+                                      !! call to register_nw2_tracers.
+
+  integer :: m
+
+  if (associated(CS)) then
+    if (associated(CS%tr)) deallocate(CS%tr)
+    deallocate(CS)
+  endif
+end subroutine nw2_tracers_end
+
+!> \namespace nw2_tracers
+!!
+!!  TBD
+
+end module nw2_tracers

From ca517073610271080c5be107d7b36c63b38295fd Mon Sep 17 00:00:00 2001
From: Alistair Adcroft <Alistair.Adcroft@noaa.gov>
Date: Sun, 3 Jan 2021 17:19:23 -0500
Subject: [PATCH 02/24] Fix FPE problems with NW2 tracers

- I had a "*" in place of a "/" for calculating a non-dimensional
  quantity that led to OOB numbers.
---
 src/tracer/nw2_tracers.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
index e71387f9c7..6b2bdf2c0d 100644
--- a/src/tracer/nw2_tracers.F90
+++ b/src/tracer/nw2_tracers.F90
@@ -244,7 +244,7 @@ subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US
   endif
 
   do m=1,CS%ntr
-    dt_x_rate = dt * ( CS%restore_rate(m)*((365.0*86400.0)*US%s_to_T) )
+    dt_x_rate = dt / ( CS%restore_rate(m)*((365.0*86400.0)*US%s_to_T) )
 !$OMP parallel do default(private) shared(CS,G,dt,dt_x_rate)
     do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
       target_value = nw2_tracer_dist(m, G, GV, eta, i, j, k)

From 77c67fcfb416e10b49eb14b2a2efe5be77032271 Mon Sep 17 00:00:00 2001
From: Alistair Adcroft <Alistair.Adcroft@noaa.gov>
Date: Sun, 3 Jan 2021 17:33:44 -0500
Subject: [PATCH 03/24] Allow NW2 tracesr to be added mid-run

- I'd previously made the NW2 tracers mandatory in the restart field
  which prohibits adding tracers on a restart. This commit allows a
  restart without tracers to be used.
---
 src/tracer/nw2_tracers.F90 | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
index 6b2bdf2c0d..b3d751d084 100644
--- a/src/tracer/nw2_tracers.F90
+++ b/src/tracer/nw2_tracers.F90
@@ -87,7 +87,7 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
     call query_vardesc(tr_desc, name=var_name, caller="register_nw2_tracers")
     ! Register the tracer for horizontal advection, diffusion, and restarts.
     call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, &
-                         registry_diags=.true., restart_CS=restart_CS, mandatory=.true.)
+                         registry_diags=.true., restart_CS=restart_CS, mandatory=.false.)
     if (m<=10) then
       CS%restore_rate(m) = 1.0 / 6.0 ! 6 years
     else
@@ -148,14 +148,13 @@ subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS)
     call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode")
   endif
 
-  if (.not.restart) then
+  ! Initialize only if this is not a restart or we are using a restart
+  ! in which the tracers were not present
+  if ((.not.restart) .or. &
+      (.not. query_initialized(CS%tr(:,:,:,m),'nw2_tracer',CS%restart_CSp))) then
     do m=1,CS%ntr
       do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
-     !  if (G%mask2dT(i,j) < 0.5) then
-     !    CS%tr(i,j,k,m) = 0.
-     !  else
           CS%tr(i,j,k,m) = nw2_tracer_dist(m, G, GV, eta, i, j, k)
-     !  endif
       enddo ; enddo ; enddo
     enddo ! Tracer loop
   endif ! restart

From b333739df4e3c95f0681078ae1db9a3e63da338c Mon Sep 17 00:00:00 2001
From: Alistair Adcroft <Alistair.Adcroft@noaa.gov>
Date: Sun, 3 Jan 2021 17:35:16 -0500
Subject: [PATCH 04/24] Cleanup NW2 registration

- I noticed an unnecessary function call which is now removed.
---
 src/tracer/nw2_tracers.F90 | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
index b3d751d084..3199a5b271 100644
--- a/src/tracer/nw2_tracers.F90
+++ b/src/tracer/nw2_tracers.F90
@@ -9,7 +9,7 @@ module nw2_tracers
 use MOM_forcing_type, only : forcing
 use MOM_grid, only : ocean_grid_type
 use MOM_hor_index, only : hor_index_type
-use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc
+use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc
 use MOM_restart, only : query_initialized, MOM_restart_CS
 use MOM_time_manager, only : time_type, time_type_to_real
 use MOM_tracer_registry, only : register_tracer, tracer_registry_type
@@ -84,7 +84,6 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
     ! This is needed to force the compiler not to do a copy in the registration
     ! calls.  Curses on the designers and implementers of Fortran90.
     tr_ptr => CS%tr(:,:,:,m)
-    call query_vardesc(tr_desc, name=var_name, caller="register_nw2_tracers")
     ! Register the tracer for horizontal advection, diffusion, and restarts.
     call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, &
                          registry_diags=.true., restart_CS=restart_CS, mandatory=.false.)

From adb3554add3c787093e8aeb0f7327f2d0feaede9 Mon Sep 17 00:00:00 2001
From: Alistair Adcroft <Alistair.Adcroft@noaa.gov>
Date: Sun, 3 Jan 2021 18:00:20 -0500
Subject: [PATCH 05/24] Fixed NW2 tracer restarts

- An inverted do/if loop was using uninitialized variables and
  was overwriting restart data.
---
 src/tracer/nw2_tracers.F90 | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
index 3199a5b271..5385771445 100644
--- a/src/tracer/nw2_tracers.F90
+++ b/src/tracer/nw2_tracers.F90
@@ -118,6 +118,7 @@ subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS)
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights
   real :: rscl ! z* scaling factor
+  character(len=8)  :: var_name ! The variable's name.
   integer :: i, j, k, m
 
   if (.not.associated(CS)) return
@@ -147,16 +148,17 @@ subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS)
     call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode")
   endif
 
-  ! Initialize only if this is not a restart or we are using a restart
-  ! in which the tracers were not present
-  if ((.not.restart) .or. &
-      (.not. query_initialized(CS%tr(:,:,:,m),'nw2_tracer',CS%restart_CSp))) then
-    do m=1,CS%ntr
+  do m=1,CS%ntr
+    ! Initialize only if this is not a restart or we are using a restart
+    ! in which the tracers were not present
+    write(var_name(1:8),'(a6,i2.2)') 'tracer',m
+    if ((.not.restart) .or. &
+        (.not. query_initialized(CS%tr(:,:,:,m),var_name,CS%restart_CSp))) then
       do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
           CS%tr(i,j,k,m) = nw2_tracer_dist(m, G, GV, eta, i, j, k)
       enddo ; enddo ; enddo
-    enddo ! Tracer loop
-  endif ! restart
+    endif ! restart
+  enddo ! Tracer loop
 
 end subroutine initialize_nw2_tracers
 

From d20ecf7b6f382f6e62a7aee940ca7bf79d5c6263 Mon Sep 17 00:00:00 2001
From: Alistair Adcroft <Alistair.Adcroft@noaa.gov>
Date: Mon, 22 Feb 2021 12:56:55 -0500
Subject: [PATCH 06/24] Remove use of parameter for NTR_MAX in NW2

- The number of tracers was hard coded via an integer parameter which
  is unnecessarily restrictive.
---
 src/tracer/nw2_tracers.F90 | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
index 5385771445..5417373b24 100644
--- a/src/tracer/nw2_tracers.F90
+++ b/src/tracer/nw2_tracers.F90
@@ -27,15 +27,13 @@ module nw2_tracers
 public nw2_tracer_column_physics
 public nw2_tracers_end
 
-integer, parameter :: NTR_MAX = 20 !< the maximum number of tracers in this module.
-
 !> The control structure for the nw2_tracers package
 type, public :: nw2_tracers_CS ; private
-  integer :: ntr = NTR_MAX  !< The number of tracers that are actually used.
+  integer :: ntr = 0  !< The number of tracers that are actually used.
   type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
   type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry
   real, pointer :: tr(:,:,:,:) => NULL()   !< The array of tracers used in this package, in g m-3?
-  real, dimension(NTR_MAX) :: restore_rate !< The exponential growth rate for restoration value [year-1].
+  real, allocatable , dimension(:) :: restore_rate !< The exponential growth rate for restoration value [year-1].
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                              !! regulate the timing of diagnostic output.
   type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart controls structure
@@ -57,7 +55,7 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
 
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
-  character(len=40)  :: mdl = "nw2_example" ! This module's name.
+  character(len=40)  :: mdl = "nw2_tracers" ! This module's name.
   character(len=200) :: inputdir ! The directory where the input files are.
   character(len=8)  :: var_name ! The variable's name.
   real, pointer :: tr_ptr(:,:,:) => NULL()
@@ -76,7 +74,9 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
   ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
 
+  CS%ntr=20
   allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%tr(:,:,:,:) = 0.0
+  allocate(CS%restore_rate(CS%ntr))
 
   do m=1,CS%ntr
     write(var_name(1:8),'(a6,i2.2)') 'tracer',m

From 8f6ec58bc0246e923a5a6430a2524a5bb693fbdc Mon Sep 17 00:00:00 2001
From: Alistair Adcroft <Alistair.Adcroft@noaa.gov>
Date: Mon, 22 Feb 2021 15:08:14 -0500
Subject: [PATCH 07/24] Replace NW2 2x10 tracers with 3x3

- The first interation of NW2 tracers implemented 20 tracesr with 10
  different spatial structures and 2 different restoration time scales.
  This has been replaced with 9 tracers with 3 spatial, corresponding
  to sin(2x), y, z, and 3 time scales.
- Number of tracer groups and time scales are now run-time configurable.
---
 src/tracer/nw2_tracers.F90 | 51 +++++++++++++++++---------------------
 1 file changed, 23 insertions(+), 28 deletions(-)

diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
index 5417373b24..e75c5c5d38 100644
--- a/src/tracer/nw2_tracers.F90
+++ b/src/tracer/nw2_tracers.F90
@@ -60,7 +60,9 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
   character(len=8)  :: var_name ! The variable's name.
   real, pointer :: tr_ptr(:,:,:) => NULL()
   logical :: do_nw2
-  integer :: isd, ied, jsd, jed, nz, m
+  integer :: isd, ied, jsd, jed, nz, m, ig
+  integer :: n_groups ! Number of groups of three tracers (i.e. # tracers/3)
+  real, allocatable, dimension(:) :: timescale_in_days
   type(vardesc) :: tr_desc ! Descriptions and metadata for the tracers
   isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
 
@@ -73,8 +75,18 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
 
   ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
-
-  CS%ntr=20
+  call get_param(param_file, mdl, "NW2_TRACER_GROUPS", n_groups, &
+                 "The number of tracer groups where a group is of three tracers "//&
+                 "initialized and restored to sin(x), y and z, respectively. Each "//&
+                 "group is restored with an independent restoration rate.", &
+                 default=3)
+  allocate(timescale_in_days(n_groups))
+  timescale_in_days = (/365., 730., 1460./)
+  call get_param(param_file, mdl, "NW2_TRACER_RESTORE_TIMESCALE", timescale_in_days, &
+                 "A list of timescales, one for each tracer group.", &
+                 units="days")
+
+  CS%ntr = 3 * n_groups
   allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%tr(:,:,:,:) = 0.0
   allocate(CS%restore_rate(CS%ntr))
 
@@ -87,11 +99,8 @@ logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS
     ! Register the tracer for horizontal advection, diffusion, and restarts.
     call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, &
                          registry_diags=.true., restart_CS=restart_CS, mandatory=.false.)
-    if (m<=10) then
-      CS%restore_rate(m) = 1.0 / 6.0 ! 6 years
-    else
-      CS%restore_rate(m) = 1.0 ! 1 year
-    endif
+    ig = int( (m+2)/3 ) ! maps (1,2,3)->1, (4,5,6)->2, ...
+    CS%restore_rate(m) = 1.0 / ( timescale_in_days(ig) * 86400.0 )
   enddo
 
   CS%tr_Reg => tr_Reg
@@ -244,7 +253,7 @@ subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US
   endif
 
   do m=1,CS%ntr
-    dt_x_rate = dt / ( CS%restore_rate(m)*((365.0*86400.0)*US%s_to_T) )
+    dt_x_rate = ( dt * CS%restore_rate(m) ) * US%T_to_s
 !$OMP parallel do default(private) shared(CS,G,dt,dt_x_rate)
     do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
       target_value = nw2_tracer_dist(m, G, GV, eta, i, j, k)
@@ -272,27 +281,13 @@ real function nw2_tracer_dist(m, G, GV, eta, i, j, k)
   x = ( G%geolonT(i,j) - G%west_lon ) / G%len_lon ! 0 ... 1
   y = -G%geolatT(i,j) / G%south_lat ! -1 ... 1
   z = - 0.5 * ( eta(i,j,K-1) + eta(i,j,K) ) / GV%max_depth ! 0 ... 1
-  select case ( mod(m-1,10) )
-  case (0) ! y/L
+  select case ( mod(m-1,3) )
+  case (0) ! sin(2 pi x/L)
+    nw2_tracer_dist = sin( 2.0 * pi * x )
+  case (1) ! y/L
     nw2_tracer_dist = y
-  case (1) ! -z/L
+  case (2) ! -z/L
     nw2_tracer_dist = -z
-  case (2) ! cos(2 pi x/L)
-    nw2_tracer_dist = cos( 2.0 * pi * x )
-  case (3) ! sin(2 pi x/L)
-    nw2_tracer_dist = sin( 2.0 * pi * x )
-  case (4) ! sin(4 pi x/L)
-    nw2_tracer_dist = sin( 2.0 * pi * mod( 2.0*x, 1.0 ) )
-  case (5) ! sin(pi y/L)
-    nw2_tracer_dist = sin( pi * y )
-  case (6) ! cos(2 pi y/L)
-    nw2_tracer_dist = cos( 2.0 * pi * y )
-  case (7) ! sin(2 pi y/L)
-    nw2_tracer_dist = sin( 2.0 * pi * y )
-  case (8) ! cos(pi z/H)
-    nw2_tracer_dist = cos( pi * z )
-  case (9) ! sin(pi z/H)
-    nw2_tracer_dist = sin( pi * z )
   case default
     stop 'This should not happen. Died in nw2_tracer_dist()!'
   end select

From 7a650261c2765d3c60ccd931a1ed609914842079 Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Sat, 24 Jul 2021 13:46:41 -0600
Subject: [PATCH 08/24] Implement visc remnant versions of PFu, CAu,
 u_accel_bt, diffu

---
 src/core/MOM_barotropic.F90                   | 12 ++-
 src/core/MOM_dynamics_split_RK2.F90           | 85 +++++++++++++++++++
 src/core/MOM_variables.F90                    | 12 ++-
 .../lateral/MOM_hor_visc.F90                  | 36 ++++++++
 4 files changed, 143 insertions(+), 2 deletions(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index a8262608f8..42f91b2d8e 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -2685,7 +2685,17 @@ 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)
     call complete_group_pass(CS%pass_ubta_uhbta, G%Domain)
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 07a302b2b0..c55dbce684 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -163,10 +163,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
@@ -175,6 +177,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
@@ -353,6 +356,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [H L T-2 ~> m2 s-2].
   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].
 
@@ -1102,6 +1111,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%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%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%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%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%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%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)
@@ -1605,6 +1669,27 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       'Depth-integral of Barotropic Anomaly Meridional Acceleration', &
       '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', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  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', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  
+  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', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  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', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  
+  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', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  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', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
 
   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)
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index f966ab2ad2..c81d1f6f8d 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -163,16 +163,24 @@ module MOM_variables
   real, pointer, dimension(:,:,:) :: &
     diffu => NULL(), &     !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
     diffv => NULL(), &     !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
+    diffu_visc_rem => NULL(), &     !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+    diffv_visc_rem => NULL(), &     !< Meridional acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     CAu => NULL(), &       !< Zonal Coriolis and momentum advection accelerations [L T-2 ~> m s-2]
     CAv => NULL(), &       !< Meridional Coriolis and momentum advection accelerations [L T-2 ~> m s-2]
+    CAu_visc_rem => NULL(), &       !< Zonal Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+    CAv_visc_rem => NULL(), &       !< Meridional Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     PFu => NULL(), &       !< Zonal acceleration due to pressure forces [L T-2 ~> m s-2]
     PFv => NULL(), &       !< Meridional acceleration due to pressure forces [L T-2 ~> m s-2]
+    PFu_visc_rem => NULL(), &       !< Zonal acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+    PFv_visc_rem => NULL(), &       !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2]
     dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2]
     du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
     dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
     u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2]
-    v_accel_bt => NULL()   !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
+    v_accel_bt => NULL(), & !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
+    u_accel_bt_visc_rem => NULL(), &!< Pointer to the zonal barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+    v_accel_bt_visc_rem => NULL()   !< Pointer to the meridional barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2]
   real, pointer, dimension(:,:,:) :: du_other => NULL()
                            !< Zonal velocity changes due to any other processes that are
                            !! not due to any explicit accelerations [L T-1 ~> m s-1].
@@ -191,6 +199,8 @@ 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/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index 2324970254..bafce4b209 100644
--- 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
@@ -280,6 +281,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 ~> m2 s-2]
+  real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [L T-2 ~> m2 s-2]
 
   real, dimension(SZIB_(G),SZJB_(G)) :: &
     dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
@@ -1703,6 +1706,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     call post_data(CS%id_h_diffv, h_diffv, CS%diag)
     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
 
@@ -2447,6 +2469,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
   if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then
     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', 'm2 s-2', &
+      conversion=GV%H_to_m*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', 'm2 s-2', &
+      conversion=GV%H_to_m*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,    &

From 2601b47285ebcc9ed6fc3742a960f43c92d38809 Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Sat, 24 Jul 2021 16:03:01 -0600
Subject: [PATCH 09/24] Add d[uv]_dt_visc_rem diagnostics

---
 src/core/MOM_variables.F90                    |  4 ++
 src/diagnostics/MOM_diagnostics.F90           |  7 +++
 .../vertical/MOM_vert_friction.F90            | 53 +++++++++++++++++++
 3 files changed, 64 insertions(+)

diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index c81d1f6f8d..7272a5994f 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -161,6 +161,8 @@ module MOM_variables
 
   ! Each of the following fields has nz layers.
   real, pointer, dimension(:,:,:) :: &
+    du_dt => NULL(), &     !< Zonal acceleration [L T-2 ~> m s-2]
+    dv_dt => NULL(), &     !< Meridional acceleration [L T-2 ~> m s-2]
     diffu => NULL(), &     !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
     diffv => NULL(), &     !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
     diffu_visc_rem => NULL(), &     !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
@@ -175,6 +177,8 @@ module MOM_variables
     PFv_visc_rem => NULL(), &       !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2]
     dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2]
+    du_dt_visc_rem => NULL(), &!< Zonal acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+    dv_dt_visc_rem => NULL(), &!< Meridional acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
     dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
     u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2]
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 1d59969011..00ed7682d7 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -290,6 +290,13 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
     if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h = diag_pre_sync%h_state)
 
     if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state)
+    
+    do k=1,nz ; do j=js,je ; do I=is-1,ie
+        ADp%du_dt(I,j,k) = CS%du_dt(I,j,k)
+    enddo ; enddo ; enddo
+    do k=1,nz ; do J=js-1,je ; do i=is,ie
+        ADp%dv_dt(i,J,k) = CS%dv_dt(i,J,k)
+    enddo ; enddo ; enddo
 
     !! Diagnostics for terms multiplied by fractional thicknesses
 
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 081179eb41..d958188a00 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -127,6 +127,7 @@ module MOM_vert_friction
   ! integer :: id_hf_du_dt_visc    = -1, id_hf_dv_dt_visc    = -1
   integer :: id_h_du_dt_visc    = -1, id_h_dv_dt_visc    = -1
   integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1
+  integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1
   !>@}
 
   type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure
@@ -216,6 +217,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
 
   real, allocatable, dimension(:,:,:) :: h_du_dt_visc ! h x du_dt_visc [H L T-2 ~> m2 s-2]
   real, allocatable, dimension(:,:,:) :: h_dv_dt_visc ! h x dv_dt_visc [H L T-2 ~> m2 s-2]
+  
+  real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u) [L T-2 ~> m2 s-2]
+  real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v) [L T-2 ~> m2 s-2]
 
   logical :: do_i(SZIB_(G))
   logical :: DoStokesMixing
@@ -335,6 +339,11 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq
       ADp%du_dt_visc(I,j,k) = (u(I,j,k) - ADp%du_dt_visc(I,j,k))*Idt
     enddo ; enddo ; endif
+    
+    !if (associated(ADp%du_dt_visc_rem)) then ; do k=1,nz ; do I=Isq,Ieq
+    !  ADp%du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + &
+    !                            (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k)
+    !enddo ; enddo ; endif
 
     if (associated(visc%taux_shelf)) then ; do I=Isq,Ieq
       visc%taux_shelf(I,j) = -GV%Rho0*CS%a1_shelf_u(I,j)*u(I,j,1) ! - u_shelf?
@@ -416,6 +425,11 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
       ADp%dv_dt_visc(i,J,k) = (v(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt
     enddo ; enddo ; endif
+    
+    !if (associated(ADp%dv_dt_visc_rem)) then ; do k=1,nz ; do I=Isq,Ieq
+    !    ADp%dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + &
+    !                            (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k)
+    !enddo ; enddo ; endif
 
     if (associated(visc%tauy_shelf)) then ; do i=is,ie
       visc%tauy_shelf(i,J) = -GV%Rho0*CS%a1_shelf_v(i,J)*v(i,J,1) ! - v_shelf?
@@ -521,6 +535,27 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     call post_data(CS%id_h_dv_dt_visc, h_dv_dt_visc, CS%diag)
     deallocate(h_dv_dt_visc)
   endif
+  
+  if (CS%id_du_dt_visc_rem > 0) then
+    allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
+    du_dt_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + &
+                                (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag)
+    deallocate(du_dt_visc_rem)
+  endif
+  if (CS%id_dv_dt_visc_rem > 0) then
+    allocate(dv_dt_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
+    dv_dt_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + &
+                                (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag)
+    deallocate(dv_dt_visc_rem)
+  endif
 
 end subroutine vertvisc
 
@@ -1876,6 +1911,24 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
     call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
+  
+  CS%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, &
+      'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if (CS%id_du_dt_visc_rem > 0) then
+    call safe_alloc_ptr(ADp%du_dt,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
+  endif
+
+  CS%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, &
+      'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if (CS%id_dv_dt_visc_rem > 0) then
+    call safe_alloc_ptr(ADp%dv_dt,isd,ied,JsdB,JedB,nz)
+    call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
+    call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+  endif
 
   if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &
     call PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS%PointAccel_CSp)

From 20b8bfc9e6cad5fd773e0d1bffe64a6d13e7bf84 Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Sun, 25 Jul 2021 16:14:11 -0600
Subject: [PATCH 10/24] Attempt to implement KE budget modified by visc rem

---
 src/core/MOM_variables.F90          |  24 ++--
 src/diagnostics/MOM_diagnostics.F90 | 178 +++++++++++++++++++++++++++-
 2 files changed, 185 insertions(+), 17 deletions(-)

diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 7272a5994f..2f9c3cf643 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -165,26 +165,16 @@ module MOM_variables
     dv_dt => NULL(), &     !< Meridional acceleration [L T-2 ~> m s-2]
     diffu => NULL(), &     !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
     diffv => NULL(), &     !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
-    diffu_visc_rem => NULL(), &     !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-    diffv_visc_rem => NULL(), &     !< Meridional acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     CAu => NULL(), &       !< Zonal Coriolis and momentum advection accelerations [L T-2 ~> m s-2]
     CAv => NULL(), &       !< Meridional Coriolis and momentum advection accelerations [L T-2 ~> m s-2]
-    CAu_visc_rem => NULL(), &       !< Zonal Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-    CAv_visc_rem => NULL(), &       !< Meridional Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     PFu => NULL(), &       !< Zonal acceleration due to pressure forces [L T-2 ~> m s-2]
     PFv => NULL(), &       !< Meridional acceleration due to pressure forces [L T-2 ~> m s-2]
-    PFu_visc_rem => NULL(), &       !< Zonal acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-    PFv_visc_rem => NULL(), &       !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2]
     dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2]
-    du_dt_visc_rem => NULL(), &!< Zonal acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-    dv_dt_visc_rem => NULL(), &!< Meridional acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
     du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
     dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
     u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2]
-    v_accel_bt => NULL(), & !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
-    u_accel_bt_visc_rem => NULL(), &!< Pointer to the zonal barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-    v_accel_bt_visc_rem => NULL()   !< Pointer to the meridional barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+    v_accel_bt => NULL()  !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
   real, pointer, dimension(:,:,:) :: du_other => NULL()
                            !< Zonal velocity changes due to any other processes that are
                            !! not due to any explicit accelerations [L T-1 ~> m s-1].
@@ -205,6 +195,18 @@ module MOM_variables
 
   real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
   real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
+  
+  real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+
 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
index 00ed7682d7..28b5764079 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -108,7 +108,16 @@ module MOM_diagnostics
     KE_adv     => NULL(), & !< KE source from along-layer advection [H L2 T-3 ~> m3 s-3]
     KE_visc    => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3]
     KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3]
-    KE_dia     => NULL()    !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]
+    KE_dia     => NULL(), & !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]
+  
+  ! The following arrays hold diagnostics in the modified layer-integrated energy budget.
+  ! Modification is through using the visc_rem_[uv]-filtered momentum equation
+    PE_to_KE_visc_rem  => NULL(), &  !< potential energy to KE term [m3 s-3]
+    KE_BT_visc_rem     => NULL(), &  !< barotropic contribution to KE term [m3 s-3]
+    KE_CorAdv_visc_rem => NULL(), &  !< KE source from the combined Coriolis and
+                            !! advection terms [H L2 T-3 ~> m3 s-3].
+    KE_visc_rem    => NULL(), &  !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3]
+    KE_horvisc_rem => NULL() !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3]
 
   !>@{ Diagnostic IDs
   integer :: id_u   = -1,   id_v   = -1, id_h = -1
@@ -124,6 +133,10 @@ module MOM_diagnostics
   integer :: id_KE_Coradv      = -1
   integer :: id_KE_adv         = -1, id_KE_visc        = -1
   integer :: id_KE_horvisc     = -1, id_KE_dia         = -1
+  integer :: id_PE_to_KE_visc_rem = -1, id_KE_BT_visc_rem = -1
+  integer :: id_KE_Coradv_visc_rem = -1
+  integer :: id_KE_visc_rem = -1
+  integer :: id_KE_horvisc_rem = -1
   integer :: id_uh_Rlay        = -1, id_vh_Rlay        = -1
   integer :: id_uhGM_Rlay      = -1, id_vhGM_Rlay      = -1
   integer :: id_h_Rlay         = -1, id_Rd1            = -1
@@ -1055,7 +1068,10 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
   if (.not.G%symmetric) then
     if (associated(CS%dKE_dt) .OR. associated(CS%PE_to_KE) .OR. associated(CS%KE_BT) .OR. &
         associated(CS%KE_CorAdv) .OR. associated(CS%KE_adv) .OR. associated(CS%KE_visc) .OR. &
-        associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) ) then
+        associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) .OR. &
+        associated(CS%PE_to_KE_visc_rem) .OR. associated(CS%KE_BT_visc_rem) .OR. &
+        associated(CS%KE_CorAdv_visc_rem) .OR. associated(CS%KE_visc_rem) .OR. &
+        associated(CS%KE_horvisc_rem)) then
       call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East)
     endif
   endif
@@ -1223,6 +1239,100 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
     enddo
     if (CS%id_KE_dia > 0) call post_data(CS%id_KE_dia, CS%KE_dia, CS%diag)
   endif
+  
+  if (associated(CS%PE_to_KE_visc_rem)) then
+    do k=1,nz
+      do j=js,je ; do I=Isq,Ieq
+        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%PFu_visc_rem(I,j,k)
+      enddo ; enddo
+      do J=Jsq,Jeq ; do i=is,ie
+        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%PFv_visc_rem(i,J,k)
+      enddo ; enddo
+      if (.not.G%symmetric) &
+         call do_group_pass(CS%pass_KE_uv, G%domain)
+      do j=js,je ; do i=is,ie
+        CS%PE_to_KE_visc_rem(i,j,k) = 0.5 * G%IareaT(i,j) &
+            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
+      enddo ; enddo
+    enddo
+    if (CS%id_PE_to_KE_visc_rem > 0) call post_data(CS%id_PE_to_KE_visc_rem, CS%PE_to_KE_visc_rem, CS%diag)
+  endif
+
+  if (associated(CS%KE_BT_visc_rem)) then
+    do k=1,nz
+      do j=js,je ; do I=Isq,Ieq
+        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%u_accel_bt_visc_rem(I,j,k)
+      enddo ; enddo
+      do J=Jsq,Jeq ; do i=is,ie
+        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%v_accel_bt_visc_rem(i,J,k)
+      enddo ; enddo
+      if (.not.G%symmetric) &
+         call do_group_pass(CS%pass_KE_uv, G%domain)
+      do j=js,je ; do i=is,ie
+        CS%KE_BT_visc_rem(i,j,k) = 0.5 * G%IareaT(i,j) &
+            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
+      enddo ; enddo
+    enddo
+    if (CS%id_KE_BT_visc_rem > 0) call post_data(CS%id_KE_BT_visc_rem, CS%KE_BT_visc_rem, CS%diag)
+  endif
+
+  if (associated(CS%KE_CorAdv_visc_rem)) then
+    do k=1,nz
+      do j=js,je ; do I=Isq,Ieq
+        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%CAu_visc_rem(I,j,k)
+      enddo ; enddo
+      do J=Jsq,Jeq ; do i=is,ie
+        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%CAv_visc_rem(i,J,k)
+      enddo ; enddo
+      do j=js,je ; do i=is,ie
+        KE_h(i,j) = -CS%KE(i,j,k) * G%IareaT(i,j) &
+            * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k))
+      enddo ; enddo
+      if (.not.G%symmetric) &
+         call do_group_pass(CS%pass_KE_uv, G%domain)
+      do j=js,je ; do i=is,ie
+        CS%KE_CorAdv_visc_rem(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) &
+            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
+      enddo ; enddo
+    enddo
+    if (CS%id_KE_CorAdv_visc_rem > 0) call post_data(CS%id_KE_Coradv_visc_rem, CS%KE_Coradv_visc_rem, CS%diag)
+  endif
+
+  if (associated(CS%KE_visc_rem)) then
+    do k=1,nz
+      do j=js,je ; do I=Isq,Ieq
+        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc_rem(I,j,k)
+      enddo ; enddo
+      do J=Jsq,Jeq ; do i=is,ie
+        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc_rem(i,J,k)
+      enddo ; enddo
+      if (.not.G%symmetric) &
+         call do_group_pass(CS%pass_KE_uv, G%domain)
+      do j=js,je ; do i=is,ie
+        CS%KE_visc(i,j,k) = 0.5 * G%IareaT(i,j) &
+            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
+      enddo ; enddo
+    enddo
+    if (CS%id_KE_visc_rem > 0) call post_data(CS%id_KE_visc_rem, CS%KE_visc_rem, CS%diag)
+  endif
+
+  if (associated(CS%KE_horvisc_rem)) then
+    do k=1,nz
+      do j=js,je ; do I=Isq,Ieq
+        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%diffu_visc_rem(I,j,k)
+      enddo ; enddo
+      do J=Jsq,Jeq ; do i=is,ie
+        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%diffv_visc_rem(i,J,k)
+      enddo ; enddo
+      if (.not.G%symmetric) &
+         call do_group_pass(CS%pass_KE_uv, G%domain)
+      do j=js,je ; do i=is,ie
+        CS%KE_horvisc_rem(i,j,k) = 0.5 * G%IareaT(i,j) &
+            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
+      enddo ; enddo
+    enddo
+    if (CS%id_KE_horvisc_rem > 0) call post_data(CS%id_KE_horvisc_rem, CS%KE_horvisc_rem, CS%diag)
+  endif
 
 end subroutine calculate_energy_diagnostics
 
@@ -1883,18 +1993,31 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'Potential to Kinetic Energy Conversion of Layer', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_PE_to_KE>0) call safe_alloc_ptr(CS%PE_to_KE,isd,ied,jsd,jed,nz)
+  
+ CS%id_PE_to_KE_visc_rem = register_diag_field('ocean_model', 'PE_to_KE_visc_rem', diag%axesTL, Time, &
+      'Potential to Kinetic Energy Conversion multiplied by viscous remnant fraction', &
+      'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
+  if (CS%id_PE_to_KE_visc_rem>0) call safe_alloc_ptr(CS%PE_to_KE_visc_rem,isd,ied,jsd,jed,nz)
 
   if (split) then
     CS%id_KE_BT = register_diag_field('ocean_model', 'KE_BT', diag%axesTL, Time, &
         'Barotropic contribution to Kinetic Energy', &
         'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
-    if (CS%id_KE_BT>0) call safe_alloc_ptr(CS%KE_BT,isd,ied,jsd,jed,nz)
+    if (CS%id_KE_BT_visc_rem>0) call safe_alloc_ptr(CS%KE_BT_visc_rem,isd,ied,jsd,jed,nz)
+    CS%id_KE_BT_visc_rem = register_diag_field('ocean_model', 'KE_BT_visc_rem', diag%axesTL, Time, &
+        'Barotropic contribution to Kinetic Energy multiplied by viscous remnant fraction', &
+        'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
+    if (CS%id_KE_BT_visc_rem>0) call safe_alloc_ptr(CS%KE_BT_visc_rem,isd,ied,jsd,jed,nz)
   endif
 
   CS%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, Time, &
       'Kinetic Energy Source from Coriolis and Advection', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_KE_Coradv>0) call safe_alloc_ptr(CS%KE_Coradv,isd,ied,jsd,jed,nz)
+  CS%id_KE_Coradv_visc_rem = register_diag_field('ocean_model', 'KE_Coradv_visc_rem', diag%axesTL, Time, &
+      'Kinetic Energy Source from Coriolis and Advection multiplied by viscous remnant fraction', &
+      'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
+  if (CS%id_KE_Coradv_visc_rem>0) call safe_alloc_ptr(CS%KE_Coradv_visc_rem,isd,ied,jsd,jed,nz)
 
   CS%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, Time, &
       'Kinetic Energy Source from Advection', &
@@ -1905,11 +2028,19 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'Kinetic Energy Source from Vertical Viscosity and Stresses', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_KE_visc>0) call safe_alloc_ptr(CS%KE_visc,isd,ied,jsd,jed,nz)
+  CS%id_KE_visc_rem = register_diag_field('ocean_model', 'KE_visc_rem', diag%axesTL, Time, &
+      'Kinetic Energy Source from Vertical Viscosity and Stresses multiplied by viscous remnant fraction', &
+      'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
+  if (CS%id_KE_visc_rem>0) call safe_alloc_ptr(CS%KE_visc_rem,isd,ied,jsd,jed,nz)
 
   CS%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, Time, &
       'Kinetic Energy Source from Horizontal Viscosity', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_KE_horvisc>0) call safe_alloc_ptr(CS%KE_horvisc,isd,ied,jsd,jed,nz)
+  CS%id_KE_horvisc_rem = register_diag_field('ocean_model', 'KE_horvisc_rem', diag%axesTL, Time, &
+      'Kinetic Energy Source from Horizontal Viscosity multiplied by viscous remnant fraction', &
+      'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
+  if (CS%id_KE_horvisc_rem>0) call safe_alloc_ptr(CS%KE_horvisc_rem,isd,ied,jsd,jed,nz)
 
   if (.not. adiabatic) then
     CS%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, Time, &
@@ -2307,8 +2438,11 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
   if (associated(CS%dKE_dt) .or. associated(CS%PE_to_KE) .or. &
       associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. &
       associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. &
-      associated(CS%KE_horvisc) .or. associated(CS%KE_dia)) &
+      associated(CS%KE_horvisc) .or. associated(CS%KE_dia) .or. & 
+      associated(CS%PE_to_KE_visc_rem) .or. & associated(CS%KE_BT_visc_rem) .or. &
+      associated(CS%KE_CorAdv_visc_rem) .or. associated(CS%KE_horvisc_rem)) then
     call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz)
+  endif
 
   if (associated(CS%dKE_dt)) then
     if (.not.associated(CS%du_dt)) then
@@ -2329,11 +2463,30 @@ 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)
   endif
+  if (associated(CS%KE_CorAdv_visc_rem)) then
+    call safe_alloc_ptr(ADp%CAu_visc_rem,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(ADp%CAv_visc_rem,isd,ied,JsdB,JedB,nz)
+  endif
+  if (associated(CS%PE_to_KE_visc_rem)) then
+    call safe_alloc_ptr(ADp%PFu_visc_rem,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(ADp%PFv_visc_rem,isd,ied,JsdB,JedB,nz)
+  endif
+  if (associated(CS%KE_BT_visc_rem)) then
+    call safe_alloc_ptr(ADp%u_accel_BT_visc_rem,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(ADp%v_accel_BT_visc_rem,isd,ied,JsdB,JedB,nz)
+  endif
+  if (associated(CS%KE_visc_rem)) then
+    call safe_alloc_ptr(ADp%du_dt_visc_rem,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(ADp%dv_dt_visc_rem,isd,ied,JsdB,JedB,nz)
+  endif
+  if (associated(CS%KE_horvisc_rem)) then
+    call safe_alloc_ptr(ADp%diffu_visc_rem,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(ADp%diffv_visc_rem,isd,ied,JsdB,JedB,nz)
+  endif
 
   if (associated(CS%KE_dia)) then
     call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz)
@@ -2359,11 +2512,16 @@ subroutine MOM_diagnostics_end(CS, ADp)
   if (associated(CS%KE))         deallocate(CS%KE)
   if (associated(CS%dKE_dt))     deallocate(CS%dKE_dt)
   if (associated(CS%PE_to_KE))   deallocate(CS%PE_to_KE)
+  if (associated(CS%PE_to_KE_visc_rem))   deallocate(CS%PE_to_KE_visc_rem)
   if (associated(CS%KE_BT))      deallocate(CS%KE_BT)
+  if (associated(CS%KE_BT_visc_rem))      deallocate(CS%KE_BT_visc_rem)
   if (associated(CS%KE_Coradv))  deallocate(CS%KE_Coradv)
+  if (associated(CS%KE_Coradv_visc_rem))  deallocate(CS%KE_Coradv_visc_rem)
   if (associated(CS%KE_adv))     deallocate(CS%KE_adv)
   if (associated(CS%KE_visc))    deallocate(CS%KE_visc)
+  if (associated(CS%KE_visc_rem))    deallocate(CS%KE_visc_rem)
   if (associated(CS%KE_horvisc)) deallocate(CS%KE_horvisc)
+  if (associated(CS%KE_horvisc_rem)) deallocate(CS%KE_horvisc_rem)
   if (associated(CS%KE_dia))     deallocate(CS%KE_dia)
   if (associated(CS%dv_dt))      deallocate(CS%dv_dt)
   if (associated(CS%dh_dt))      deallocate(CS%dh_dt)
@@ -2375,13 +2533,21 @@ subroutine MOM_diagnostics_end(CS, ADp)
   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_dia))  deallocate(ADp%du_dt_dia)
   if (associated(ADp%dv_dt_dia))  deallocate(ADp%dv_dt_dia)
   if (associated(ADp%du_other))   deallocate(ADp%du_other)
   if (associated(ADp%dv_other))   deallocate(ADp%dv_other)
+  if (associated(ADp%CAu_visc_rem))    deallocate(ADp%CAu_visc_rem)
+  if (associated(ADp%CAv_visc_rem))    deallocate(ADp%CAv_visc_rem)
+  if (associated(ADp%PFu_visc_rem))    deallocate(ADp%PFu_visc_rem)
+  if (associated(ADp%PFv_visc_rem))    deallocate(ADp%PFv_visc_rem)
+  if (associated(ADp%u_accel_bt_visc_rem))    deallocate(ADp%u_accel_bt_visc_rem)
+  if (associated(ADp%v_accel_bt_visc_rem))    deallocate(ADp%v_accel_bt_visc_rem)
+  if (associated(ADp%du_dt_visc_rem))    deallocate(ADp%du_dt_visc_rem)
+  if (associated(ADp%dv_dt_visc_rem))    deallocate(ADp%dv_dt_visc_rem)
 
   if (associated(ADp%diag_hfrac_u)) deallocate(ADp%diag_hfrac_u)
   if (associated(ADp%diag_hfrac_v)) deallocate(ADp%diag_hfrac_v)

From 4208525a1400641f941ac98953520cf8a4875372 Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Sat, 31 Jul 2021 12:56:31 -0600
Subject: [PATCH 11/24] Fix syntax error

---
 src/diagnostics/MOM_diagnostics.F90 | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 28b5764079..549cf5d0e2 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -2438,9 +2438,12 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
   if (associated(CS%dKE_dt) .or. associated(CS%PE_to_KE) .or. &
       associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. &
       associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. &
-      associated(CS%KE_horvisc) .or. associated(CS%KE_dia) .or. & 
-      associated(CS%PE_to_KE_visc_rem) .or. & associated(CS%KE_BT_visc_rem) .or. &
-      associated(CS%KE_CorAdv_visc_rem) .or. associated(CS%KE_horvisc_rem)) then
+      associated(CS%KE_horvisc) .or. &
+      associated(CS%KE_dia) .or. & 
+      associated(CS%PE_to_KE_visc_rem) .or. & 
+      associated(CS%KE_BT_visc_rem) .or. &
+      associated(CS%KE_CorAdv_visc_rem) .or. &
+      associated(CS%KE_horvisc_rem)) then
     call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz)
   endif
 

From f4d48393a7cd160bbf086e30f20e900af6205073 Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Sat, 31 Jul 2021 13:46:28 -0600
Subject: [PATCH 12/24] Fix some style errors

---
 src/core/MOM_barotropic.F90                   |  2 +-
 src/core/MOM_dynamics_split_RK2.F90           |  8 ++---
 src/core/MOM_variables.F90                    | 31 ++++++++++++-------
 src/diagnostics/MOM_diagnostics.F90           |  8 ++---
 .../lateral/MOM_hor_visc.F90                  |  4 +--
 .../vertical/MOM_vert_friction.F90            | 19 +++---------
 6 files changed, 35 insertions(+), 37 deletions(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 42f91b2d8e..5b6fa03bb8 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -2695,7 +2695,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
       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)
     call complete_group_pass(CS%pass_ubta_uhbta, G%Domain)
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index c55dbce684..52cd711ed9 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -356,7 +356,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [H L T-2 ~> m2 s-2].
   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].
@@ -1669,21 +1669,21 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       'Depth-integral of Barotropic Anomaly Meridional Acceleration', &
       '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', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
   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', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  
+
   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', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
   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', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  
+
   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', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 2f9c3cf643..b33eb6f92f 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -195,17 +195,26 @@ module MOM_variables
 
   real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
   real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
-  
-  real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity 
+                                                   ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity 
+                                                   ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations 
+                                                 ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations 
+                                                 ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied 
+                                                 ! by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces 
+                                                 ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied
+                                                   ! by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied 
+                                                   ! by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied 
+                                                        ! by viscous remnant fraction [L T-2 ~> m s-2]
+  real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied 
+                                                        ! by viscous remnant fraction [L T-2 ~> m s-2]
 
 end type accel_diag_ptrs
 
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 549cf5d0e2..0c5412bd72 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -109,7 +109,6 @@ module MOM_diagnostics
     KE_visc    => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3]
     KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3]
     KE_dia     => NULL(), & !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]
-  
   ! The following arrays hold diagnostics in the modified layer-integrated energy budget.
   ! Modification is through using the visc_rem_[uv]-filtered momentum equation
     PE_to_KE_visc_rem  => NULL(), &  !< potential energy to KE term [m3 s-3]
@@ -303,7 +302,6 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
     if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h = diag_pre_sync%h_state)
 
     if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state)
-    
     do k=1,nz ; do j=js,je ; do I=is-1,ie
         ADp%du_dt(I,j,k) = CS%du_dt(I,j,k)
     enddo ; enddo ; enddo
@@ -1239,7 +1237,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
     enddo
     if (CS%id_KE_dia > 0) call post_data(CS%id_KE_dia, CS%KE_dia, CS%diag)
   endif
-  
+
   if (associated(CS%PE_to_KE_visc_rem)) then
     do k=1,nz
       do j=js,je ; do I=Isq,Ieq
@@ -2439,8 +2437,8 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
       associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. &
       associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. &
       associated(CS%KE_horvisc) .or. &
-      associated(CS%KE_dia) .or. & 
-      associated(CS%PE_to_KE_visc_rem) .or. & 
+      associated(CS%KE_dia) .or. &
+      associated(CS%PE_to_KE_visc_rem) .or. &
       associated(CS%KE_BT_visc_rem) .or. &
       associated(CS%KE_CorAdv_visc_rem) .or. &
       associated(CS%KE_horvisc_rem)) then
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index bafce4b209..ad32f81f71 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -1706,7 +1706,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     call post_data(CS%id_h_diffv, h_diffv, CS%diag)
     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
@@ -2469,7 +2469,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
   if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then
     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', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index d958188a00..4549a5f190 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -217,9 +217,11 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
 
   real, allocatable, dimension(:,:,:) :: h_du_dt_visc ! h x du_dt_visc [H L T-2 ~> m2 s-2]
   real, allocatable, dimension(:,:,:) :: h_dv_dt_visc ! h x dv_dt_visc [H L T-2 ~> m2 s-2]
-  
-  real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u) [L T-2 ~> m2 s-2]
-  real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v) [L T-2 ~> m2 s-2]
+
+  real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u)
+                                                        ! [L T-2 ~> m2 s-2]
+  real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v)
+                                                        ! [L T-2 ~> m2 s-2]
 
   logical :: do_i(SZIB_(G))
   logical :: DoStokesMixing
@@ -339,11 +341,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq
       ADp%du_dt_visc(I,j,k) = (u(I,j,k) - ADp%du_dt_visc(I,j,k))*Idt
     enddo ; enddo ; endif
-    
-    !if (associated(ADp%du_dt_visc_rem)) then ; do k=1,nz ; do I=Isq,Ieq
-    !  ADp%du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + &
-    !                            (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k)
-    !enddo ; enddo ; endif
 
     if (associated(visc%taux_shelf)) then ; do I=Isq,Ieq
       visc%taux_shelf(I,j) = -GV%Rho0*CS%a1_shelf_u(I,j)*u(I,j,1) ! - u_shelf?
@@ -425,11 +422,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
       ADp%dv_dt_visc(i,J,k) = (v(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt
     enddo ; enddo ; endif
-    
-    !if (associated(ADp%dv_dt_visc_rem)) then ; do k=1,nz ; do I=Isq,Ieq
-    !    ADp%dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + &
-    !                            (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k)
-    !enddo ; enddo ; endif
 
     if (associated(visc%tauy_shelf)) then ; do i=is,ie
       visc%tauy_shelf(i,J) = -GV%Rho0*CS%a1_shelf_v(i,J)*v(i,J,1) ! - v_shelf?
@@ -1920,7 +1912,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
     call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
   endif
-
   CS%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, &
       'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)

From b1357bc2933e42fee71f64d5224fe6fb3a29d098 Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Sat, 31 Jul 2021 14:00:47 -0600
Subject: [PATCH 13/24] Fix more style errors

---
 src/core/MOM_variables.F90                     | 18 +++++++++---------
 src/diagnostics/MOM_diagnostics.F90            |  2 +-
 .../vertical/MOM_vert_friction.F90             |  4 ++--
 3 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index b33eb6f92f..8a1d0c8c60 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -195,25 +195,25 @@ module MOM_variables
 
   real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
   real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
-  real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity 
+  real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity
                                                    ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity 
+  real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity
                                                    ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations 
+  real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations
                                                  ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations 
+  real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations
                                                  ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied 
+  real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied
                                                  ! by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces 
+  real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces
                                                  ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
   real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied
                                                    ! by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied 
+  real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied
                                                    ! by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied 
+  real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied
                                                         ! by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied 
+  real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied
                                                         ! by viscous remnant fraction [L T-2 ~> m s-2]
 
 end type accel_diag_ptrs
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 0c5412bd72..a81fbca4da 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -1991,7 +1991,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'Potential to Kinetic Energy Conversion of Layer', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_PE_to_KE>0) call safe_alloc_ptr(CS%PE_to_KE,isd,ied,jsd,jed,nz)
-  
+
  CS%id_PE_to_KE_visc_rem = register_diag_field('ocean_model', 'PE_to_KE_visc_rem', diag%axesTL, Time, &
       'Potential to Kinetic Energy Conversion multiplied by viscous remnant fraction', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 4549a5f190..83e5a391fa 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -527,7 +527,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     call post_data(CS%id_h_dv_dt_visc, h_dv_dt_visc, CS%diag)
     deallocate(h_dv_dt_visc)
   endif
-  
+
   if (CS%id_du_dt_visc_rem > 0) then
     allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
     du_dt_visc_rem(:,:,:) = 0.0
@@ -1903,7 +1903,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
     call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
-  
+
   CS%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, &
       'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)

From be09de03ebd114b43bfca5f2c1ae2687f849366c Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Tue, 3 Aug 2021 18:58:58 -0600
Subject: [PATCH 14/24] Remove added KE budget diagnostics

---
 src/core/MOM_variables.F90          |  20 ----
 src/diagnostics/MOM_diagnostics.F90 | 175 +---------------------------
 2 files changed, 4 insertions(+), 191 deletions(-)

diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 8a1d0c8c60..6f6cf12214 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -195,26 +195,6 @@ module MOM_variables
 
   real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
   real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
-  real, pointer :: diffu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to along isopycnal viscosity
-                                                   ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: diffv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to along isopycnal viscosity
-                                                   ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: CAu_visc_rem(:,:,:) => NULL() !< Zonal Coriolis and momentum advection accelerations
-                                                 ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: CAv_visc_rem(:,:,:) => NULL() !< Meridional Coriolis and momentum advection accelerations
-                                                 ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: PFu_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to pressure forces multiplied
-                                                 ! by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: PFv_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to pressure forces
-                                                 ! multiplied by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: du_dt_visc_rem(:,:,:) => NULL() !< Zonal acceleration due to vertical viscosity multiplied
-                                                   ! by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: dv_dt_visc_rem(:,:,:) => NULL() !< Meridional acceleration due to vertical viscosity multiplied
-                                                   ! by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: u_accel_bt_visc_rem(:,:,:) => NULL() !< Zonal barotropic-solver acceleration multiplied
-                                                        ! by viscous remnant fraction [L T-2 ~> m s-2]
-  real, pointer :: v_accel_bt_visc_rem(:,:,:) => NULL() !< Meridional barotropic-solver acceleration multiplied
-                                                        ! by viscous remnant fraction [L T-2 ~> m s-2]
 
 end type accel_diag_ptrs
 
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index a81fbca4da..88a805fd04 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -108,15 +108,7 @@ module MOM_diagnostics
     KE_adv     => NULL(), & !< KE source from along-layer advection [H L2 T-3 ~> m3 s-3]
     KE_visc    => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3]
     KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3]
-    KE_dia     => NULL(), & !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]
-  ! The following arrays hold diagnostics in the modified layer-integrated energy budget.
-  ! Modification is through using the visc_rem_[uv]-filtered momentum equation
-    PE_to_KE_visc_rem  => NULL(), &  !< potential energy to KE term [m3 s-3]
-    KE_BT_visc_rem     => NULL(), &  !< barotropic contribution to KE term [m3 s-3]
-    KE_CorAdv_visc_rem => NULL(), &  !< KE source from the combined Coriolis and
-                            !! advection terms [H L2 T-3 ~> m3 s-3].
-    KE_visc_rem    => NULL(), &  !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3]
-    KE_horvisc_rem => NULL() !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3]
+    KE_dia     => NULL(),   !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]
 
   !>@{ Diagnostic IDs
   integer :: id_u   = -1,   id_v   = -1, id_h = -1
@@ -132,10 +124,6 @@ module MOM_diagnostics
   integer :: id_KE_Coradv      = -1
   integer :: id_KE_adv         = -1, id_KE_visc        = -1
   integer :: id_KE_horvisc     = -1, id_KE_dia         = -1
-  integer :: id_PE_to_KE_visc_rem = -1, id_KE_BT_visc_rem = -1
-  integer :: id_KE_Coradv_visc_rem = -1
-  integer :: id_KE_visc_rem = -1
-  integer :: id_KE_horvisc_rem = -1
   integer :: id_uh_Rlay        = -1, id_vh_Rlay        = -1
   integer :: id_uhGM_Rlay      = -1, id_vhGM_Rlay      = -1
   integer :: id_h_Rlay         = -1, id_Rd1            = -1
@@ -1066,10 +1054,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
   if (.not.G%symmetric) then
     if (associated(CS%dKE_dt) .OR. associated(CS%PE_to_KE) .OR. associated(CS%KE_BT) .OR. &
         associated(CS%KE_CorAdv) .OR. associated(CS%KE_adv) .OR. associated(CS%KE_visc) .OR. &
-        associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) .OR. &
-        associated(CS%PE_to_KE_visc_rem) .OR. associated(CS%KE_BT_visc_rem) .OR. &
-        associated(CS%KE_CorAdv_visc_rem) .OR. associated(CS%KE_visc_rem) .OR. &
-        associated(CS%KE_horvisc_rem)) then
+        associated(CS%KE_horvisc) .OR. associated(CS%KE_dia)) then
       call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East)
     endif
   endif
@@ -1238,100 +1223,6 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
     if (CS%id_KE_dia > 0) call post_data(CS%id_KE_dia, CS%KE_dia, CS%diag)
   endif
 
-  if (associated(CS%PE_to_KE_visc_rem)) then
-    do k=1,nz
-      do j=js,je ; do I=Isq,Ieq
-        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%PFu_visc_rem(I,j,k)
-      enddo ; enddo
-      do J=Jsq,Jeq ; do i=is,ie
-        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%PFv_visc_rem(i,J,k)
-      enddo ; enddo
-      if (.not.G%symmetric) &
-         call do_group_pass(CS%pass_KE_uv, G%domain)
-      do j=js,je ; do i=is,ie
-        CS%PE_to_KE_visc_rem(i,j,k) = 0.5 * G%IareaT(i,j) &
-            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
-      enddo ; enddo
-    enddo
-    if (CS%id_PE_to_KE_visc_rem > 0) call post_data(CS%id_PE_to_KE_visc_rem, CS%PE_to_KE_visc_rem, CS%diag)
-  endif
-
-  if (associated(CS%KE_BT_visc_rem)) then
-    do k=1,nz
-      do j=js,je ; do I=Isq,Ieq
-        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%u_accel_bt_visc_rem(I,j,k)
-      enddo ; enddo
-      do J=Jsq,Jeq ; do i=is,ie
-        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%v_accel_bt_visc_rem(i,J,k)
-      enddo ; enddo
-      if (.not.G%symmetric) &
-         call do_group_pass(CS%pass_KE_uv, G%domain)
-      do j=js,je ; do i=is,ie
-        CS%KE_BT_visc_rem(i,j,k) = 0.5 * G%IareaT(i,j) &
-            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
-      enddo ; enddo
-    enddo
-    if (CS%id_KE_BT_visc_rem > 0) call post_data(CS%id_KE_BT_visc_rem, CS%KE_BT_visc_rem, CS%diag)
-  endif
-
-  if (associated(CS%KE_CorAdv_visc_rem)) then
-    do k=1,nz
-      do j=js,je ; do I=Isq,Ieq
-        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%CAu_visc_rem(I,j,k)
-      enddo ; enddo
-      do J=Jsq,Jeq ; do i=is,ie
-        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%CAv_visc_rem(i,J,k)
-      enddo ; enddo
-      do j=js,je ; do i=is,ie
-        KE_h(i,j) = -CS%KE(i,j,k) * G%IareaT(i,j) &
-            * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k))
-      enddo ; enddo
-      if (.not.G%symmetric) &
-         call do_group_pass(CS%pass_KE_uv, G%domain)
-      do j=js,je ; do i=is,ie
-        CS%KE_CorAdv_visc_rem(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) &
-            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
-      enddo ; enddo
-    enddo
-    if (CS%id_KE_CorAdv_visc_rem > 0) call post_data(CS%id_KE_Coradv_visc_rem, CS%KE_Coradv_visc_rem, CS%diag)
-  endif
-
-  if (associated(CS%KE_visc_rem)) then
-    do k=1,nz
-      do j=js,je ; do I=Isq,Ieq
-        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc_rem(I,j,k)
-      enddo ; enddo
-      do J=Jsq,Jeq ; do i=is,ie
-        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc_rem(i,J,k)
-      enddo ; enddo
-      if (.not.G%symmetric) &
-         call do_group_pass(CS%pass_KE_uv, G%domain)
-      do j=js,je ; do i=is,ie
-        CS%KE_visc(i,j,k) = 0.5 * G%IareaT(i,j) &
-            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
-      enddo ; enddo
-    enddo
-    if (CS%id_KE_visc_rem > 0) call post_data(CS%id_KE_visc_rem, CS%KE_visc_rem, CS%diag)
-  endif
-
-  if (associated(CS%KE_horvisc_rem)) then
-    do k=1,nz
-      do j=js,je ; do I=Isq,Ieq
-        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%diffu_visc_rem(I,j,k)
-      enddo ; enddo
-      do J=Jsq,Jeq ; do i=is,ie
-        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%diffv_visc_rem(i,J,k)
-      enddo ; enddo
-      if (.not.G%symmetric) &
-         call do_group_pass(CS%pass_KE_uv, G%domain)
-      do j=js,je ; do i=is,ie
-        CS%KE_horvisc_rem(i,j,k) = 0.5 * G%IareaT(i,j) &
-            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
-      enddo ; enddo
-    enddo
-    if (CS%id_KE_horvisc_rem > 0) call post_data(CS%id_KE_horvisc_rem, CS%KE_horvisc_rem, CS%diag)
-  endif
-
 end subroutine calculate_energy_diagnostics
 
 !> This subroutine registers fields to calculate a diagnostic time derivative.
@@ -1992,30 +1883,17 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_PE_to_KE>0) call safe_alloc_ptr(CS%PE_to_KE,isd,ied,jsd,jed,nz)
 
- CS%id_PE_to_KE_visc_rem = register_diag_field('ocean_model', 'PE_to_KE_visc_rem', diag%axesTL, Time, &
-      'Potential to Kinetic Energy Conversion multiplied by viscous remnant fraction', &
-      'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
-  if (CS%id_PE_to_KE_visc_rem>0) call safe_alloc_ptr(CS%PE_to_KE_visc_rem,isd,ied,jsd,jed,nz)
-
   if (split) then
     CS%id_KE_BT = register_diag_field('ocean_model', 'KE_BT', diag%axesTL, Time, &
         'Barotropic contribution to Kinetic Energy', &
         'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
-    if (CS%id_KE_BT_visc_rem>0) call safe_alloc_ptr(CS%KE_BT_visc_rem,isd,ied,jsd,jed,nz)
-    CS%id_KE_BT_visc_rem = register_diag_field('ocean_model', 'KE_BT_visc_rem', diag%axesTL, Time, &
-        'Barotropic contribution to Kinetic Energy multiplied by viscous remnant fraction', &
-        'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
-    if (CS%id_KE_BT_visc_rem>0) call safe_alloc_ptr(CS%KE_BT_visc_rem,isd,ied,jsd,jed,nz)
+    if (CS%id_KE_BT>0) call safe_alloc_ptr(CS%KE_BT,isd,ied,jsd,jed,nz)
   endif
 
   CS%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, Time, &
       'Kinetic Energy Source from Coriolis and Advection', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_KE_Coradv>0) call safe_alloc_ptr(CS%KE_Coradv,isd,ied,jsd,jed,nz)
-  CS%id_KE_Coradv_visc_rem = register_diag_field('ocean_model', 'KE_Coradv_visc_rem', diag%axesTL, Time, &
-      'Kinetic Energy Source from Coriolis and Advection multiplied by viscous remnant fraction', &
-      'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
-  if (CS%id_KE_Coradv_visc_rem>0) call safe_alloc_ptr(CS%KE_Coradv_visc_rem,isd,ied,jsd,jed,nz)
 
   CS%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, Time, &
       'Kinetic Energy Source from Advection', &
@@ -2026,19 +1904,11 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'Kinetic Energy Source from Vertical Viscosity and Stresses', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_KE_visc>0) call safe_alloc_ptr(CS%KE_visc,isd,ied,jsd,jed,nz)
-  CS%id_KE_visc_rem = register_diag_field('ocean_model', 'KE_visc_rem', diag%axesTL, Time, &
-      'Kinetic Energy Source from Vertical Viscosity and Stresses multiplied by viscous remnant fraction', &
-      'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
-  if (CS%id_KE_visc_rem>0) call safe_alloc_ptr(CS%KE_visc_rem,isd,ied,jsd,jed,nz)
 
   CS%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, Time, &
       'Kinetic Energy Source from Horizontal Viscosity', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_KE_horvisc>0) call safe_alloc_ptr(CS%KE_horvisc,isd,ied,jsd,jed,nz)
-  CS%id_KE_horvisc_rem = register_diag_field('ocean_model', 'KE_horvisc_rem', diag%axesTL, Time, &
-      'Kinetic Energy Source from Horizontal Viscosity multiplied by viscous remnant fraction', &
-      'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
-  if (CS%id_KE_horvisc_rem>0) call safe_alloc_ptr(CS%KE_horvisc_rem,isd,ied,jsd,jed,nz)
 
   if (.not. adiabatic) then
     CS%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, Time, &
@@ -2437,11 +2307,7 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
       associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. &
       associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. &
       associated(CS%KE_horvisc) .or. &
-      associated(CS%KE_dia) .or. &
-      associated(CS%PE_to_KE_visc_rem) .or. &
-      associated(CS%KE_BT_visc_rem) .or. &
-      associated(CS%KE_CorAdv_visc_rem) .or. &
-      associated(CS%KE_horvisc_rem)) then
+      associated(CS%KE_dia)) then
     call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz)
   endif
 
@@ -2468,26 +2334,6 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
     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)
   endif
-  if (associated(CS%KE_CorAdv_visc_rem)) then
-    call safe_alloc_ptr(ADp%CAu_visc_rem,IsdB,IedB,jsd,jed,nz)
-    call safe_alloc_ptr(ADp%CAv_visc_rem,isd,ied,JsdB,JedB,nz)
-  endif
-  if (associated(CS%PE_to_KE_visc_rem)) then
-    call safe_alloc_ptr(ADp%PFu_visc_rem,IsdB,IedB,jsd,jed,nz)
-    call safe_alloc_ptr(ADp%PFv_visc_rem,isd,ied,JsdB,JedB,nz)
-  endif
-  if (associated(CS%KE_BT_visc_rem)) then
-    call safe_alloc_ptr(ADp%u_accel_BT_visc_rem,IsdB,IedB,jsd,jed,nz)
-    call safe_alloc_ptr(ADp%v_accel_BT_visc_rem,isd,ied,JsdB,JedB,nz)
-  endif
-  if (associated(CS%KE_visc_rem)) then
-    call safe_alloc_ptr(ADp%du_dt_visc_rem,IsdB,IedB,jsd,jed,nz)
-    call safe_alloc_ptr(ADp%dv_dt_visc_rem,isd,ied,JsdB,JedB,nz)
-  endif
-  if (associated(CS%KE_horvisc_rem)) then
-    call safe_alloc_ptr(ADp%diffu_visc_rem,IsdB,IedB,jsd,jed,nz)
-    call safe_alloc_ptr(ADp%diffv_visc_rem,isd,ied,JsdB,JedB,nz)
-  endif
 
   if (associated(CS%KE_dia)) then
     call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz)
@@ -2513,16 +2359,11 @@ subroutine MOM_diagnostics_end(CS, ADp)
   if (associated(CS%KE))         deallocate(CS%KE)
   if (associated(CS%dKE_dt))     deallocate(CS%dKE_dt)
   if (associated(CS%PE_to_KE))   deallocate(CS%PE_to_KE)
-  if (associated(CS%PE_to_KE_visc_rem))   deallocate(CS%PE_to_KE_visc_rem)
   if (associated(CS%KE_BT))      deallocate(CS%KE_BT)
-  if (associated(CS%KE_BT_visc_rem))      deallocate(CS%KE_BT_visc_rem)
   if (associated(CS%KE_Coradv))  deallocate(CS%KE_Coradv)
-  if (associated(CS%KE_Coradv_visc_rem))  deallocate(CS%KE_Coradv_visc_rem)
   if (associated(CS%KE_adv))     deallocate(CS%KE_adv)
   if (associated(CS%KE_visc))    deallocate(CS%KE_visc)
-  if (associated(CS%KE_visc_rem))    deallocate(CS%KE_visc_rem)
   if (associated(CS%KE_horvisc)) deallocate(CS%KE_horvisc)
-  if (associated(CS%KE_horvisc_rem)) deallocate(CS%KE_horvisc_rem)
   if (associated(CS%KE_dia))     deallocate(CS%KE_dia)
   if (associated(CS%dv_dt))      deallocate(CS%dv_dt)
   if (associated(CS%dh_dt))      deallocate(CS%dh_dt)
@@ -2541,14 +2382,6 @@ subroutine MOM_diagnostics_end(CS, ADp)
   if (associated(ADp%dv_dt_dia))  deallocate(ADp%dv_dt_dia)
   if (associated(ADp%du_other))   deallocate(ADp%du_other)
   if (associated(ADp%dv_other))   deallocate(ADp%dv_other)
-  if (associated(ADp%CAu_visc_rem))    deallocate(ADp%CAu_visc_rem)
-  if (associated(ADp%CAv_visc_rem))    deallocate(ADp%CAv_visc_rem)
-  if (associated(ADp%PFu_visc_rem))    deallocate(ADp%PFu_visc_rem)
-  if (associated(ADp%PFv_visc_rem))    deallocate(ADp%PFv_visc_rem)
-  if (associated(ADp%u_accel_bt_visc_rem))    deallocate(ADp%u_accel_bt_visc_rem)
-  if (associated(ADp%v_accel_bt_visc_rem))    deallocate(ADp%v_accel_bt_visc_rem)
-  if (associated(ADp%du_dt_visc_rem))    deallocate(ADp%du_dt_visc_rem)
-  if (associated(ADp%dv_dt_visc_rem))    deallocate(ADp%dv_dt_visc_rem)
 
   if (associated(ADp%diag_hfrac_u)) deallocate(ADp%diag_hfrac_u)
   if (associated(ADp%diag_hfrac_v)) deallocate(ADp%diag_hfrac_v)

From fbdac20c99dfadd1f1dbba8824755c37ec8372a0 Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Tue, 3 Aug 2021 20:30:10 -0600
Subject: [PATCH 15/24] Fix small bug

---
 src/diagnostics/MOM_diagnostics.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 88a805fd04..af2b357c3c 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -108,7 +108,7 @@ module MOM_diagnostics
     KE_adv     => NULL(), & !< KE source from along-layer advection [H L2 T-3 ~> m3 s-3]
     KE_visc    => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3]
     KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3]
-    KE_dia     => NULL(),   !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]
+    KE_dia     => NULL()    !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]
 
   !>@{ Diagnostic IDs
   integer :: id_u   = -1,   id_v   = -1, id_h = -1

From 86741cd1d3b049303db325d61089263eda90979b Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Tue, 3 Aug 2021 21:51:46 -0600
Subject: [PATCH 16/24] Small style fix

---
 src/diagnostics/MOM_diagnostics.F90 | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index af2b357c3c..2c567ccceb 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -1054,7 +1054,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
   if (.not.G%symmetric) then
     if (associated(CS%dKE_dt) .OR. associated(CS%PE_to_KE) .OR. associated(CS%KE_BT) .OR. &
         associated(CS%KE_CorAdv) .OR. associated(CS%KE_adv) .OR. associated(CS%KE_visc) .OR. &
-        associated(CS%KE_horvisc) .OR. associated(CS%KE_dia)) then
+        associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) ) then
       call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East)
     endif
   endif
@@ -2306,10 +2306,8 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
   if (associated(CS%dKE_dt) .or. associated(CS%PE_to_KE) .or. &
       associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. &
       associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. &
-      associated(CS%KE_horvisc) .or. &
-      associated(CS%KE_dia)) then
+      associated(CS%KE_horvisc) .or. associated(CS%KE_dia)) &
     call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz)
-  endif
 
   if (associated(CS%dKE_dt)) then
     if (.not.associated(CS%du_dt)) then

From a40a90f3f0a18b9f2b5bd71672ef835a8f478de1 Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Wed, 4 Aug 2021 15:43:43 -0600
Subject: [PATCH 17/24] Fix allocation problems

---
 src/core/MOM_dynamics_split_RK2.F90            | 16 +++++++++++-----
 src/core/MOM_variables.F90                     |  2 ++
 src/diagnostics/MOM_diagnostics.F90            |  2 +-
 .../vertical/MOM_vert_friction.F90             | 18 +++++++++---------
 4 files changed, 23 insertions(+), 15 deletions(-)

diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 52cd711ed9..924e9614d4 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -1115,7 +1115,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     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%visc_rem_u(I,j,k)
+      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)
@@ -1124,7 +1124,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     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%visc_rem_v(i,J,k)
+      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)
@@ -1133,7 +1133,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     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%visc_rem_u(I,j,k)
+      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)
@@ -1142,7 +1142,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     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%visc_rem_v(i,J,k)
+      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)
@@ -1151,7 +1151,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     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%visc_rem_u(I,j,k)
+      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)
@@ -1673,23 +1673,29 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   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', 'm2 s-2', &
       conversion=GV%H_to_m*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', 'm2 s-2', &
       conversion=GV%H_to_m*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', 'm2 s-2', &
       conversion=GV%H_to_m*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', 'm2 s-2', &
       conversion=GV%H_to_m*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', 'm2 s-2', &
       conversion=GV%H_to_m*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', 'm2 s-2', &
       conversion=GV%H_to_m*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)
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 6f6cf12214..387ef98492 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -195,6 +195,8 @@ module MOM_variables
 
   real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
   real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
+ 
+  integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1
 
 end type accel_diag_ptrs
 
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 2c567ccceb..e715edcaf0 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -215,7 +215,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
                                                  !! [H L2 T-1 ~> m3 s-1 or kg s-1].
   type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various
                                                  !! thermodynamic variables.
-  type(accel_diag_ptrs),   intent(in)    :: ADp  !< structure with pointers to
+  type(accel_diag_ptrs),   intent(inout)    :: ADp  !< structure with pointers to
                                                  !! accelerations in momentum equation.
   type(cont_diag_ptrs),    intent(in)    :: CDp  !< structure with pointers to
                                                  !! terms in continuity equation.
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 83e5a391fa..0182e14f90 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -127,7 +127,7 @@ module MOM_vert_friction
   ! integer :: id_hf_du_dt_visc    = -1, id_hf_dv_dt_visc    = -1
   integer :: id_h_du_dt_visc    = -1, id_h_dv_dt_visc    = -1
   integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1
-  integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1
+  ! integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1  ! moved to MOM_variables
   !>@}
 
   type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure
@@ -528,24 +528,24 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     deallocate(h_dv_dt_visc)
   endif
 
-  if (CS%id_du_dt_visc_rem > 0) then
+  if (ADp%id_du_dt_visc_rem > 0) then
     allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
     du_dt_visc_rem(:,:,:) = 0.0
     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
       du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + &
                                 (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k)
     enddo ; enddo ; enddo
-    call post_data(CS%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag)
+    call post_data(ADp%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag)
     deallocate(du_dt_visc_rem)
   endif
-  if (CS%id_dv_dt_visc_rem > 0) then
+  if (ADp%id_dv_dt_visc_rem > 0) then
     allocate(dv_dt_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
     dv_dt_visc_rem(:,:,:) = 0.0
     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
       dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + &
                                 (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k)
     enddo ; enddo ; enddo
-    call post_data(CS%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag)
+    call post_data(ADp%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag)
     deallocate(dv_dt_visc_rem)
   endif
 
@@ -1904,18 +1904,18 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
     call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
-  CS%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, &
+  ADp%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, &
       'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if (CS%id_du_dt_visc_rem > 0) then
+  if (ADp%id_du_dt_visc_rem > 0) then
     call safe_alloc_ptr(ADp%du_dt,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
   endif
-  CS%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, &
+  ADp%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, &
       'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if (CS%id_dv_dt_visc_rem > 0) then
+  if (ADp%id_dv_dt_visc_rem > 0) then
     call safe_alloc_ptr(ADp%dv_dt,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)

From 156d02c122207cb273d286db3dfdac8242cd5fea Mon Sep 17 00:00:00 2001
From: Gustavo Marques <gmarques@ucar.edu>
Date: Wed, 4 Aug 2021 16:26:24 -0600
Subject: [PATCH 18/24] Make use_drag_rate independent of MEKE_damping

Remove MEKE_damping from use_drag_rate so we can use linear
dissipation without bottom drag. This is needed for the
GEOMETRIC scheme.
---
 src/parameterizations/lateral/MOM_MEKE.F90 | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
index 850f94cff2..50fd35bae9 100644
--- a/src/parameterizations/lateral/MOM_MEKE.F90
+++ b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -175,8 +175,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.

From f31e5e64b3992bea0ff6f4b54abfccb7226f92cc Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Wed, 4 Aug 2021 16:39:51 -0600
Subject: [PATCH 19/24] Remove du_dt_visc_rem diagnostic

---
 src/core/MOM_variables.F90                    |  6 +--
 src/diagnostics/MOM_diagnostics.F90           |  6 ---
 .../vertical/MOM_vert_friction.F90            | 44 -------------------
 3 files changed, 1 insertion(+), 55 deletions(-)

diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index 387ef98492..f7f35ed2d1 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -161,8 +161,6 @@ module MOM_variables
 
   ! Each of the following fields has nz layers.
   real, pointer, dimension(:,:,:) :: &
-    du_dt => NULL(), &     !< Zonal acceleration [L T-2 ~> m s-2]
-    dv_dt => NULL(), &     !< Meridional acceleration [L T-2 ~> m s-2]
     diffu => NULL(), &     !< Zonal acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
     diffv => NULL(), &     !< Meridional acceleration due to along isopycnal viscosity [L T-2 ~> m s-2]
     CAu => NULL(), &       !< Zonal Coriolis and momentum advection accelerations [L T-2 ~> m s-2]
@@ -174,7 +172,7 @@ module MOM_variables
     du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
     dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
     u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2]
-    v_accel_bt => NULL()  !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
+    v_accel_bt => NULL()   !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
   real, pointer, dimension(:,:,:) :: du_other => NULL()
                            !< Zonal velocity changes due to any other processes that are
                            !! not due to any explicit accelerations [L T-1 ~> m s-1].
@@ -195,8 +193,6 @@ module MOM_variables
 
   real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
   real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
- 
-  integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1
 
 end type accel_diag_ptrs
 
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index e715edcaf0..cb0f10e6ff 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -290,12 +290,6 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
     if (CS%id_dv_dt>0) call post_data(CS%id_dv_dt, CS%dv_dt, CS%diag, alt_h = diag_pre_sync%h_state)
 
     if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state)
-    do k=1,nz ; do j=js,je ; do I=is-1,ie
-        ADp%du_dt(I,j,k) = CS%du_dt(I,j,k)
-    enddo ; enddo ; enddo
-    do k=1,nz ; do J=js-1,je ; do i=is,ie
-        ADp%dv_dt(i,J,k) = CS%dv_dt(i,J,k)
-    enddo ; enddo ; enddo
 
     !! Diagnostics for terms multiplied by fractional thicknesses
 
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 0182e14f90..081179eb41 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -127,7 +127,6 @@ module MOM_vert_friction
   ! integer :: id_hf_du_dt_visc    = -1, id_hf_dv_dt_visc    = -1
   integer :: id_h_du_dt_visc    = -1, id_h_dv_dt_visc    = -1
   integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1
-  ! integer :: id_du_dt_visc_rem = -1, id_dv_dt_visc_rem = -1  ! moved to MOM_variables
   !>@}
 
   type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure
@@ -218,11 +217,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
   real, allocatable, dimension(:,:,:) :: h_du_dt_visc ! h x du_dt_visc [H L T-2 ~> m2 s-2]
   real, allocatable, dimension(:,:,:) :: h_dv_dt_visc ! h x dv_dt_visc [H L T-2 ~> m2 s-2]
 
-  real, allocatable, dimension(:,:,:) :: du_dt_visc_rem ! du_dt_visc x visc_rem_u + du_dt x (1-visc_rem_u)
-                                                        ! [L T-2 ~> m2 s-2]
-  real, allocatable, dimension(:,:,:) :: dv_dt_visc_rem ! dv_dt_visc x visc_rem_v + dv_dt x (1-visc_rem_v)
-                                                        ! [L T-2 ~> m2 s-2]
-
   logical :: do_i(SZIB_(G))
   logical :: DoStokesMixing
 
@@ -528,27 +522,6 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     deallocate(h_dv_dt_visc)
   endif
 
-  if (ADp%id_du_dt_visc_rem > 0) then
-    allocate(du_dt_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
-    du_dt_visc_rem(:,:,:) = 0.0
-    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
-      du_dt_visc_rem(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%visc_rem_u(I,j,k) + &
-                                (1-ADp%visc_rem_u(I,j,k)) * ADp%du_dt(I,j,k)
-    enddo ; enddo ; enddo
-    call post_data(ADp%id_du_dt_visc_rem, du_dt_visc_rem, CS%diag)
-    deallocate(du_dt_visc_rem)
-  endif
-  if (ADp%id_dv_dt_visc_rem > 0) then
-    allocate(dv_dt_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
-    dv_dt_visc_rem(:,:,:) = 0.0
-    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
-      dv_dt_visc_rem(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%visc_rem_v(i,J,k) + &
-                                (1-ADp%visc_rem_v(i,J,k)) * ADp%dv_dt(i,J,k)
-    enddo ; enddo ; enddo
-    call post_data(ADp%id_dv_dt_visc_rem, dv_dt_visc_rem, CS%diag)
-    deallocate(dv_dt_visc_rem)
-  endif
-
 end subroutine vertvisc
 
 !> Calculate the fraction of momentum originally in a layer that remains
@@ -1904,23 +1877,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
     call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
-  ADp%id_du_dt_visc_rem = register_diag_field('ocean_model', 'du_dt_visc_rem', diag%axesCuL, Time, &
-      'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if (ADp%id_du_dt_visc_rem > 0) then
-    call safe_alloc_ptr(ADp%du_dt,IsdB,IedB,jsd,jed,nz)
-    call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
-    call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
-  endif
-  ADp%id_dv_dt_visc_rem = register_diag_field('ocean_model', 'dv_dt_visc_rem', diag%axesCvL, Time, &
-      'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant fraction', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if (ADp%id_dv_dt_visc_rem > 0) then
-    call safe_alloc_ptr(ADp%dv_dt,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
-  endif
-
   if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &
     call PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS%PointAccel_CSp)
 

From ca0fb8950851b235a18e83ae0a254b6670abb3bc Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Wed, 4 Aug 2021 17:38:57 -0600
Subject: [PATCH 20/24] Fix small inconsistency: CS%visc_rem_v -->
 CS%ADp%visc_rem_v

---
 src/core/MOM_dynamics_split_RK2.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 924e9614d4..ebe53fc908 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -1160,7 +1160,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     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%visc_rem_v(i,J,k)
+      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)

From a874d6df79e79b84a7706dd7b77a10a97f4c35ea Mon Sep 17 00:00:00 2001
From: NoraLoose <nora.loose@gmail.com>
Date: Wed, 4 Aug 2021 17:40:44 -0600
Subject: [PATCH 21/24] Don't need inout for ADp anymore

---
 src/diagnostics/MOM_diagnostics.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index cb0f10e6ff..f874d08a12 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -215,7 +215,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
                                                  !! [H L2 T-1 ~> m3 s-1 or kg s-1].
   type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various
                                                  !! thermodynamic variables.
-  type(accel_diag_ptrs),   intent(inout)    :: ADp  !< structure with pointers to
+  type(accel_diag_ptrs),   intent(in)    :: ADp  !< structure with pointers to
                                                  !! accelerations in momentum equation.
   type(cont_diag_ptrs),    intent(in)    :: CDp  !< structure with pointers to
                                                  !! terms in continuity equation.

From 8ee65a6ff45efb6e6c5deb780e0b12f87361e6a9 Mon Sep 17 00:00:00 2001
From: Gustavo Marques <gmarques@ucar.edu>
Date: Mon, 23 Aug 2021 11:53:49 -0600
Subject: [PATCH 22/24] Fix units of *_visc_rem terms

---
 src/core/MOM_dynamics_split_RK2.F90           | 24 +++++++++----------
 .../lateral/MOM_hor_visc.F90                  | 12 +++++-----
 2 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 96a0a5f92f..a168fe1319 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -1677,30 +1677,30 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   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', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      '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', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      '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', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      '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', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      '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', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      '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', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      '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)
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index 592fca4d24..61127c745f 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -281,8 +281,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 ~> m2 s-2]
-  real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [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]
@@ -2471,15 +2471,15 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
   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', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      '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', 'm2 s-2', &
-      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+      '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

From 868c11a2a91603b6673cf31596457507c19f14e7 Mon Sep 17 00:00:00 2001
From: Gustavo Marques <gmarques@ucar.edu>
Date: Mon, 30 Aug 2021 10:40:58 -0600
Subject: [PATCH 23/24] Always initialize drag_rate

This patch now always initializes drag_rate, either to zero or to
the bottom drag formula.
It also changes the logic to always execute the second Strang
splitting calculation.

Removed unused code, including drag_rate_J15 and some
commented code.
---
 src/parameterizations/lateral/MOM_MEKE.F90 | 41 ++++++++--------------
 1 file changed, 14 insertions(+), 27 deletions(-)

diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
index 50fd35bae9..7b92704651 100644
--- a/src/parameterizations/lateral/MOM_MEKE.F90
+++ b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -130,8 +130,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].
@@ -230,12 +228,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
@@ -339,12 +331,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       enddo ; enddo
     endif
 
-    ! Increase EKE by a full time-steps worth of source
-    !$OMP parallel do default(shared)
-    do j=js,je ; do i=is,ie
-      MEKE%MEKE(i,j) = (MEKE%MEKE(i,j) + sdt*src(i,j))*G%mask2dT(i,j)
-    enddo ; enddo
-
     if (use_drag_rate) then
       ! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)
       !$OMP parallel do default(shared)
@@ -352,6 +338,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
@@ -520,23 +511,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)

From ba6f56d5f002ebf2093ecb8e11b615acd65d827b Mon Sep 17 00:00:00 2001
From: Gustavo Marques <gmarques@ucar.edu>
Date: Mon, 30 Aug 2021 13:33:50 -0600
Subject: [PATCH 24/24] Add code that was deleted unintentionally

---
 src/parameterizations/lateral/MOM_MEKE.F90 | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
index 66f75c82e3..6ce0e58471 100644
--- a/src/parameterizations/lateral/MOM_MEKE.F90
+++ b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -331,6 +331,12 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       enddo ; enddo
     endif
 
+    ! Increase EKE by a full time-steps worth of source
+    !$OMP parallel do default(shared)
+    do j=js,je ; do i=is,ie
+      MEKE%MEKE(i,j) = (MEKE%MEKE(i,j) + sdt*src(i,j))*G%mask2dT(i,j)
+    enddo ; enddo
+
     if (use_drag_rate) then
       ! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)
       !$OMP parallel do default(shared)