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/20] 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/20] 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/20] 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/20] 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/20] 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/20] 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/20] 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 48239bb37844a51b7539178618247502b42dc9e5 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 11:12:20 -0500
Subject: [PATCH 08/20] MOM_variables merge

---
 src/core/MOM_variables.F90 | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index d81cf28e17..886ee77510 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -188,6 +188,9 @@ module MOM_variables
 
   real, pointer :: diag_hfrac_u(:,:,:) => NULL() !< Fractional layer thickness at u points
   real, pointer :: diag_hfrac_v(:,:,:) => NULL() !< Fractional layer thickness at v points
+  real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points
+  real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points
+
 end type accel_diag_ptrs
 
 !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances.

From 373da29dc45c49684439b785acfeb9cff081183e Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 11:19:31 -0500
Subject: [PATCH 09/20] Return layer thicknesses from barotropic module

---
 src/core/MOM_barotropic.F90 | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 628f8e1c39..1b609d0981 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -2632,6 +2632,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     enddo ; enddo ; enddo
   endif
 
+  if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_hu))) then
+    do k=1,nz ; do j=js,je ; do I=is-1,ie
+      ADp%diag_hu(I,j,k) = BT_cont%h_u(I,j,k)
+    enddo ; enddo ; enddo
+  endif
+  if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_hv))) then
+    do k=1,nz ; do J=js-1,je ; do i=is,ie
+      ADp%diag_hv(i,J,k) = BT_cont%h_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)

From a6055e1652b44df8ec7325ed5b2c43dcd690ad44 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 11:56:22 -0500
Subject: [PATCH 10/20] Some thickness multiplied 3D diagnostics

---
 src/core/MOM_CoriolisAdv.F90                  | 70 +++++++++++++++++++
 .../lateral/MOM_hor_visc.F90                  | 32 +++++++++
 2 files changed, 102 insertions(+)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 8dab711d32..8a477a1548 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -78,6 +78,8 @@ module MOM_CoriolisAdv
   integer :: id_hf_gKEu_2d = -1, id_hf_gKEv_2d = -1
   ! integer :: id_hf_rvxu    = -1, id_hf_rvxv    = -1
   integer :: id_hf_rvxu_2d = -1, id_hf_rvxv_2d = -1
+  integer :: id_h_gKEu = -1, id_h_gKEv = -1
+  integer :: id_h_rvxu = -1, id_h_rvxv = -1
   !>@}
 end type CoriolisAdv_CS
 
@@ -225,6 +227,10 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
   ! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option.
   ! The code is retained for degugging purposes in the future.
 
+! Diagnostics for thickness multiplied momentum budget terms
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_gKEu, h_rvxv ! [L2 T-2 ~> m2 s-2].
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_gKEv, h_rvxu ! [L2 T-2 ~> m2 s-2].
+
 ! To work, the following fields must be set outside of the usual
 ! is to ie range before this subroutine is called:
 !   v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2),
@@ -918,6 +924,38 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
       call post_data(CS%id_hf_rvxu_2d, hf_rvxu_2d, CS%diag)
       deallocate(hf_rvxu_2d)
     endif
+
+    ! Diagnostics for layer thickness multiplied terms
+    if (CS%id_h_gKEu > 0) then
+      h_gKEu(:,:,:) = 0.0
+      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+        h_gKEu(I,j,k) = AD%gradKEu(I,j,k) * AD%diag_hu(I,j,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_h_gKEu, h_gKEu, CS%diag)
+    endif
+    if (CS%id_h_gKEv > 0) then
+      h_gKEv(:,:,:) = 0.0
+      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+        h_gKEv(i,J,k) = AD%gradKEv(i,J,k) * AD%diag_hv(i,J,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_h_gKEv, h_gKEv, CS%diag)
+    endif
+
+    if (CS%id_h_rvxv > 0) then
+      h_rvxv(:,:,:) = 0.0
+      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+        h_rvxv(I,j,k) = AD%rv_x_v(I,j,k) * AD%diag_hu(I,j,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_h_rvxv, h_rvxv, CS%diag)
+    endif
+    if (CS%id_h_rvxu > 0) then
+      h_rvxu(:,:,:) = 0.0
+      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+        h_rvxu(i,J,k) = AD%rv_x_u(i,J,k) * AD%diag_hv(i,J,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_h_rvxu, h_rvxu, CS%diag)
+    endif
+
   endif
 
 end subroutine CorAdCalc
