MOM_dynamics_legacy_split.F90

Go to the documentation of this file.

module mom_dynamics_legacy_split

!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg and Alistair Adcroft, 1996-2012                 *
!*                                                                     *
!*    This file contains code that does the time-stepping of the       *
!*  adiabatic dynamic core, in this case with mode-splitting between   *
!*  the baroclinic and barotropic modes and a pseudo-second order      *
!*  Runge-Kutta time stepping scheme for the baroclinic momentum       *
!*  equation and a forward-backward coupling between the baroclinic    *
!*  momentum and continuity equations.  This split time-stepping       *
!*  scheme is described in detail in Hallberg (JCP, 1997), with the    *
!*  primary issues related to exact tracer conservation and how to     *
!*  ensure consistency between the barotropic and layered estimates    *
!*  of the free surface height described carefully in Hallberg and     *
!*  Adcroft (Ocean Modelling, 2009).  This was the time stepping code  *
!*  that is used for most GOLD applications, including GFDL's ESM2G    *
!*  Earth system model, and all of the examples provided with the      *
!*  MOM code (although several of these solutions are routinely        *
!*  verified by comparison with the slower unsplit schemes).           *
!*                                                                     *
!*    The subroutine step_MOM_dyn_legacy_split actually does the time  *
!*  stepping, while register_restarts_dyn_legacy_split sets the fields *
!*  that are found in a full restart file with this scheme, and        *
!*  initialize_dyn_legacy_split initializes the cpu clocks that are    *                                      *
!*  used in this module.  For largely historical reasons, this module  *
!*  does not have its own control structure, but shares the same       *
!*  control structure with MOM.F90 and the other MOM_dynamics_...      *
!*  modules.                                                           *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h. *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, CoriolisBu                            *
!*    j+1  > o > o >   At ^:  v, PFv, CAv, vh, diffv, tauy, vbt, vhtr  *
!*    j    x ^ x ^ x   At >:  u, PFu, CAu, uh, diffu, taux, ubt, uhtr  *
!*    j    > o > o >   At o:  h, bathyT, eta, T, S, tr                 *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1                                                  *
!*           i  i+1                                                    *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**


use mom_variables, only : vertvisc_type, thermo_var_ptrs
use mom_variables, only : bt_cont_type, alloc_bt_cont_type, dealloc_bt_cont_type
use mom_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs
use mom_forcing_type, only : forcing

use mom_checksum_packages, only : mom_thermo_chksum, mom_state_chksum, mom_accel_chksum
use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock, only : clock_component, clock_subcomponent
use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
use mom_diag_mediator, only : diag_mediator_init, enable_averaging
use mom_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
use mom_diag_mediator, only : register_diag_field, register_static_field
use mom_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids
use mom_domains, only : mom_domains_init, pass_var, pass_vector
use mom_domains, only : pass_var_start, pass_var_complete
use mom_domains, only : pass_vector_start, pass_vector_complete
use mom_domains, only : to_south, to_west, to_all, cgrid_ne, scalar_pair
use mom_debugging, only : hchksum, uvchksum
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
use mom_error_handler, only : mom_set_verbosity
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_get_input, only : directories
use mom_io, only : mom_io_init, vardesc, var_desc
use mom_restart, only : register_restart_field, query_initialized, save_restart
use mom_restart, only : restart_init, mom_restart_cs
use mom_time_manager, only : time_type, set_time, time_type_to_real, operator(+)
use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)

use mom_ale, only : ale_cs
use mom_legacy_barotropic, only : legacy_barotropic_init, legacy_btstep
use mom_legacy_barotropic, only : legacy_btcalc, legacy_bt_mass_source
use mom_legacy_barotropic, only : register_legacy_barotropic_restarts
use mom_legacy_barotropic, only : legacy_set_dtbt, legacy_barotropic_cs
use mom_continuity, only : continuity, continuity_init, continuity_cs
use mom_coriolisadv, only : coradcalc, coriolisadv_init, coriolisadv_cs
use mom_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_cs
use mom_debugging, only : check_redundant
use mom_grid, only : ocean_grid_type
use mom_hor_index, only : hor_index_type
use mom_hor_visc, only : horizontal_viscosity, hor_visc_init, hor_visc_cs
use mom_interface_heights, only : find_eta
use mom_lateral_mixing_coeffs, only : varmix_cs
use mom_meke_types, only : meke_type
use mom_open_boundary, only : ocean_obc_type
use mom_open_boundary, only : radiation_open_bdry_conds
use mom_open_boundary, only : open_boundary_zero_normal_flow
use mom_boundary_update, only : update_obc_data, update_obc_cs
use mom_pressureforce, only : pressureforce, pressureforce_init, pressureforce_cs
use mom_tidal_forcing, only : tidal_forcing_init, tidal_forcing_cs
use mom_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant
use mom_vert_friction, only : vertvisc_limit_vel, vertvisc_init, vertvisc_cs
use mom_set_visc, only : set_viscous_ml, set_visc_cs
use mom_verticalgrid, only : verticalgrid_type, get_thickness_units
use mom_verticalgrid, only : get_flux_units, get_tr_flux_units

implicit none ; private

#include <MOM_memory.h>

type, public :: mom_dyn_legacy_split_cs ; private
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
    cau, &    ! CAu = f*v - u.grad(u) in m s-2.
    pfu, &    ! PFu = -dM/dx, in m s-2.
    diffu, &  ! Zonal acceleration due to convergence of the along-isopycnal
              ! stress tensor, in m s-2.
    visc_rem_u, & ! Both the fraction of the zonal momentum originally in a
              ! layer that remains after a time-step of viscosity, and the
              ! fraction of a time-step's worth of a barotropic acceleration
              ! that a layer experiences after viscosity is applied.
              ! Nondimensional between 0 (at the bottom) and 1 (far above).
    u_accel_bt ! The layers' zonal accelerations due to the difference between
              ! the barotropic accelerations and the baroclinic accelerations
              ! that were fed into the barotopic calculation, in m s-2.
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
    cav, &    ! CAv = -f*u - u.grad(v) in m s-2.
    pfv, &    ! PFv = -dM/dy, in m s-2.
    diffv, &  ! Meridional acceleration due to convergence of the
              ! along-isopycnal stress tensor, in m s-2.
    visc_rem_v, & ! Both the fraction of the meridional momentum originally in
              ! a layer that remains after a time-step of viscosity, and the
              ! fraction of a time-step's worth of a barotropic acceleration
              ! that a layer experiences after viscosity is applied.
              ! Nondimensional between 0 (at the bottom) and 1 (far above).
    v_accel_bt ! The layers' meridional accelerations due to the difference between
              ! the barotropic accelerations and the baroclinic accelerations
              ! that were fed into the barotopic calculation, in m s-2.

! The following variables are only used with the split time stepping scheme.
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    eta       ! Instantaneous free surface height, in m.
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_av
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_av
    ! u_av and v_av are the layer velocities with the vertical mean replaced by
    ! the time-mean barotropic velocity over a baroclinic timestep, in m s-1.
  real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_)  :: h_av
    ! The arithmetic mean of two successive layer thicknesses, in m or kg m-2.
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    eta_pf    ! The instantaneous SSH used in calculating PFu and PFv, in m.
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: uhbt
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: vhbt
    ! uhbt and vhbt are the average volume or mass fluxes determined by the
    ! barotropic solver in m3 s-1 or kg s-1.  uhbt and vhbt should (roughly?)
    ! equal the verticals sum of uh and vh, respectively.
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: uhbt_in
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: vhbt_in
    ! uhbt_in and vhbt_in are the vertically summed transports from based on
    ! the final thicknessses and velocities from the previous dynamics time
    ! step, both in units of m3 s-1 or kg s-1.
  real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce
      ! pbce times eta gives the baroclinic pressure anomaly in each layer due
      ! to free surface height anomalies.  pbce has units of m2 H-1 s-2.

  real, pointer, dimension(:,:) :: taux_bot => null(), tauy_bot => null()
    ! The frictional bottom stresses from the ocean to the seafloor, in Pa.
  type(bt_cont_type), pointer :: bt_cont => null()
                              ! A structure with elements that describe the
                              ! effective summed open face areas as a function
                              ! of barotropic flow.

  ! This is to allow the previous, velocity-based coupling with between the
  ! baroclinic and barotropic modes.
  logical :: flux_bt_coupling  ! If true, use volume fluxes, not velocities,
                               ! to couple the baroclinic and barotropic modes.
  logical :: bt_use_layer_fluxes ! If true, use the summed layered fluxes plus
                               ! an adjustment due to a changed barotropic
                               ! velocity in the barotropic continuity equation.
  logical :: split_bottom_stress  ! If true, provide the bottom stress
                               ! calculated by the vertical viscosity to the
                               ! barotropic solver.
  logical :: readjust_bt_trans ! If true, readjust the barotropic transport of
                               ! the input velocities to agree with CS%uhbt_in
                               ! and CS%vhbt_in after the diabatic step.
  logical :: readjust_velocity ! A flag that varies with time that determines
                               ! whether the velocities currently need to be
                               ! readjusted to agree with CS%uhbt_in and
                               ! CS%vhbt_in.  This is only used if
                               ! CS%readjust_BT_trans is true.
  logical :: calc_dtbt         ! If true, calculate the barotropic time-step
                               ! dynamically.

  real    :: be              ! A nondimensional number from 0.5 to 1 that controls
                             ! the backward weighting of the time stepping scheme.
  real    :: begw            ! A nondimensional number from 0 to 1 that controls
                             ! the extent to which the treatment of gravity waves
                             ! is forward-backward (0) or simulated backward
                             ! Euler (1).  0 is almost always used.
  logical :: debug           ! If true, write verbose checksums for debugging purposes.

  logical :: module_is_initialized = .false.

  integer :: id_uh = -1, id_vh = -1
  integer :: id_pfu = -1, id_pfv = -1, id_cau = -1, id_cav = -1