@@ -1210,6 +1248,22 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
     call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
   endif
 
+  CS%id_h_gKEu = register_diag_field('ocean_model', 'h_gKEu', diag%axesCuL, Time, &
+     'Thickness Multiplied Zonal Acceleration from Grad. Kinetic Energy', &
+     'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if (CS%id_h_gKEu > 0) then
+    call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
+  endif
+
+  CS%id_h_gKEv = register_diag_field('ocean_model', 'h_gKEv', diag%axesCvL, Time, &
+     'Thickness Multiplied Meridional Acceleration from Grad. Kinetic Energy', &
+     'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if (CS%id_h_gKEv > 0) then
+    call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
+  endif
+
   !CS%id_hf_rvxu = register_diag_field('ocean_model', 'hf_rvxu', diag%axesCvL, Time, &
   !   'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
   !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
@@ -1242,6 +1296,22 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
     call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
   endif
 
+  CS%id_h_rvxu = register_diag_field('ocean_model', 'h_rvxu', diag%axesCvL, Time, &
+     'Thickness Multiplied Meridional Acceleration from Relative Vorticity', &
+     'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if (CS%id_h_rvxu > 0) then
+    call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
+  endif
+
+  CS%id_h_rvxv = register_diag_field('ocean_model', 'h_rvxv', diag%axesCuL, Time, &
+     'Thickness Multiplied Zonal Acceleration from Relative Vorticity', &
+     'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if (CS%id_h_rvxv > 0) then
+    call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
+  endif
+
 end subroutine CoriolisAdv_init
 
 !> Destructor for coriolisadv_cs
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index fc6d97040a..3f71691078 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -187,6 +187,7 @@ module MOM_hor_visc
   integer :: id_grid_Re_Ah = -1, id_grid_Re_Kh   = -1
   integer :: id_diffu     = -1, id_diffv         = -1
   ! integer :: id_hf_diffu  = -1, id_hf_diffv      = -1
+  integer :: id_h_diffu  = -1, id_h_diffv      = -1
   integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1
   integer :: id_Ah_h      = -1, id_Ah_q          = -1
   integer :: id_Kh_h      = -1, id_Kh_q          = -1
@@ -278,6 +279,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
   real, allocatable, dimension(:,:) :: hf_diffu_2d ! Depth sum of hf_diffu [L T-2 ~> m s-2]
   real, allocatable, dimension(:,:) :: hf_diffv_2d ! Depth sum of hf_diffv [L T-2 ~> m s-2]
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_diffu ! h x diffu [L2 T-2 ~> m2 s-2]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_diffv ! h x diffv [L2 T-2 ~> m2 s-2]
 
   real, dimension(SZIB_(G),SZJB_(G)) :: &
     dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
@@ -1662,6 +1665,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     deallocate(hf_diffv_2d)
   endif
 
+  if (present(ADp) .and. (CS%id_h_diffu > 0)) then
+    h_diffu(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      h_diffu(I,j,k) = diffu(I,j,k) * ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_diffu, h_diffu, CS%diag)
+  endif
+  if (present(ADp) .and. (CS%id_h_diffv > 0)) then
+    h_diffv(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      h_diffv(i,J,k) = diffv(i,J,k) * ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_diffv, h_diffv, CS%diag)
+  endif
+
 end subroutine horizontal_viscosity
 
 !> Allocates space for and calculates static variables used by horizontal_viscosity().
@@ -2378,6 +2396,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
     call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
   endif
 
+  CS%id_h_diffu = register_diag_field('ocean_model', 'h_diffu', diag%axesCuL, Time, &
+      'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if ((CS%id_h_diffu > 0) .and. (present(ADp))) then
+    call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
+  endif
+
+  CS%id_h_diffv = register_diag_field('ocean_model', 'h_diffv', diag%axesCvL, Time, &
+      'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if ((CS%id_h_diffv > 0) .and. (present(ADp))) then
+    call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
+  endif
+
   if (CS%biharmonic) then
     CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time,    &
         'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, &

From 1616df4c25157ccda1aad183bb029998da994871 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 12:24:06 -0500
Subject: [PATCH 11/20] Thickness multiplied diagnostics for stress terms

---
 .../vertical/MOM_vert_friction.F90            | 37 ++++++++++++++++++-
 1 file changed, 36 insertions(+), 1 deletion(-)

diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 7e2c2c6926..1742dfd9d5 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -125,6 +125,7 @@ module MOM_vert_friction
   integer :: id_taux_bot = -1, id_tauy_bot = -1
   integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1
   ! 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
   !>@}
 
@@ -213,6 +214,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
   real, allocatable, dimension(:,:) :: hf_du_dt_visc_2d ! Depth sum of hf_du_dt_visc [L T-2 ~> m s-2]
   real, allocatable, dimension(:,:) :: hf_dv_dt_visc_2d ! Depth sum of hf_dv_dt_visc [L T-2 ~> m s-2]
 
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_du_dt_visc ! h x du_dt_visc [L2 T-2 ~> m2 s-2]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_dv_dt_visc ! h x h_dv_dt_visc [L2 T-2 ~> m2 s-2]
+
   logical :: do_i(SZIB_(G))
   logical :: DoStokesMixing
 
@@ -499,6 +503,21 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
     deallocate(hf_dv_dt_visc_2d)
   endif
 
+  if (CS%id_h_du_dt_visc > 0) then
+    h_du_dt_visc(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      h_du_dt_visc(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_du_dt_visc, h_du_dt_visc, CS%diag)
+  endif
+  if (CS%id_h_dv_dt_visc > 0) then
+    h_dv_dt_visc(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      h_dv_dt_visc(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_dv_dt_visc, h_dv_dt_visc, CS%diag)
+  endif
+
 end subroutine vertvisc
 
 !> Calculate the fraction of momentum originally in a layer that remains
@@ -1833,7 +1852,23 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
       conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_dv_dt_visc_2d > 0) then
     call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
+  endif
+
+  CS%id_h_du_dt_visc = register_diag_field('ocean_model', 'h_du_dt_visc', diag%axesCuL, Time, &
+      'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if ((CS%id_intz_diffu_2d > 0) .and. (present(ADp))) then
+    call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+  endif
+
+  CS%id_h_dv_dt_visc = register_diag_field('ocean_model', 'h_dv_dt_visc', diag%axesCvL, Time, &
+      'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if ((CS%id_h_dv_dt_visc > 0) .and. (present(ADp))) then
+    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
 
   if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &

From c987de7fcb701dbb6dfb7ed7eb26727e0eca4037 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 12:30:34 -0500
Subject: [PATCH 12/20] typo correct

---
 src/parameterizations/vertical/MOM_vert_friction.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 1742dfd9d5..b0949eb60a 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -1858,7 +1858,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   CS%id_h_du_dt_visc = register_diag_field('ocean_model', 'h_du_dt_visc', diag%axesCuL, Time, &
       'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if ((CS%id_intz_diffu_2d > 0) .and. (present(ADp))) then
+  if (CS%id_intz_diffu_2d > 0) then
     call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
   endif
@@ -1866,7 +1866,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   CS%id_h_dv_dt_visc = register_diag_field('ocean_model', 'h_dv_dt_visc', diag%axesCvL, Time, &
       'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if ((CS%id_h_dv_dt_visc > 0) .and. (present(ADp))) then
+  if (CS%id_h_dv_dt_visc > 0) then
     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

From 6779f46ca1af0530fd2bd90a22dff35dc7339868 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 12:54:35 -0500
Subject: [PATCH 13/20] More diagnostics

---
 src/core/MOM_dynamics_split_RK2.F90 | 53 +++++++++++++++++++++++++++++
 1 file changed, 53 insertions(+)

diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index a89514cde4..a82ddaf95b 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -160,14 +160,17 @@ module MOM_dynamics_split_RK2
   integer :: id_PFu    = -1, id_PFv    = -1
   integer :: id_CAu    = -1, id_CAv    = -1
   ! integer :: id_hf_PFu    = -1, id_hf_PFv    = -1
+  integer :: id_h_PFu    = -1, id_h_PFv    = -1
   integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -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
 
   ! Split scheme only.
   integer :: id_uav        = -1, id_vav        = -1
   integer :: id_u_BT_accel = -1, id_v_BT_accel = -1
   ! integer :: id_hf_u_BT_accel    = -1, id_hf_v_BT_accel    = -1
+  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
   !>@}
 
@@ -336,6 +339,11 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     hf_CAu_2d, hf_CAv_2d, & ! Depth integeral of hf_CAu, hf_CAv [L T-2 ~> m s-2].
     hf_u_BT_accel_2d, hf_v_BT_accel_2d ! Depth integeral of hf_u_BT_accel, hf_v_BT_accel
 
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
+    h_PFu, h_CAu, h_u_BT_accel ! [L2 T-2 ~> m2 s-2].
+  real, dimension(SZI_(G),SZJB_(G)),SZK_(GV) :: &
+    h_PFv, h_CAv, h_v_BT_accel ! [L2 T-2 ~> m2 s-2].
+
   real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
 
   logical :: dyn_p_surf
@@ -916,6 +924,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     deallocate(hf_PFv_2d)
   endif
 
+  if (CS%id_h_PFu > 0) then
+    h_PFu(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      h_PFu(I,j,k) = CS%PFu(I,j,k) * CS%ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_PFu, h_PFu, CS%diag)
+  endif
+  if (CS%id_h_PFv > 0) then
+    h_PFv(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      h_PFv(i,J,k) = CS%PFv(i,J,k) * CS%ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_PFv, h_PFv, CS%diag)
+  endif
+
   !if (CS%id_hf_CAu > 0) then
   !  allocate(hf_CAu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
   !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
@@ -949,6 +972,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     deallocate(hf_CAv_2d)
   endif
 
+  if (CS%id_h_CAu > 0) then
+    h_CAu(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      h_CAu(I,j,k) = CS%CAu(I,j,k) * CS%ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_CAu, h_CAu, CS%diag)
+  endif
+  if (CS%id_h_CAv > 0) then
+    h_CAv(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      h_CAv(i,J,k) = CS%CAv(i,J,k) * CS%ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_CAv, h_CAv, CS%diag)
+  endif
+
   !if (CS%id_hf_u_BT_accel > 0) then
   !  allocate(hf_u_BT_accel(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
   !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
@@ -982,6 +1020,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     deallocate(hf_v_BT_accel_2d)
   endif
 
+  if (CS%id_h_u_BT_accel > 0) then
+    h_u_BT_accel(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      h_u_BT_accel(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_u_BT_accel, h_u_BT_accel, CS%diag)
+  endif
+  if (CS%id_h_v_BT_accel > 0) then
+    h_v_BT_accel(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      h_v_BT_accel(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_h_v_BT_accel, h_v_BT_accel, CS%diag)
+  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)

From b8c1717e5b7987cd3e4816d3c8320e1bc8e06200 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 18:42:43 -0500
Subject: [PATCH 14/20] Thickness multiplied momentum budget terms

---
 src/core/MOM_CoriolisAdv.F90                  | 12 ++---
 src/core/MOM_dynamics_split_RK2.F90           | 44 ++++++++++++++++---
 src/diagnostics/MOM_diagnostics.F90           | 43 +++++++++++++++++-
 .../vertical/MOM_vert_friction.F90            |  4 +-
 4 files changed, 86 insertions(+), 17 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 8a477a1548..f081754515 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -1229,7 +1229,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
   !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
   !if (CS%id_hf_gKEv > 0) then
   !  call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
-  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
 
   CS%id_hf_gKEu_2d = register_diag_field('ocean_model', 'hf_gKEu_2d', diag%axesCu1, Time, &
@@ -1245,7 +1245,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'm-1 s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_gKEv_2d > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   endif
 
   CS%id_h_gKEu = register_diag_field('ocean_model', 'h_gKEu', diag%axesCuL, Time, &
@@ -1261,7 +1261,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_h_gKEv > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
   !CS%id_hf_rvxu = register_diag_field('ocean_model', 'hf_rvxu', diag%axesCvL, Time, &
@@ -1269,7 +1269,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
   !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
   !if (CS%id_hf_rvxu > 0) then
   !  call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
-  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
 
   !CS%id_hf_rvxv = register_diag_field('ocean_model', 'hf_rvxv', diag%axesCuL, Time, &
@@ -1285,7 +1285,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'm-1 s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   endif
 
   CS%id_hf_rvxv_2d = register_diag_field('ocean_model', 'hf_rvxv_2d', diag%axesCu1, Time, &
@@ -1301,7 +1301,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_h_rvxu > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
   CS%id_h_rvxv = register_diag_field('ocean_model', 'h_rvxv', diag%axesCuL, Time, &
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index a82ddaf95b..2a1ef5e04e 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -341,7 +341,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
 
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
     h_PFu, h_CAu, h_u_BT_accel ! [L2 T-2 ~> m2 s-2].
-  real, dimension(SZI_(G),SZJB_(G)),SZK_(GV) :: &
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
     h_PFv, h_CAv, h_v_BT_accel ! [L2 T-2 ~> m2 s-2].
 
   real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
@@ -1419,7 +1419,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, &
   !    'Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', &
   !    v_extensive=.true., conversion=US%L_T2_to_m_s2)
-  !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, &
   !    'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
@@ -1429,7 +1429,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, &
   !    'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', &
   !    v_extensive=.true., conversion=US%L_T2_to_m_s2)
-  !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, Time, &
       'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', &
@@ -1439,7 +1439,17 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, Time, &
       'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', &
       conversion=US%L_T2_to_m_s2)
-  if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
+
+  CS%id_h_PFu = register_diag_field('ocean_model', 'h_PFu', diag%axesCuL, Time, &
+      'Thickness Multiplied Zonal Pressure Force Acceleration', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if(CS%id_h_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+
+  CS%id_h_PFv = register_diag_field('ocean_model', 'h_PFv', diag%axesCvL, Time, &
+      'Thickness Multiplied Meridional Pressure Force Acceleration', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if(CS%id_h_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, &
       'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
@@ -1449,7 +1459,17 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, Time, &
       'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', &
       conversion=US%L_T2_to_m_s2)
-  if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
+
+  CS%id_h_CAu = register_diag_field('ocean_model', 'h_CAu', diag%axesCuL, Time, &
+      'Thickness Multiplied Zonal Coriolis and Advective Acceleration', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if(CS%id_h_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+
+  CS%id_h_CAv = register_diag_field('ocean_model', 'h_CAv', diag%axesCvL, Time, &
+      'Thickness Multiplied Meridional Coriolis and Advective Acceleration', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if(CS%id_h_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, &
       'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s)
@@ -1469,7 +1489,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, &
   !    'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', &
   !    v_extensive=.true., conversion=US%L_T2_to_m_s2)
-  !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, Time, &
       'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', &
@@ -1479,7 +1499,17 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, Time, &
       'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', &
       conversion=US%L_T2_to_m_s2)
-  if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
+
+  CS%id_h_u_BT_accel = register_diag_field('ocean_model', 'h_u_BT_accel', diag%axesCuL, Time, &
+      'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if(CS%id_h_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+
+  CS%id_h_v_BT_accel = register_diag_field('ocean_model', 'h_v_BT_accel', diag%axesCvL, Time, &
+      'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', 'm2 s-2', &
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if(CS%id_h_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hv,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/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 47d322dfa0..c19a64fac2 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -116,6 +116,7 @@ module MOM_diagnostics
   integer :: id_e              = -1, id_e_D            = -1
   integer :: id_du_dt          = -1, id_dv_dt          = -1
   ! integer :: id_hf_du_dt       = -1, id_hf_dv_dt       = -1
+  integer :: id_h_du_dt       = -1, id_h_dv_dt       = -1
   integer :: id_hf_du_dt_2d    = -1, id_hf_dv_dt_2d    = -1
   integer :: id_col_ht         = -1, id_dh_dt          = -1
   integer :: id_KE             = -1, id_dKEdt          = -1
@@ -248,6 +249,9 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
   real, allocatable, dimension(:,:) :: &
     hf_du_dt_2d, hf_dv_dt_2d ! z integeral of hf_du_dt, hf_dv_dt [L T-2 ~> m s-2].
 
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: h_du_dt ! h x dudt  [L2 T-2 ~> m2 s-2]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: h_dv_dt ! h x dvdt [L2 T-2 ~> m2 s-2]
+
   ! tmp array for surface properties
   real :: surface_field(SZI_(G),SZJ_(G))
   real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS [R L2 T-2 ~> Pa]
@@ -325,6 +329,21 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
       deallocate(hf_dv_dt_2d)
     endif
 
+    if (CS%id_h_du_dt > 0) then
+      h_du_dt(:,:,:) = 0.0
+      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+        h_du_dt(I,j,k) = CS%du_dt(I,j,k) * ADp%diag_hu(I,j,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_h_du_dt, h_du_dt, CS%diag)
+    endif
+    if (CS%id_h_dv_dt > 0) then
+      h_dv_dt(:,:,:) = 0.0
+      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+        h_dv_dt(i,J,k) = CS%dv_dt(i,J,k) * ADp%diag_hv(i,J,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_h_dv_dt, h_dv_dt, CS%diag)
+    endif
+
     call diag_restore_grids(CS%diag)
 
     call calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS)
@@ -1767,7 +1786,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
   !    call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
   !    call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
   !  endif
-  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
 
   CS%id_hf_du_dt_2d = register_diag_field('ocean_model', 'hf_dudt_2d', diag%axesCu1, Time, &
@@ -1787,7 +1806,27 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
       call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
     endif
-    call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
+  endif
+
+  CS%id_h_du_dt = register_diag_field('ocean_model', 'h_du_dt', diag%axesCuL, Time, &
+      'Thickness Multiplied Zonal Acceleration', 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if (CS%id_h_du_dt > 0) then
+    if (.not.associated(CS%du_dt)) then
+      call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz)
+      call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS)
+    endif
+    call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+  endif
+
+  CS%id_h_dv_dt = register_diag_field('ocean_model', 'h_dv_dt', diag%axesCvL, Time, &
+      'Thickness Multiplied Meridional Acceleration', 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
+  if (CS%id_h_dv_dt > 0) then
+    if (.not.associated(CS%dv_dt)) then
+      call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
+      call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
+    endif
+    call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
   ! layer thickness variables
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index b0949eb60a..15fb6f238b 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -1836,7 +1836,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   !if (CS%id_hf_dv_dt_visc > 0) then
   !  call safe_alloc_ptr(CS%hf_dv_dt_visc,isd,ied,JsdB,JedB,nz)
   !  call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
-  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
 
   CS%id_hf_du_dt_visc_2d = register_diag_field('ocean_model', 'hf_du_dt_visc_2d', diag%axesCu1, Time, &
@@ -1858,7 +1858,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   CS%id_h_du_dt_visc = register_diag_field('ocean_model', 'h_du_dt_visc', diag%axesCuL, Time, &
       'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if (CS%id_intz_diffu_2d > 0) then
+  if (CS%id_h_du_dt_visc > 0) then
     call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
   endif

From 3dd66f5630c9cf18be3d8ceb683a3795665863d1 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Mon, 26 Apr 2021 19:57:09 -0400
Subject: [PATCH 15/20] remove duplicate initialization of array

---
 src/parameterizations/lateral/MOM_hor_visc.F90 | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index d9aac0cc6c..75b26de700 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -278,8 +278,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2]
     boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim]
 
-  real, allocatable, dimension(:,:) :: hf_diffu_2d ! Depth sum of hf_diffu [L T-2 ~> m s-2]
-  real, allocatable, dimension(:,:) :: hf_diffv_2d ! Depth sum of hf_diffv [L T-2 ~> m s-2]
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_diffu ! h x diffu [L2 T-2 ~> m2 s-2]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_diffv ! h x diffv [L2 T-2 ~> m2 s-2]
 

From 3f5682fa8381b3b6728162b904f848ef5a4672de Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Tue, 27 Apr 2021 19:44:05 -0400
Subject: [PATCH 16/20] remove conflicts

---
 src/tracer/nw2_tracers.F90 | 314 -------------------------------------
 1 file changed, 314 deletions(-)
 delete mode 100644 src/tracer/nw2_tracers.F90

diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90
deleted file mode 100644
index e75c5c5d38..0000000000
--- a/src/tracer/nw2_tracers.F90
+++ /dev/null
@@ -1,314 +0,0 @@
-!> 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
-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
-
-!> The control structure for the nw2_tracers package
-type, public :: nw2_tracers_CS ; private
-  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, 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
-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_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()
-  logical :: do_nw2
-  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
-
-  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, "")
-  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))
-
-  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)
-    ! 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.)
-    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
-  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
-  character(len=8)  :: var_name ! The variable's name.
-  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
-
-  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
-    endif ! restart
-  enddo ! Tracer loop
-
-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) ) * 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)
-      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,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 (2) ! -z/L
-    nw2_tracer_dist = -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 709165b62dd2d6b0a1541411a9db9cea51674ca4 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Tue, 27 Apr 2021 20:23:44 -0400
Subject: [PATCH 17/20] remove conflicts

---
 src/tracer/MOM_tracer_flow_control.F90 | 24 +++++-------------------
 1 file changed, 5 insertions(+), 19 deletions(-)

diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90
index 60ff9c981c..26ef197ae2 100644
--- a/src/tracer/MOM_tracer_flow_control.F90
+++ b/src/tracer/MOM_tracer_flow_control.F90
@@ -62,8 +62,6 @@ 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
 
@@ -86,7 +84,6 @@ 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()
@@ -101,7 +98,6 @@ 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
 
@@ -210,9 +206,6 @@ 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
@@ -256,8 +249,7 @@ 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
 
@@ -336,8 +328,6 @@ 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
 
@@ -496,11 +486,8 @@ 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, &
@@ -545,10 +532,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
@@ -789,7 +776,6 @@ 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

From f62e5b7907dc44d0fd8c3235767e58f0bb54ff41 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Wed, 28 Apr 2021 14:04:57 -0400
Subject: [PATCH 18/20] conflict resolve

---
 .gitmodules   | 3 ---
 pkg/CVMix-src | 1 -
 2 files changed, 4 deletions(-)
 delete mode 160000 pkg/CVMix-src

diff --git a/.gitmodules b/.gitmodules
index 637f1188ed..e06031129d 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,6 +1,3 @@
-[submodule "pkg/CVMix-src"]
-	path = pkg/CVMix-src
-	url = https://github.com/CVMix/CVMix-src.git
 [submodule "pkg/GSW-Fortran"]
 	path = pkg/GSW-Fortran
 	url = https://github.com/TEOS-10/GSW-Fortran.git
diff --git a/pkg/CVMix-src b/pkg/CVMix-src
deleted file mode 160000
index 534fc38e75..0000000000
--- a/pkg/CVMix-src
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 534fc38e759fcb8a2563fa0dc4c0bbf81f758db9

From 26ce31d1edc4d796bce96108db0df76bc6d2c603 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Tue, 4 May 2021 13:45:14 -0400
Subject: [PATCH 19/20] submodule correct

---
 .gitmodules | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/.gitmodules b/.gitmodules
index e06031129d..637f1188ed 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,3 +1,6 @@
+[submodule "pkg/CVMix-src"]
+	path = pkg/CVMix-src
+	url = https://github.com/CVMix/CVMix-src.git
 [submodule "pkg/GSW-Fortran"]
 	path = pkg/GSW-Fortran
 	url = https://github.com/TEOS-10/GSW-Fortran.git

From 0b9404111d16ce3216e5a089c04ddf19518f93cd Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Tue, 4 May 2021 13:48:06 -0400
Subject: [PATCH 20/20] cvmix correct

---
 pkg/CVMix-src | 1 +
 1 file changed, 1 insertion(+)
 create mode 160000 pkg/CVMix-src

diff --git a/pkg/CVMix-src b/pkg/CVMix-src
new file mode 160000
index 0000000000..9423197f89
--- /dev/null
+++ b/pkg/CVMix-src
@@ -0,0 +1 @@
+Subproject commit 9423197f894112edfcb1502245f7d7b873d551f9