! Split scheme only.
  integer :: id_uav = -1, id_vav = -1
  integer :: id_du_adj = -1, id_dv_adj = -1, id_du_adj2 = -1, id_dv_adj2 = -1
  integer :: id_u_bt_accel = -1, id_v_bt_accel = -1

  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                                   ! timing of diagnostic output.
  type(accel_diag_ptrs), pointer :: adp ! A structure pointing to the various
                                   ! accelerations in the momentum equations,
                                   ! which can later be used to calculate
                                   ! derived diagnostics like energy budgets.
  type(cont_diag_ptrs), pointer :: cdp ! A structure with pointers to various
                                   ! terms in the continuity equations,
                                   ! which can later be used to calculate
                                   ! derived diagnostics like energy budgets.
! The remainder of the structure is pointers to child subroutines' control strings.
  type(hor_visc_cs), pointer :: hor_visc_csp => null()
  type(continuity_cs), pointer :: continuity_csp => null()
  type(coriolisadv_cs), pointer :: coriolisadv_csp => null()
  type(pressureforce_cs), pointer :: pressureforce_csp => null()
  type(legacy_barotropic_cs), pointer :: barotropic_csp => null()
  type(vertvisc_cs), pointer :: vertvisc_csp => null()
  type(set_visc_cs), pointer :: set_visc_csp => null()
  type(ocean_obc_type), pointer :: obc => null() ! A pointer to an open boundary
     ! condition type that specifies whether, where, and  what open boundary
     ! conditions are used.  If no open BCs are used, this pointer stays
     ! nullified.  Flather OBCs use open boundary_CS as well.
  type(update_obc_cs),    pointer :: update_obc_csp => null()
  type(tidal_forcing_cs), pointer :: tides_csp => null()

! This is a copy of the pointer in the top-level control structure.
  type(ale_cs), pointer :: ale_csp => null()

end type mom_dyn_legacy_split_cs

public step_mom_dyn_legacy_split, register_restarts_dyn_legacy_split
public adjustments_dyn_legacy_split
public initialize_dyn_legacy_split, end_dyn_legacy_split

integer :: id_clock_cor, id_clock_pres, id_clock_vertvisc
integer :: id_clock_horvisc, id_clock_mom_update
integer :: id_clock_continuity, id_clock_thick_diff
integer :: id_clock_btstep, id_clock_btcalc, id_clock_btforce
integer :: id_clock_pass, id_clock_pass_init

contains

! =============================================================================

subroutine step_mom_dyn_legacy_split(u, v, h, tv, visc, &
                 Time_local, dt, fluxes, p_surf_begin, p_surf_end, &
                 dt_since_flux, dt_therm, uh, vh, uhtr, vhtr, eta_av, &
                 G, GV, CS, calc_dtbt, VarMix, MEKE)
  type(ocean_grid_type),      intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),    intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                      target, intent(inout) :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                      target, intent(inout) :: v    !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                              intent(inout) :: h   !< Layer thicknesses, in H (usually m or kg m-2).
  type(thermo_var_ptrs),      intent(in)    :: tv   !< A structure pointing to various.
                                                    !! thermodynamic variables.
  type(vertvisc_type),        intent(inout) :: visc !< A structure containing vertical viscosities,
                                                    !! bottom drag viscosities, and related fields.
  type(time_type),            intent(in)    :: Time_local !< The model time at the end
                                                          !! of the time step.
  real,                       intent(in)    :: dt   !< The baroclinic dynamics time step, in s
  type(forcing),              intent(in)    :: fluxes !< A structure containing pointers to any
                                                      !! possible forcing fields.  Unused fields
                                                      !! have NULL ptrs.
  real, dimension(:,:),       pointer       :: p_surf_begin !< A pointer (perhaps NULL) to the
                                                            !! surface pressure at the beginning
                                                            !! of this dynamic step, in Pa.
  real, dimension(:,:),       pointer       :: p_surf_end !< A pointer (perhaps NULL) to the
                                                          !! surface pressure at the end of
                                                          !! this dynamic step, in Pa.
  real,                       intent(in)    :: dt_since_flux !< The elapsed time since fluxes
                                                             !! were applied, in s.
  real,                       intent(in)    :: dt_therm !< The thermodynamic time step, in s.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                      target, intent(inout) :: uh !< The zonal volume or mass transport,
                                                  !! in m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                      target, intent(inout) :: vh !< The meridional volume or mass transport,
                                                  !! in m3 s-1 or kg s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                              intent(inout) :: uhtr !< The accumulated zonal volume or mass
                                                    !! transport since the last tracer advection,
                                                    !! in m3 or kg.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                              intent(inout) :: vhtr !< The accumulated meridional volume or mass
                                                    !! transport since the last tracer advection,
                                                    !! in m3 or kg.
  real, dimension(SZI_(G),SZJ_(G)),          &
                              intent(out)   :: eta_av !< The free surface height or column mass
                                                      !! time-averaged over a time step,
                                                      !! in m or kg m-2.
  type(mom_dyn_legacy_split_cs),             &
                              pointer       :: CS !< The control structure set up by
                                                  !! initialize_dyn_legacy_split.
  logical,                    intent(in)    :: calc_dtbt !< If true, recalculate the
                                                         !! barotropic time step.
  type(varmix_cs),            pointer       :: VarMix !<A pointer to a structure with fields that
                                                      !! specify the spatially variable viscosities.
  type(meke_type),            pointer       :: MEKE   !< A pointer to a structure containing fields
                                                    !! related to the Mesoscale Eddy Kinetic Energy.
! Arguments: u - The zonal velocity, in m s-1.
!  (inout)   v - The meridional velocity, in m s-1.
!  (inout)   h - The layer thicknesses, in m or kg m-2, depending on
!                whether the Boussinesq approximation is made.
!  (in)      tv - a structure pointing to various thermodynamic variables.
!  (inout)   visc - A structure containing vertical viscosities, bottom drag
!                   viscosities, and related fields.
!  (in)      Time_local - The model time at the end of the time step.
!  (in)      dt - The time step in s.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      p_surf_begin - A pointer (perhaps NULL) to the surface pressure
!                     at the beginning of this dynamic step, in Pa.
!  (in)      p_surf_end - A pointer (perhaps NULL) to the surface pressure
!                     at the end of this dynamic step, in Pa.
!  (in)      dt_since_flux - The elapsed time since fluxes were applied, in s.
!  (in)      dt_therm - The thermodynamic time step, in s.
!  (inout)   uh - The zonal volume or mass transport, in m3 s-1 or kg s-1.
!  (inout)   vh - The meridional volume or mass transport, in m3 s-1 or kg s-1.
!  (inout)   uhtr - The accumulated zonal volume or mass transport since the last
!                   tracer advection, in m3 or kg.
!  (inout)   vhtr - The accumulated meridional volume or mass transport since the last
!                   tracer advection, in m3 or kg.
!  (out)     eta_av - The free surface height or column mass time-averaged
!                     over a time step, in m or kg m-2.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure set up by initialize_dyn_legacy_split.
!  (in)      calc_dtbt - If true, recalculate the barotropic time step.
!  (in)      VarMix - A pointer to a structure with fields that specify the
!                     spatially variable viscosities.
!  (inout)   MEKE - A pointer to a structure containing fields related to
!                   the Mesoscale Eddy Kinetic Energy.

  real :: dt_pred   ! The time step for the predictor part of the baroclinic
                    ! time stepping.

  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
    up   ! Predicted zonal velocitiy in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
    vp   ! Predicted meridional velocitiy in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: &
    hp   ! Predicted thickness in m or kg m-2 (H).

  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_bc_accel
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_bc_accel
    ! u_bc_accel and v_bc_accel are the summed baroclinic accelerations of each
    ! layer calculated by the non-barotropic part of the model, both in m s-2.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target :: uh_in
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target :: vh_in
    ! uh_in and vh_in are the zonal or meridional mass transports that would be
    ! obtained using the initial velocities, both in m3 s-1 or kg s-1.
  real, dimension(SZIB_(G),SZJ_(G)) :: uhbt_out
  real, dimension(SZI_(G),SZJB_(G)) :: vhbt_out
    ! uhbt_out and vhbt_out are the vertically summed transports from the
    ! barotropic solver based on its final velocities, both in m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJ_(G)) :: eta_pred
    ! eta_pred is the predictor value of the free surface height or column mass,
    ! in m or kg m-2.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target :: u_adj
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target :: v_adj
    ! u_adj and v_adj are the zonal or meridional velocities after u and v
    ! have been barotropically adjusted so the resulting transports match
    ! uhbt_out and vhbt_out, both in m s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_tmp
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_tmp
    ! u_tmp and v_tmp are temporary velocities used in diagnostic calculations.

  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_old_rad_OBC
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_old_rad_OBC
    ! u_old_rad_OBC and v_old_rad_OBC are the starting velocities, which are
    ! saved for use in the Flather open boundary condition code, both in m s-1.

  real :: Pa_to_eta ! A factor that converts pressures to the units of eta.
  real, pointer, dimension(:,:) :: &
    p_surf => null(), eta_pf_start => null(), &
    taux_bot => null(), tauy_bot => null(), &
    uhbt_in, vhbt_in, eta
  real, pointer, dimension(:,:,:) :: &
    uh_ptr => null(), u_ptr => null(),  vh_ptr => null(), v_ptr => null(), &
    u_init => null(), v_init => null(), & ! Pointers to u and v or u_adj and v_adj.
    u_av, & ! The zonal velocity time-averaged over a time step, in m s-1.
    v_av, & ! The meridional velocity time-averaged over a time step, in m s-1.
    h_av    ! The layer thickness time-averaged over a time step, in m or
            ! kg m-2.
  real :: Idt
  logical :: dyn_p_surf
  logical :: BT_cont_BT_thick ! If true, use the BT_cont_type to estimate the
                              ! relative weightings of the layers in calculating
                              ! the barotropic accelerations.
  integer :: pid_bbl_h, pid_eta_PF, pid_eta, pid_visc
  integer :: pid_h, pid_u, pid_u_av, pid_uh, pid_uhbt_in
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  u_av => cs%u_av ; v_av => cs%v_av ; h_av => cs%h_av
  eta => cs%eta ; uhbt_in => cs%uhbt_in ; vhbt_in => cs%vhbt_in
  idt = 1.0 / dt

  up(:,:,:) = 0.0 ; vp(:,:,:) = 0.0 ; hp(:,:,:) = h(:,:,:)

  if (cs%debug) then
    call mom_state_chksum("Start predictor ", u, v, h, uh, vh, g, gv)
    call check_redundant("Start predictor u ", u, v, g)
    call check_redundant("Start predictor uh ", uh, vh, g)
  endif

  dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
  if (dyn_p_surf) then
    p_surf => p_surf_end
    call safe_alloc_ptr(eta_pf_start,g%isd,g%ied,g%jsd,g%jed)
    eta_pf_start(:,:) = 0.0
  else
    p_surf => fluxes%p_surf
  endif

  if (associated(cs%OBC)) then
    do k=1,nz ; do j=js,je ; do i=is-2,ie+1
      u_old_rad_obc(i,j,k) = u(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=js-2,je+1 ; do i=is,ie
      v_old_rad_obc(i,j,k) = v(i,j,k)
    enddo ; enddo ; enddo
  endif

  if (ASSOCIATED(cs%ADp%du_other)) cs%ADp%du_other(:,:,:) = 0.0
  if (ASSOCIATED(cs%ADp%dv_other)) cs%ADp%dv_other(:,:,:) = 0.0

  bt_cont_bt_thick = .false.
  if (associated(cs%BT_cont)) bt_cont_bt_thick = &
    (associated(cs%BT_cont%h_u) .and. associated(cs%BT_cont%h_v))

  if (cs%split_bottom_stress) then
    taux_bot => cs%taux_bot ; tauy_bot => cs%tauy_bot
  endif

! PFu = d/dx M(h,T,S)
! pbce = dM/deta
  if (cs%begw == 0.0) call enable_averaging(dt, time_local, cs%diag)
  call cpu_clock_begin(id_clock_pres)
  call pressureforce(h, tv, cs%PFu, cs%PFv, g, gv, cs%PressureForce_CSp, &
                     cs%ALE_CSp, p_surf, cs%pbce, cs%eta_PF)
  if (dyn_p_surf) then
    if (gv%Boussinesq) then
      pa_to_eta = 1.0 / (gv%Rho0*gv%g_Earth)
    else
      pa_to_eta = 1.0 / gv%H_to_Pa
    endif
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      eta_pf_start(i,j) = cs%eta_PF(i,j) - pa_to_eta * &
                          (p_surf_begin(i,j) - p_surf_end(i,j))
    enddo ; enddo
  endif
  call cpu_clock_end(id_clock_pres)
  call disable_averaging(cs%diag)

  if (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    pid_eta_pf = pass_var_start(cs%eta_PF, g%Domain)
    pid_eta = pass_var_start(eta, g%Domain)
    if (cs%readjust_velocity) &
      pid_uhbt_in = pass_vector_start(uhbt_in, vhbt_in, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, tv, h, cs%update_OBC_CSp, time_local)
  endif; endif

! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, g, gv, &
                 cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)

! u_bc_accel = CAu + PFu + diffu(u[n-1])
  call cpu_clock_begin(id_clock_btforce)
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
  enddo ; enddo ; enddo
  if (associated(cs%OBC)) then
    call open_boundary_zero_normal_flow(cs%OBC, g, u_bc_accel, v_bc_accel)
  endif
  call cpu_clock_end(id_clock_btforce)

  if (cs%debug) then
    call check_redundant("pre-btstep CS%Ca ", cs%Cau, cs%Cav, g)
    call check_redundant("pre-btstep CS%PF ", cs%PFu, cs%PFv, g)
    call check_redundant("pre-btstep CS%diff ", cs%diffu, cs%diffv, g)
    call check_redundant("pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g)
  endif

  if (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    call pass_var_complete(pid_eta_pf, cs%eta_PF, g%Domain)
    call pass_var_complete(pid_eta, eta, g%Domain)
    if (cs%readjust_velocity) &
      call pass_vector_complete(pid_uhbt_in, uhbt_in, vhbt_in, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  call cpu_clock_begin(id_clock_vertvisc)
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    up(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * u_bc_accel(i,j,k))
  enddo ; enddo ;  enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    vp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * v_bc_accel(i,j,k))
  enddo ; enddo ;  enddo
  call enable_averaging(dt, time_local, cs%diag)
  call set_viscous_ml(u, v, h, tv, fluxes, visc, dt, g, gv, &
                      cs%set_visc_CSp)
  call disable_averaging(cs%diag)

  call vertvisc_coef(up, vp, h, fluxes, visc, dt, g, gv, cs%vertvisc_CSp, cs%OBC)
  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)

  call cpu_clock_begin(id_clock_pass)
  if (g%nonblocking_updates) then
    pid_visc = pass_vector_start(cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                                 to_all+scalar_pair, cgrid_ne)
  else
    call pass_var(cs%eta_PF, g%Domain, complete=.false.)
    call pass_var(eta, g%Domain)
    if (cs%readjust_velocity) call pass_vector(uhbt_in, vhbt_in, g%Domain)
    call pass_vector(cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                     to_all+scalar_pair, cgrid_ne)
  endif
  call cpu_clock_end(id_clock_pass)

  call cpu_clock_begin(id_clock_btcalc)
  ! Calculate the relative layer weights for determining barotropic quantities.
  if (.not.bt_cont_bt_thick) &
    call legacy_btcalc(h, g, gv, cs%barotropic_CSp)
  call legacy_bt_mass_source(h, eta, fluxes, .true., dt_therm, dt_since_flux, &
                      g, gv, cs%barotropic_CSp)
  call cpu_clock_end(id_clock_btcalc)

  if (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    call pass_vector_complete(pid_visc, cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                     to_all+scalar_pair, cgrid_ne)
    call cpu_clock_end(id_clock_pass)
  endif

! u_accel_bt = layer accelerations due to barotropic solver
  if (cs%flux_BT_coupling) then
    call cpu_clock_begin(id_clock_continuity)
    if (cs%readjust_velocity) then
      ! Adjust the input velocites so that their transports match uhbt_out & vhbt_out.
      call continuity(u, v, h, hp, uh_in, vh_in, dt, g, gv, &
                      cs%continuity_CSp, uhbt_in, vhbt_in, cs%OBC, &
                      cs%visc_rem_u, cs%visc_rem_v, u_adj, v_adj, &
                      bt_cont=cs%BT_cont)
      u_init => u_adj ; v_init => v_adj
      if (ASSOCIATED(cs%ADp%du_other)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq
        cs%ADp%du_other(i,j,k) = u_adj(i,j,k) - u(i,j,k)
      enddo ; enddo ; enddo ; endif
      if (ASSOCIATED(cs%ADp%dv_other)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie
        cs%ADp%dv_other(i,j,k) = v_adj(i,j,k) - v(i,j,k)
      enddo ; enddo ; enddo ; endif
      cs%readjust_velocity = .false.
    else
      call continuity(u, v, h, hp, uh_in, vh_in, dt, g, gv, &
                      cs%continuity_CSp, obc=cs%OBC, bt_cont=cs%BT_cont)
!###   call continuity(u, v, h, hp, uh_in, vh_in, dt, G, GV, &
!###                   CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, &
!###                      visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont)
      u_init => u ; v_init => v
    endif
    call cpu_clock_end(id_clock_continuity)

    if (bt_cont_bt_thick) then
      call cpu_clock_begin(id_clock_pass)
      call pass_vector(cs%BT_cont%h_u, cs%BT_cont%h_v, g%Domain, &
                       to_all+scalar_pair, cgrid_ne)
      call cpu_clock_end(id_clock_pass)
      call legacy_btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v)
    endif
    call cpu_clock_begin(id_clock_btstep)
    if (calc_dtbt) call legacy_set_dtbt(g, gv, cs%barotropic_CSp, eta, cs%pbce, cs%BT_cont)
    call legacy_btstep(.true., uh_in, vh_in, eta, dt, u_bc_accel, v_bc_accel, &
                fluxes, cs%pbce, cs%eta_PF, uh, vh, cs%u_accel_bt, &
                cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, &
                cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, &
                uhbt_out = uhbt_out, vhbt_out = vhbt_out, obc = cs%OBC, &
                bt_cont = cs%BT_cont, eta_pf_start = eta_pf_start, &
                taux_bot=taux_bot, tauy_bot=tauy_bot)
    call cpu_clock_end(id_clock_btstep)
  else

    if (associated(cs%BT_cont) .or. cs%BT_use_layer_fluxes) then
      call cpu_clock_begin(id_clock_continuity)
      call continuity(u, v, h, hp, uh_in, vh_in, dt, g, gv, &
                      cs%continuity_CSp, obc=cs%OBC, &
                      visc_rem_u=cs%visc_rem_u, visc_rem_v=cs%visc_rem_v, &
                      bt_cont=cs%BT_cont)
      call cpu_clock_end(id_clock_continuity)
      if (bt_cont_bt_thick) then
        call cpu_clock_begin(id_clock_pass)
        call pass_vector(cs%BT_cont%h_u, cs%BT_cont%h_v, g%Domain, &
                         to_all+scalar_pair, cgrid_ne)
        call cpu_clock_end(id_clock_pass)
        call legacy_btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v)
      endif
    endif

    if (cs%BT_use_layer_fluxes) then
      uh_ptr => uh_in; vh_ptr => vh_in; u_ptr => u; v_ptr => v
    endif

    u_init => u ; v_init => v
    call cpu_clock_begin(id_clock_btstep)
    if (calc_dtbt) call legacy_set_dtbt(g, gv, cs%barotropic_CSp, eta, cs%pbce)
    call legacy_btstep(.false., u, v, eta, dt, u_bc_accel, v_bc_accel, &
                fluxes, cs%pbce, cs%eta_PF, u_av, v_av, cs%u_accel_bt, &
                cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, cs%barotropic_CSp,&
                cs%visc_rem_u, cs%visc_rem_v, obc=cs%OBC, &
                bt_cont = cs%BT_cont, eta_pf_start=eta_pf_start, &
                taux_bot=taux_bot, tauy_bot=tauy_bot, &
                uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr)
    call cpu_clock_end(id_clock_btstep)
  endif

! up = u + dt_pred*( u_bc_accel + u_accel_bt )
  dt_pred = dt * cs%be
  call cpu_clock_begin(id_clock_mom_update)
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    vp(i,j,k) = g%mask2dCv(i,j) * (v_init(i,j,k) + dt_pred * &
                    (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
  enddo ; enddo ;  enddo

  do k=1,nz ; do j=js,je ; do i=isq,ieq
    up(i,j,k) = g%mask2dCu(i,j) * (u_init(i,j,k) + dt_pred  * &
                    (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
  enddo ; enddo ;  enddo
  call cpu_clock_end(id_clock_mom_update)

  if (cs%debug) then
    call uvchksum("Predictor 1 [uv]", up, vp, g%HI,haloshift=0)
    call hchksum(h, "Predictor 1 h", g%HI, haloshift=1, scale=gv%H_to_m)
    call uvchksum("Predictor 1 [uv]h", uh, vh, &
                  g%HI, haloshift=2, scale=gv%H_to_m)
!   call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, haloshift=1)
    call mom_accel_chksum("Predictor accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
             cs%diffu, cs%diffv, g, gv, cs%pbce, cs%u_accel_bt, cs%v_accel_bt)
    call mom_state_chksum("Predictor 1 init", u_init, v_init, h, uh, vh, g, gv, haloshift=2)
    call check_redundant("Predictor 1 up", up, vp, g)
    call check_redundant("Predictor 1 uh", uh, vh, g)
  endif

! up <- up + dt_pred d/dz visc d/dz up
! u_av  <- u_av  + dt_pred d/dz visc d/dz u_av
  call cpu_clock_begin(id_clock_vertvisc)
  call vertvisc_coef(up, vp, h, fluxes, visc, dt_pred, g, gv, cs%vertvisc_CSp, &
                     cs%OBC)
  call vertvisc(up, vp, h, fluxes, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, g, &
                gv, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
  if (g%nonblocking_updates) then
    call cpu_clock_end(id_clock_vertvisc) ; call cpu_clock_begin(id_clock_pass)
    pid_u = pass_vector_start(up, vp, g%Domain)
    call cpu_clock_end(id_clock_pass) ; call cpu_clock_begin(id_clock_vertvisc)
  endif
  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt_pred, g, gv, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)

  call cpu_clock_begin(id_clock_pass)
  call pass_vector(cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                   to_all+scalar_pair, cgrid_ne)
  if (g%nonblocking_updates) then
    call pass_vector_complete(pid_u, up, vp, g%Domain)
  else
    call pass_vector(up, vp, g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

! uh = u_av * h
! hp = h + dt * div . uh
  call cpu_clock_begin(id_clock_continuity)
  call continuity(up, vp, h, hp, uh, vh, dt, g, gv, cs%continuity_CSp, &
                  cs%uhbt, cs%vhbt, cs%OBC, cs%visc_rem_u, &
                  cs%visc_rem_v, u_av, v_av, bt_cont=cs%BT_cont)
  call cpu_clock_end(id_clock_continuity)

  call cpu_clock_begin(id_clock_pass)
  call pass_var(hp, g%Domain)
  if (g%nonblocking_updates) then
    pid_u_av = pass_vector_start(u_av, v_av, g%Domain)
    pid_uh = pass_vector_start(uh(:,:,:), vh(:,:,:), g%Domain)
  else
    call pass_vector(u_av, v_av, g%Domain, complete=.false.)
    call pass_vector(uh(:,:,:), vh(:,:,:), g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

  if (associated(cs%OBC)) then
    call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, dt)
  endif

! h_av = (h + hp)/2
  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
    h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
  enddo ; enddo ; enddo

! The correction phase of the time step starts here.
  call enable_averaging(dt, time_local, cs%diag)

!   Calculate a revised estimate of the free-surface height correction to be
! used in the next call to btstep.  This call is at this point so that
! hp can be changed if CS%begw /= 0.
! eta_cor = ...                 (hidden inside CS%barotropic_CSp)
  call cpu_clock_begin(id_clock_btcalc)
  call legacy_bt_mass_source(hp, eta_pred, fluxes, .false., dt_therm, &
                      dt_since_flux+dt, g, gv, cs%barotropic_CSp)
  call cpu_clock_end(id_clock_btcalc)

  if (cs%begw /= 0.0) then
    ! hp <- (1-begw)*h_in + begw*hp
    ! Back up hp to the value it would have had after a time-step of
    ! begw*dt.  hp is not used again until recalculated by continuity.
    do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
      hp(i,j,k) = (1.0-cs%begw)*h(i,j,k) + cs%begw*hp(i,j,k)
    enddo ; enddo ; enddo

! PFu = d/dx M(hp,T,S)
! pbce = dM/deta
    call cpu_clock_begin(id_clock_pres)
    call pressureforce(hp, tv, cs%PFu, cs%PFv, g, gv, &
                       cs%PressureForce_CSp, cs%ALE_CSp, &
                       p_surf, cs%pbce, cs%eta_PF)
    call cpu_clock_end(id_clock_pres)
    call cpu_clock_begin(id_clock_pass)
    call pass_var(cs%eta_PF, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  if (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    call pass_vector_complete(pid_u_av, u_av, v_av, g%Domain)
    call pass_vector_complete(pid_uh, uh(:,:,:), vh(:,:,:), g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, tv, h, cs%update_OBC_CSp, time_local)
  endif; endif

  if (bt_cont_bt_thick) then
    call cpu_clock_begin(id_clock_pass)
    call pass_vector(cs%BT_cont%h_u, cs%BT_cont%h_v, g%Domain, &
                     to_all+scalar_pair, cgrid_ne)
    call cpu_clock_end(id_clock_pass)
    call legacy_btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v)
  endif

  if (cs%debug) then
    call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv)
    call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1)
    call hchksum(h_av,"Predictor avg h",g%HI,haloshift=0, scale=gv%H_to_m)
  ! call MOM_state_chksum("Predictor avg ", u_av, v_av,  h_av,uh, vh, G, GV)
    call check_redundant("Predictor up ", up, vp, g)
    call check_redundant("Predictor uh ", uh, vh, g)
  endif

! diffu = horizontal viscosity terms (u_av)
  call cpu_clock_begin(id_clock_horvisc)
  call horizontal_viscosity(u_av, v_av, h_av, cs%diffu, cs%diffv, &
                            meke, varmix, g, gv, cs%hor_visc_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_horvisc)

! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, g, gv, &
                 cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)

! Calculate the momentum forcing terms for the barotropic equations.

! u_bc_accel = CAu + PFu + diffu(u[n-1])
  call cpu_clock_begin(id_clock_btforce)
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
  enddo ; enddo ; enddo
  if (associated(cs%OBC)) then
    call open_boundary_zero_normal_flow(cs%OBC, g, u_bc_accel, v_bc_accel)
  endif
  call cpu_clock_end(id_clock_btforce)

  if (cs%debug) then
    call check_redundant("corr pre-btstep CS%Ca ", cs%Cau, cs%Cav, g)
    call check_redundant("corr pre-btstep CS%PF ", cs%PFu, cs%PFv, g)
    call check_redundant("corr pre-btstep CS%diff ", cs%diffu, cs%diffv, g)
    call check_redundant("corr pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g)
  endif

! u_accel_bt = layer accelerations due to barotropic solver
! pbce = dM/deta
  call cpu_clock_begin(id_clock_btstep)
  if (cs%flux_BT_coupling) then
    call legacy_btstep(.true., uh_in, vh_in, eta, dt, u_bc_accel, v_bc_accel, &
                fluxes, cs%pbce, cs%eta_PF, uh, vh, cs%u_accel_bt, &
                cs%v_accel_bt, eta, cs%uhbt, cs%vhbt, g, gv, &
                cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, etaav=eta_av, &
                uhbt_out = uhbt_out, vhbt_out = vhbt_out, obc=cs%OBC, &
                bt_cont = cs%BT_cont, eta_pf_start = eta_pf_start, &
                taux_bot=taux_bot, tauy_bot=tauy_bot)
  else
    if (cs%BT_use_layer_fluxes) then
      uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
    endif

    call legacy_btstep(.false., u, v, eta, dt, u_bc_accel, v_bc_accel, &
                fluxes, cs%pbce, cs%eta_PF, u_av, v_av, cs%u_accel_bt, &
                cs%v_accel_bt, eta, cs%uhbt, cs%vhbt, g, gv, &
                cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, &
                etaav=eta_av, obc=cs%OBC, &
                bt_cont = cs%BT_cont, eta_pf_start=eta_pf_start, &
                taux_bot=taux_bot, tauy_bot=tauy_bot, &
                uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr)
  endif
  call cpu_clock_end(id_clock_btstep)

  if (cs%debug) then
    call check_redundant("u_accel_bt ", cs%u_accel_bt, cs%v_accel_bt, g)
  endif

! u = u + dt*( u_bc_accel + u_accel_bt )
  call cpu_clock_begin(id_clock_mom_update)
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u(i,j,k) = g%mask2dCu(i,j) * (u_init(i,j,k) + dt * &
                    (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
  enddo ; enddo ; enddo

  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v(i,j,k) = g%mask2dCv(i,j) * (v_init(i,j,k) + dt * &
                    (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
  enddo ; enddo ; enddo
  call cpu_clock_end(id_clock_mom_update)

  if (cs%debug) then
    call uvchksum("Corrector 1 [uv]", u, v, g%HI, haloshift=0)
    call hchksum(h,"Corrector 1 h",g%HI,haloshift=2, scale=gv%H_to_m)
    call uvchksum("Corrector 1 [uv]h", &
                  uh, vh, g%HI, haloshift=2, scale=gv%H_to_m)
  ! call MOM_state_chksum("Corrector 1", u, v, h, uh, vh, G, GV, haloshift=1)
    call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
             cs%diffu, cs%diffv, g, gv, cs%pbce, cs%u_accel_bt, cs%v_accel_bt)
  endif

! u <- u + dt d/dz visc d/dz u
! u_av <- u_av + dt d/dz visc d/dz u_av
  call cpu_clock_begin(id_clock_vertvisc)
  call vertvisc_coef(u, v, h, fluxes, visc, dt, g, gv, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(u, v, h, fluxes, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, &
                cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
  if (g%nonblocking_updates) then
    call cpu_clock_end(id_clock_vertvisc) ; call cpu_clock_begin(id_clock_pass)
    pid_u = pass_vector_start(u, v, g%Domain)
    call cpu_clock_end(id_clock_pass) ; call cpu_clock_begin(id_clock_vertvisc)
  endif
  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)

! Later, h_av = (h_in + h_out)/2, but for now use h_av to store h_in.
  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
    h_av(i,j,k) = h(i,j,k)
  enddo ; enddo ; enddo

  call cpu_clock_begin(id_clock_pass)
  call pass_vector(cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                   to_all+scalar_pair, cgrid_ne)
  if (g%nonblocking_updates) then
    call pass_vector_complete(pid_u, u, v, g%Domain)
  else
    call pass_vector(u, v, g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

! uh = u_av * h
! h = h + dt * div . uh
  if (cs%flux_BT_coupling) then
    ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
    ! Also, determine the values of u and v so that their transports
    ! that agree with uhbt_out and vhbt_out.
    if (ASSOCIATED(cs%ADp%du_other)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq
      u_tmp(i,j,k) = u(i,j,k)
    enddo ; enddo ; enddo ; endif
    if (ASSOCIATED(cs%ADp%dv_other)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie
      v_tmp(i,j,k) = v(i,j,k)
    enddo ; enddo ; enddo ; endif
    call cpu_clock_begin(id_clock_continuity)
    call continuity(u, v, h, h, uh, vh, dt, g, gv, &
                    cs%continuity_CSp, cs%uhbt, cs%vhbt, cs%OBC, &
                    cs%visc_rem_u, cs%visc_rem_v, u_av, v_av, &
                    uhbt_out, vhbt_out, u, v)
    call cpu_clock_end(id_clock_continuity)
    ! Whenever thickness changes let the diag manager know, target grids
    ! for vertical remapping may need to be regenerated.
    call diag_update_remap_grids(cs%diag)
    if (g%nonblocking_updates) then
      call cpu_clock_begin(id_clock_pass)
      pid_h = pass_var_start(h, g%Domain)
      call cpu_clock_end(id_clock_pass)
    endif
    if (ASSOCIATED(cs%ADp%du_other)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq
      cs%ADp%du_other(i,j,k) = cs%ADp%du_other(i,j,k) + (u(i,j,k) - u_tmp(i,j,k))
    enddo ; enddo ; enddo ; endif
    if (ASSOCIATED(cs%ADp%dv_other)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie
      cs%ADp%dv_other(i,j,k) = cs%ADp%dv_other(i,j,k) + (v(i,j,k) - v_tmp(i,j,k))
    enddo ; enddo ; enddo ; endif

    call cpu_clock_begin(id_clock_vertvisc)
    call vertvisc_limit_vel(u, v, h_av, cs%ADp, cs%CDp, fluxes, visc, dt, g, gv, cs%vertvisc_CSp)
    if (g%nonblocking_updates) then
      call cpu_clock_end(id_clock_vertvisc) ; call cpu_clock_begin(id_clock_pass)
      pid_u = pass_vector_start(u, v, g%Domain)
      call cpu_clock_end(id_clock_pass) ; call cpu_clock_begin(id_clock_vertvisc)
    endif
    call vertvisc_limit_vel(u_av, v_av, h_av, cs%ADp, cs%CDp, fluxes, visc, dt, g, gv, cs%vertvisc_CSp)
    call cpu_clock_end(id_clock_vertvisc)

    call cpu_clock_begin(id_clock_pass)
    if (g%nonblocking_updates) then
      call pass_var_complete(pid_h, h, g%Domain)
      call pass_vector_complete(pid_u, u, v, g%Domain)
    else
      call pass_var(h, g%Domain)
      call pass_vector(u, v, g%Domain, complete=.false.)
    endif
    call cpu_clock_end(id_clock_pass)
  else
    ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
    call cpu_clock_begin(id_clock_continuity)
    call continuity(u, v, h, h, uh, vh, dt, g, gv, &
                    cs%continuity_CSp, cs%uhbt, cs%vhbt, cs%OBC, &
                    cs%visc_rem_u, cs%visc_rem_v, u_av, v_av)
    call cpu_clock_end(id_clock_continuity)
    ! Whenever thickness changes let the diag manager know, target grids
    ! for vertical remapping may need to be regenerated.
    call diag_update_remap_grids(cs%diag)
    call cpu_clock_begin(id_clock_pass)
    call pass_var(h, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  call cpu_clock_begin(id_clock_pass)
  if (g%nonblocking_updates) then
    pid_uh = pass_vector_start(uh(:,:,:), vh(:,:,:), g%Domain)
    pid_u_av = pass_vector_start(u_av, v_av, g%Domain)
  else
    call pass_vector(u_av, v_av, g%Domain, complete=.false.)
    call pass_vector(uh(:,:,:), vh(:,:,:), g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

  if (associated(cs%OBC)) then
    call radiation_open_bdry_conds(cs%OBC, u, u_old_rad_obc, v, v_old_rad_obc, g, dt)
  endif

! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in.
  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
    h_av(i,j,k) = 0.5*(h_av(i,j,k) + h(i,j,k))
  enddo ; enddo ; enddo

  if (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    call pass_vector_complete(pid_uh, uh(:,:,:), vh(:,:,:), g%Domain)
    call pass_vector_complete(pid_u_av, u_av, v_av, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  do k=1,nz ; do j=js-2,je+2 ; do i=isq-2,ieq+2
    uhtr(i,j,k) = uhtr(i,j,k) + uh(i,j,k)*dt
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq-2,jeq+2 ; do i=is-2,ie+2
    vhtr(i,j,k) = vhtr(i,j,k) + vh(i,j,k)*dt
  enddo ; enddo ; enddo

!   The time-averaged free surface height has already been set by the last
!  call to btstep.

!   Here various terms used in to update the momentum equations are
! offered for averaging.
  if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)
  if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)
  if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)
  if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)

!   Here the thickness fluxes are offered for averaging.
  if (cs%id_uh > 0) call post_data(cs%id_uh, uh, cs%diag)
  if (cs%id_vh > 0) call post_data(cs%id_vh, vh, cs%diag)
  if (cs%id_uav > 0) call post_data(cs%id_uav, u_av, cs%diag)
  if (cs%id_vav > 0) call post_data(cs%id_vav, v_av, cs%diag)
  if (cs%id_u_BT_accel > 0) call post_data(cs%id_u_BT_accel, cs%u_accel_bt, cs%diag)
  if (cs%id_v_BT_accel > 0) call post_data(cs%id_v_BT_accel, cs%v_accel_bt, cs%diag)
  if (cs%id_du_adj > 0) call post_data(cs%id_du_adj, cs%ADp%du_other, cs%diag)
  if (cs%id_dv_adj > 0) call post_data(cs%id_dv_adj, cs%ADp%dv_other, cs%diag)
  if (cs%debug) then
    call mom_state_chksum("Corrector ", u, v, h, uh, vh, g, gv)
    call uvchksum("Corrector avg [uv]", u_av, v_av, g%HI, haloshift=1)
    call hchksum(h_av,"Corrector avg h",g%HI,haloshift=1, scale=gv%H_to_m)
 !  call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV)
  endif

end subroutine step_mom_dyn_legacy_split

! =============================================================================

subroutine adjustments_dyn_legacy_split(u, v, h, dt, G, GV, CS)
  type(ocean_grid_type),          intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),        intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                  intent(inout) :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                  intent(inout) :: v    !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                                  intent(inout) :: h    !< Layer thicknesses, in H
                                                        !! (usually m or kg m-2).
  real,                           intent(in)    :: dt   !< The baroclinic dynamics time step, in s.
  type(mom_dyn_legacy_split_cs),  pointer       :: CS   !< The control structure set up by
                                                        !! initialize_dyn_legacy_split.

! Arguments: u - The zonal velocity, in m s-1.
!  (in)      v - The meridional velocity, in m s-1.
!  (in)      h - The layer thicknesses, in m or kg m-2, depending on
!                whether the Boussinesq approximation is made.
!  (in)      dt - The time step in s.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure set up by initialize_dyn_legacy_split.

  ! Temporary arrays to contain layer thickness fluxes in m3 s-1 or kg s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uh_temp
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vh_temp
  ! A temporary array to contain layer projected thicknesses in m or kg m-2.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: h_temp
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  if (cs%readjust_BT_trans) then
    call cpu_clock_begin(id_clock_continuity)
    call continuity(u, v, h, h_temp, uh_temp, vh_temp, dt, g, gv, &
                    cs%continuity_CSp, obc=cs%OBC)
    call cpu_clock_end(id_clock_continuity)
!$OMP parallel default(none) shared(is,ie,js,je,nz,CS,uh_temp,vh_temp)
!$OMP do
    do j=js,je
      do i=is-1,ie ; cs%uhbt_in(i,j) = uh_temp(i,j,1) ; enddo
      do k=2,nz ; do i=is-1,ie
        cs%uhbt_in(i,j) = cs%uhbt_in(i,j) + uh_temp(i,j,k)
      enddo ; enddo
    enddo
!$OMP do
    do j=js-1,je
      do i=is,ie ; cs%vhbt_in(i,j) = vh_temp(i,j,1) ; enddo
      do k=2,nz ; do i=is,ie
        cs%vhbt_in(i,j) = cs%vhbt_in(i,j) + vh_temp(i,j,k)
      enddo ; enddo
    enddo
!$OMP end parallel
    cs%readjust_velocity = .true.
  endif

end subroutine adjustments_dyn_legacy_split

! =============================================================================

subroutine register_restarts_dyn_legacy_split(HI, GV, param_file, CS, restart_CS, uh, vh)
  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(mom_dyn_legacy_split_cs), pointer      :: CS         !< The control structure set up by
                                                            !! initialize_dyn_legacy_split.
  type(mom_restart_cs),          pointer      :: restart_CS !< A pointer to the restart
                                                            !! control structure.
  real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
                        target, intent(inout) :: uh         !< The zonal volume or mass transport,
                                                            !! in m3 s-1 or kg s-1.
  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
                        target, intent(inout) :: vh         !< The meridional volume or mass
                                                            !! transport, in m3 s-1 or kg s-1.
!   This subroutine sets up any auxiliary restart variables that are specific
! to the unsplit time stepping scheme.  All variables registered here should
! have the ability to be recreated if they are not present in a restart file.

! Arguments: HI - A horizontal index type structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (inout)   CS - The control structure set up by initialize_dyn_legacy_split.
!  (in)      restart_CS - A pointer to the restart control structure.
!  (inout)   uh - The zonal volume or mass transport, in m3 s-1 or kg s-1.
!  (inout)   vh - The meridional volume or mass transport, in m3 s-1 or kg s-1.

  type(vardesc) :: vd
  character(len=40)  :: mdl = "MOM_dynamics_legacy_split" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  logical :: adiabatic, flux_BT_coupling, readjust_BT_trans
  integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB
  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
  isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB

! This is where a control structure that is specific to this module would be allocated.
  if (associated(cs)) then
    call mom_error(warning, "register_restarts_dyn_legacy_split called with an associated "// &
                             "control structure.")
    return
  endif
  allocate(cs)

  call get_param(param_file, "MOM_legacy_split", "FLUX_BT_COUPLING", flux_bt_coupling, &
       default=.false., do_not_log=.true.)
  call get_param(param_file, "MOM_legacy_split", "READJUST_BT_TRANS", readjust_bt_trans, &
       default=.false., do_not_log=.true.)
  call get_param(param_file, "MOM_legacy_split", "ADIABATIC", adiabatic, &
       default=.false., do_not_log=.true.)
  if (.not.flux_bt_coupling .or. adiabatic) readjust_bt_trans = .false.

  alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
  alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
  alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
  alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
  alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
  alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0

  alloc_(cs%eta(isd:ied,jsd:jed))       ; cs%eta(:,:) = 0.0
  alloc_(cs%u_av(isdb:iedb,jsd:jed,nz)) ; cs%u_av(:,:,:) = 0.0
  alloc_(cs%v_av(isd:ied,jsdb:jedb,nz)) ; cs%v_av(:,:,:) = 0.0
  alloc_(cs%h_av(isd:ied,jsd:jed,nz))   ; cs%h_av(:,:,:) = gv%Angstrom
  alloc_(cs%uhbt_in(isdb:iedb,jsd:jed)) ; cs%uhbt_in(:,:) = 0.0
  alloc_(cs%vhbt_in(isd:ied,jsdb:jedb)) ; cs%vhbt_in(:,:) = 0.0

  thickness_units = get_thickness_units(gv)
  flux_units = get_flux_units(gv)

  vd = var_desc("sfc",thickness_units,"Free surface Height",'h','1')
  call register_restart_field(cs%eta, vd, .false., restart_cs)

  vd = var_desc("u2","meter second-1","Auxiliary Zonal velocity",'u','L')
  call register_restart_field(cs%u_av, vd, .false., restart_cs)

  vd = var_desc("v2","meter second-1","Auxiliary Meridional velocity",'v','L')
  call register_restart_field(cs%v_av, vd, .false., restart_cs)

  vd = var_desc("h2",thickness_units,"Auxiliary Layer Thickness",'h','L')
  call register_restart_field(cs%h_av, vd, .false., restart_cs)

  vd = var_desc("uh",flux_units,"Zonal thickness flux",'u','L')
  call register_restart_field(uh, vd, .false., restart_cs)

  vd = var_desc("vh",flux_units,"Meridional thickness flux",'v','L')
  call register_restart_field(vh, vd, .false., restart_cs)

  vd = var_desc("diffu","meter second-2","Zonal horizontal viscous acceleration",'u','L')
  call register_restart_field(cs%diffu, vd, .false., restart_cs)

  vd = var_desc("diffv","meter second-2","Meridional horizontal viscous acceleration",'v','L')
  call register_restart_field(cs%diffv, vd, .false., restart_cs)

  call register_legacy_barotropic_restarts(hi, gv, param_file, &
           cs%barotropic_CSp, restart_cs)

  if (readjust_bt_trans) then
    vd = var_desc("uhbt_in",flux_units,"Final instantaneous barotropic zonal thickness flux",&
                 'u','1')
    call register_restart_field(cs%uhbt_in, vd, .false., restart_cs)

    vd = var_desc("vhbt_in",flux_units,"Final instantaneous barotropic meridional thickness flux",&
                 'v','1')
    call register_restart_field(cs%vhbt_in, vd, .false., restart_cs)
  endif

end subroutine register_restarts_dyn_legacy_split

subroutine initialize_dyn_legacy_split(u, v, h, uh, vh, eta, Time, G, GV, param_file, &
                      diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
                      VarMix, MEKE, OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, visc, &
                      dirs, ntrunc)
  type(ocean_grid_type),         intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),       intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                 intent(inout) :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                 intent(inout) :: v    !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , &
                                 intent(inout) :: h    !< Layer thicknesses, in H
                                                       !! (usually m or kg m-2).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                         target, intent(inout) :: uh   !< The zonal volume or mass transport,
                                                       !! in m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                         target, intent(inout) :: vh   !< The meridional volume or mass transport,
                                                       !! in m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJ_(G)),          &
                                 intent(inout) :: eta  !< The free surface height or column mass,
                                                       !! in m or kg m-2.
  type(time_type),       target, intent(in)    :: Time        !< The current model time.
  type(param_file_type),         intent(in)    :: param_file  !< A structure to parse for run-time
                                                              !! parameters.
  type(diag_ctrl),       target, intent(inout) :: diag        !< A structure that is used to
                                                              !! regulate diagnostic output.
  type(mom_dyn_legacy_split_cs), pointer       :: CS          !< The control structure set up by
                                                              !! initialize_dyn_legacy_split.
  type(mom_restart_cs),          pointer       :: restart_CS  !< A pointer to the restart control
                                                              !! structure.
  real,                          intent(in)    :: dt          !< The baroclinic dynamics time step,
                                                              !! in s.
  type(accel_diag_ptrs), target, intent(inout) :: Accel_diag  !<A set of pointers to the various
                                        !! accelerations in the momentum equations, which can be
                                        !! used for later derived diagnostics, like energy budgets.
  type(cont_diag_ptrs),  target, intent(inout) :: Cont_diag  !< A structure with pointers to various
                                                              !! terms in the continuity equations.
  type(ocean_internal_state),    intent(inout) :: MIS         !< The "MOM6 Internal State"
                                              !! structure, used to pass around pointers to
                                              !! various arrays for diagnostic purposes.
  type(varmix_cs),               pointer       :: VarMix      !< A pointer to a structure with
                                     !! fields that specify the spatially variable viscosities.
  type(meke_type),               pointer       :: MEKE        !< A pointer to a structure containing
                                              !! fields related to the Mesoscale Eddy Kinetic Energy.
  type(ocean_obc_type),          pointer       :: OBC         !< If open boundary conditions are
                                                          !! used, this points to the ocean_OBC_type
                                                          !! that was set up in MOM_initialization.
  type(update_obc_cs),           pointer       :: update_OBC_CSp !< If open boundary condition
                                                              !! updates are used, this points to
                                                              !! the appropriate control structure.
  type(ale_cs),                  pointer       :: ALE_CSp     !< This points to the ALE control
                                                              !! structure.
  type(set_visc_cs),             pointer       :: setVisc_CSp !< This points to the set_visc control
                                                              !! structure.
  type(vertvisc_type),           intent(inout) :: visc        !< A structure containing vertical
                                         !! viscosities, bottom drag viscosities, and related fields.
  type(directories),             intent(in)    :: dirs        !< A structure containing several
                                                              !! relevant directory paths.
  integer, target,               intent(inout) :: ntrunc      !< A target for the variable that
                        !! records the number of times the velocity is truncated (this should be 0).

! Arguments: u - The zonal velocity, in m s-1.
!  (inout)   v - The meridional velocity, in m s-1.
!  (inout)   h - The layer thicknesses, in m or kg m-2, depending on whether
!                the Boussinesq approximation is made.
!  (inout)   uh - The zonal volume or mass transport, in m3 s-1 or kg s-1.
!  (inout)   vh - The meridional volume or mass transport, in m3 s-1 or kg s-1.
!  (out)     eta - The free surface height or column mass, in m or kg m-2.
!  (in)      Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (inout)   CS - The control structure set up by initialize_dyn_legacy_split.
!  (in)      restart_CS - A pointer to the restart control structure.
!  (inout)   Accel_diag - A set of pointers to the various accelerations in
!                  the momentum equations, which can be used for later derived
!                  diagnostics, like energy budgets.
!  (inout)   Cont_diag - A structure with pointers to various terms in the
!                   continuity equations.
!  (inout)   MIS - The "MOM6 Internal State" structure, used to pass around
!                  pointers to various arrays for diagnostic purposes.
!  (in)      VarMix - A pointer to a structure with fields that specify the
!                     spatially variable viscosities.
!  (inout)   MEKE - A pointer to a structure containing fields related to
!                   the Mesoscale Eddy Kinetic Energy.
!  (in)      OBC - If open boundary conditions are used, this points to the
!                  ocean_OBC_type that was set up in MOM_initialization.
!  (in)      update_OBC_CSp - If open boundary condition updates are used,
!                  this points to the appropriate control structure.
!  (in)      ALE_CS - This points to the ALE control structure.
!  (in)      setVisc_CSp - This points to the set_visc control structure.
!  (inout)   visc - A structure containing vertical viscosities, bottom drag
!                   viscosities, and related fields.
!  (in)      dirs - A structure containing several relevant directory paths.
!  (in)      ntrunc - A target for the variable that records the number of times
!                     the velocity is truncated (this should be 0).

  !   This subroutine initializes all of the variables that are used by this
  ! dynamic core, including diagnostics and the cpu clocks.

  !   This subroutine initializes any variables that are specific to this time
  ! stepping scheme, including the cpu clocks.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp
  character(len=40)  :: mdl = "MOM_dynamics_legacy_split" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  logical :: adiabatic, use_tides, debug_truncations
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
  integer :: IsdB, IedB, JsdB, JedB
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (.not.associated(cs)) call mom_error(fatal, &
      "initialize_dyn_legacy_split called with an unassociated control structure.")
  if (cs%module_is_initialized) then
    call mom_error(warning, "initialize_dyn_legacy_split called with a control "// &
                            "structure that has already been initialized.")
    return
  endif
  cs%module_is_initialized = .true.

  cs%diag => diag

  call get_param(param_file, mdl, "TIDES", use_tides, &
                 "If true, apply tidal momentum forcing.", default=.false.)
  call get_param(param_file, mdl, "BE", cs%be, &
                 "If SPLIT is true, BE determines the relative weighting \n"//&
                 "of a  2nd-order Runga-Kutta baroclinic time stepping \n"//&
                 "scheme (0.5) and a backward Euler scheme (1) that is \n"//&
                 "used for the Coriolis and inertial terms.  BE may be \n"//&
                 "from 0.5 to 1, but instability may occur near 0.5. \n"//&
                 "BE is also applicable if SPLIT is false and USE_RK2 \n"//&
                 "is true.", units="nondim", default=0.6)
  call get_param(param_file, mdl, "BEGW", cs%begw, &
                 "If SPLIT is true, BEGW is a number from 0 to 1 that \n"//&
                 "controls the extent to which the treatment of gravity \n"//&
                 "waves is forward-backward (0) or simulated backward \n"//&
                 "Euler (1).  0 is almost always used.\n"//&
                 "If SPLIT is false and USE_RK2 is true, BEGW can be \n"//&
                 "between 0 and 0.5 to damp gravity waves.", &
                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "FLUX_BT_COUPLING", cs%flux_BT_coupling, &
                 "If true, use mass fluxes to ensure consistency between \n"//&
                 "the baroclinic and barotropic modes. This is only used \n"//&
                 "if SPLIT is true.", default=.false.)
  call get_param(param_file, mdl, "READJUST_BT_TRANS", cs%readjust_BT_trans, &
                 "If true, make a barotropic adjustment to the layer \n"//&
                 "velocities after the thermodynamic part of the step \n"//&
                 "to ensure that the interaction between the thermodynamics \n"//&
                 "and the continuity solver do not change the barotropic \n"//&
                 "transport.  This is only used if FLUX_BT_COUPLING and \n"//&
                 "SPLIT are true.", default=.false.)
  call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
                 "If true, provide the bottom stress calculated by the \n"//&
                 "vertical viscosity to the barotropic solver.", default=.false.)
  call get_param(param_file, mdl, "BT_USE_LAYER_FLUXES", cs%BT_use_layer_fluxes, &
                 "If true, use the summed layered fluxes plus an \n"//&
                 "adjustment due to the change in the barotropic velocity \n"//&
                 "in the barotropic continuity equation.", default=.true.)
  call get_param(param_file, mdl, "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", default=.false.)
  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., do_not_log=.true.)
  if (.not.cs%flux_BT_coupling .or. adiabatic) cs%readjust_BT_trans = .false.
  call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
                 default=.false.)

  allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
  allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0

  alloc_(cs%uhbt(isdb:iedb,jsd:jed))   ; cs%uhbt(:,:) = 0.0
  alloc_(cs%vhbt(isd:ied,jsdb:jedb))   ; cs%vhbt(:,:) = 0.0
  alloc_(cs%visc_rem_u(isdb:iedb,jsd:jed,nz)) ; cs%visc_rem_u(:,:,:) = 0.0
  alloc_(cs%visc_rem_v(isd:ied,jsdb:jedb,nz)) ; cs%visc_rem_v(:,:,:) = 0.0
  alloc_(cs%eta_PF(isd:ied,jsd:jed))   ; cs%eta_PF(:,:) = 0.0
  alloc_(cs%pbce(isd:ied,jsd:jed,nz))  ; cs%pbce(:,:,:) = 0.0

  alloc_(cs%u_accel_bt(isdb:iedb,jsd:jed,nz)) ; cs%u_accel_bt(:,:,:) = 0.0
  alloc_(cs%v_accel_bt(isd:ied,jsdb:jedb,nz)) ; cs%v_accel_bt(:,:,:) = 0.0

  mis%diffu => cs%diffu ; mis%diffv => cs%diffv
  mis%PFu => cs%PFu ; mis%PFv => cs%PFv
  mis%CAu => cs%CAu ; mis%CAv => cs%CAv
  mis%pbce => cs%pbce
  mis%u_accel_bt => cs%u_accel_bt ; mis%v_accel_bt => cs%v_accel_bt
  mis%u_av => cs%u_av ; mis%v_av => cs%v_av

  cs%ADp => accel_diag ; cs%CDp => cont_diag
  accel_diag%diffu => cs%diffu ; accel_diag%diffv => cs%diffv
  accel_diag%PFu => cs%PFu ; accel_diag%PFv => cs%PFv
  accel_diag%CAu => cs%CAu ; accel_diag%CAv => cs%CAv
!  Accel_diag%pbce => CS%pbce
!  Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
!  Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av

  call continuity_init(time, g, gv, param_file, diag, cs%continuity_CSp)
  call coriolisadv_init(time, g,  param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
  if (use_tides) call tidal_forcing_init(time, g, param_file, cs%tides_CSp)
  call pressureforce_init(time, g, gv, param_file, diag, cs%PressureForce_CSp, &
                          cs%tides_CSp)
  call hor_visc_init(time, g, param_file, diag, cs%hor_visc_CSp)
  call vertvisc_init(mis, time, g, gv, param_file, diag, cs%ADp, dirs, &
                     ntrunc, cs%vertvisc_CSp)
  if (.not.associated(setvisc_csp)) call mom_error(fatal, &
    "initialize_dyn_legacy_split called with setVisc_CSp unassociated.")
  cs%set_visc_CSp => setvisc_csp

  if (associated(ale_csp)) cs%ALE_CSp => ale_csp
  if (associated(obc)) cs%OBC => obc
  if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp

  if (.not. query_initialized(cs%eta,"sfc",restart_cs))  then
    ! Estimate eta based on the layer thicknesses - h.  With the Boussinesq
    ! approximation, eta is the free surface height anomaly, while without it
    ! eta is the mass of ocean per unit area.  eta always has the same
    ! dimensions as h, either m or kg m-3.
    !   CS%eta(:,:) = 0.0 already from initialization.
    if (gv%Boussinesq) then
      do j=js,je ; do i=is,ie ; cs%eta(i,j) = -g%bathyT(i,j) ; enddo ; enddo
    endif
    do k=1,nz ; do j=js,je ; do i=is,ie
       cs%eta(i,j) = cs%eta(i,j) + h(i,j,k)
    enddo ; enddo ; enddo
  endif
  ! Copy eta into an output array.
  do j=js,je ; do i=is,ie ; eta(i,j) = cs%eta(i,j) ; enddo ; enddo

  call legacy_barotropic_init(u, v, h, cs%eta, time, g, gv, param_file, diag, &
                       cs%barotropic_CSp, restart_cs, cs%BT_cont, cs%tides_CSp)

  if (.not. query_initialized(cs%diffu,"diffu",restart_cs) .or. &
      .not. query_initialized(cs%diffv,"diffv",restart_cs)) &
    call horizontal_viscosity(u, v, h, cs%diffu, cs%diffv, meke, varmix, &
                              g, gv, cs%hor_visc_CSp)
  if (.not. query_initialized(cs%u_av,"u2", restart_cs) .or. &
      .not. query_initialized(cs%u_av,"v2", restart_cs)) then
    cs%u_av(:,:,:) = u(:,:,:)
    cs%v_av(:,:,:) = v(:,:,:)
  endif
! This call is just here to initialize uh and vh.
  if (.not. query_initialized(uh,"uh",restart_cs) .or. &
      .not. query_initialized(vh,"vh",restart_cs)) then
    h_tmp(:,:,:) = h(:,:,:)
    call continuity(u, v, h, h_tmp, uh, vh, dt, g, gv, cs%continuity_CSp, obc=cs%OBC)
    call cpu_clock_begin(id_clock_pass_init)
    call pass_var(h_tmp, g%Domain)
    call cpu_clock_end(id_clock_pass_init)
    cs%h_av(:,:,:) = 0.5*(h(:,:,:) + h_tmp(:,:,:))
  else
    if (.not. query_initialized(cs%h_av,"h2",restart_cs)) &
      cs%h_av(:,:,:) = h(:,:,:)
  endif

  !   Determine whether there is a barotropic transport that is to be used
  ! to adjust the layers' velocities.
  cs%readjust_velocity = .false.
  if (cs%readjust_BT_trans) then
    if (query_initialized(cs%uhbt_in,"uhbt_in",restart_cs) .and. &
        query_initialized(cs%vhbt_in,"vhbt_in",restart_cs)) then
      cs%readjust_velocity = .true.
      call cpu_clock_begin(id_clock_pass_init)
      call pass_vector(cs%uhbt_in, cs%vhbt_in, g%Domain)
      call cpu_clock_end(id_clock_pass_init)
    endif
  endif

  call cpu_clock_begin(id_clock_pass_init)
  call pass_vector(cs%u_av,cs%v_av, g%Domain)
  call pass_var(cs%h_av, g%Domain)
  call pass_vector(uh, vh, g%Domain)
  call cpu_clock_end(id_clock_pass_init)

  flux_units = get_flux_units(gv)
  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
      'Zonal Thickness Flux', flux_units)
  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
      'Meridional Thickness Flux', flux_units)
  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
      'Zonal Coriolis and Advective Acceleration', 'meter second-2')
  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
      'Meridional Coriolis and Advective Acceleration', 'meter second-2')
  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
      'Zonal Pressure Force Acceleration', 'meter second-2')
  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
      'Meridional Pressure Force Acceleration', 'meter second-2')

  cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
      'Barotropic-step Averaged Zonal Velocity', 'meter second-1')
  cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
      'Barotropic-step Averaged Meridional Velocity', 'meter second-1')


  if (cs%flux_BT_coupling) then
    cs%id_du_adj = register_diag_field('ocean_model', 'du_adj', diag%axesCuL, time, &
        'Zonal velocity adjustments due to nonstandard terms', 'meter second-1')
    cs%id_dv_adj = register_diag_field('ocean_model', 'dv_adj', diag%axesCvL, time, &
        'Meridional velocity adjustments due to nonstandard terms', 'meter second-1')
    if (cs%id_du_adj > 0) call safe_alloc_ptr(cs%ADp%du_other,isdb,iedb,jsd,jed,nz)
    if (cs%id_dv_adj > 0) call safe_alloc_ptr(cs%ADp%dv_other,isd,ied,jsdb,jedb,nz)
  endif
  cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
    'Barotropic Anomaly Zonal Acceleration', 'meter second-1')
  cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
    'Barotropic Anomaly Meridional Acceleration', 'meter second-1')

  if (debug_truncations) then
    if (cs%flux_BT_coupling) then
      call safe_alloc_ptr(cs%ADp%du_other,isdb,iedb,jsd,jed,nz)
      call safe_alloc_ptr(cs%ADp%dv_other,isd,ied,jsdb,jedb,nz)
    endif
  endif

  id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
  id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
  id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
  id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
  id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
  id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
  id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
  id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', grain=clock_routine)
  id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=clock_module)
  id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=clock_module)
  id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=clock_module)


end subroutine initialize_dyn_legacy_split

subroutine end_dyn_legacy_split(CS)
  type(mom_dyn_legacy_split_cs), pointer :: CS
!  (inout)   CS - The control structure set up by initialize_dyn_legacy_split.

  dealloc_(cs%diffu) ; dealloc_(cs%diffv)
  dealloc_(cs%CAu)   ; dealloc_(cs%CAv)
  dealloc_(cs%PFu)   ; dealloc_(cs%PFv)

  if (associated(cs%taux_bot)) deallocate(cs%taux_bot)
  if (associated(cs%tauy_bot)) deallocate(cs%tauy_bot)
  dealloc_(cs%uhbt) ; dealloc_(cs%vhbt)
  dealloc_(cs%u_accel_bt) ; dealloc_(cs%v_accel_bt)
  dealloc_(cs%visc_rem_u) ; dealloc_(cs%visc_rem_v)

  dealloc_(cs%eta) ; dealloc_(cs%eta_PF) ; dealloc_(cs%pbce)
  dealloc_(cs%h_av) ; dealloc_(cs%u_av) ; dealloc_(cs%v_av)
  dealloc_(cs%uhbt_in) ; dealloc_(cs%vhbt_in)

  call dealloc_bt_cont_type(cs%BT_cont)

  deallocate(cs)
end subroutine end_dyn_legacy_split

end module mom_dynamics_legacy_split