MOM_dynamics_split_RK2.F90

Go to the documentation of this file.

!> Time step the adiabatic dynamic core of MOM using RK2 method.
module mom_dynamics_split_rk2

! This file is part of MOM6. See LICENSE.md for the license.

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
use mom_domains,           only : to_south, to_west, to_all, cgrid_ne, scalar_pair
use mom_domains,           only : to_north, to_east, omit_corners
use mom_domains,           only : create_group_pass, do_group_pass, group_pass_type
use mom_domains,           only : start_group_pass, complete_group_pass
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, calltree_showquery
use mom_error_handler,     only : calltree_enter, calltree_leave, calltree_waypoint
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_barotropic,            only : barotropic_init, btstep, btcalc, bt_mass_source
use mom_barotropic,            only : register_barotropic_restarts, set_dtbt, barotropic_cs
use mom_boundary_update,       only : update_obc_data, update_obc_cs
use mom_continuity,            only : continuity, continuity_init, continuity_cs
use mom_continuity,            only : continuity_stencil
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, radiation_open_bdry_conds
use mom_open_boundary,         only : open_boundary_zero_normal_flow
use mom_pressureforce,         only : pressureforce, pressureforce_init, pressureforce_cs
use mom_set_visc,              only : set_viscous_ml, set_visc_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_vert_friction,         only : updatecfltruncationvalue
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>

!> Module control structure
type, public :: mom_dyn_split_rk2_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    !< layer x-velocity with vertical mean replaced by
                                                                   !! time-mean barotropic velocity over a baroclinic timestep (m s-1)
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_av    !< layer y-velocity with vertical mean replaced by
                                                                   !! time-mean barotropic velocity over a baroclinic timestep (m s-1)
  real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_)      :: h_av    !< arithmetic mean of two successive layer thicknesses (m or kg m-2)
  real allocable_, dimension(NIMEM_,NJMEM_)             :: eta_pf  !< instantaneous SSH used in calculating PFu and PFv (meter)
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_)        :: uhbt    !< average x-volume or mass flux determined by barotropic solver
                                                                   !! (m3 s-1 or kg s-1). uhbt should (roughly?) equal to vertical sum of uh.
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_)        :: vhbt    !< average y-volume or mass flux determined by barotropic solver
                                                                   !! (m3 s-1 or kg s-1). vhbt should (roughly?) equal to vertical sum of vh.
  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()  !<  frictional x-bottom stress from the ocean to the seafloor (Pa)
  real, pointer, dimension(:,:) :: tauy_bot => null()  !<  frictional y-bottom stress from the ocean to the seafloor (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 :: 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 :: 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_umo    = -1, id_vmo    = -1
  integer :: id_umo_2d = -1, id_vmo_2d = -1
  integer :: id_pfu    = -1, id_pfv    = -1
  integer :: id_cau    = -1, id_cav    = -1

  ! Split scheme only.
  integer :: id_uav        = -1, id_vav        = -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.

  ! Remainder of the structure points 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(barotropic_cs),    pointer :: barotropic_csp    => null()
  type(vertvisc_cs),      pointer :: vertvisc_csp      => null()
  type(set_visc_cs),      pointer :: set_visc_csp      => null()
  type(tidal_forcing_cs), pointer :: tides_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()

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

  ! for group halo pass
  type(group_pass_type) :: pass_eta
  type(group_pass_type) :: pass_visc_rem, pass_uvp
  type(group_pass_type) :: pass_hp_uv
  type(group_pass_type) :: pass_uv
  type(group_pass_type) :: pass_h, pass_av_uvh

end type mom_dyn_split_rk2_cs


public step_mom_dyn_split_rk2
public register_restarts_dyn_split_rk2
public initialize_dyn_split_rk2
public end_dyn_split_rk2

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

!> RK2 splitting for time stepping MOM adiabatic dynamics
subroutine step_mom_dyn_split_rk2(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             !< ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: GV            !< ocean vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(inout) :: u     !< zonal velocity (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: v     !< merid velocity (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h             !< layer thickness (m or kg/m2)
  type(thermo_var_ptrs),                     intent(in)    :: tv            !< thermodynamic type
  type(vertvisc_type),                       intent(inout) :: visc          !< vertical visc, bottom drag, and related
  type(time_type),                           intent(in)    :: Time_local    !< model time at end of time step
  real,                                      intent(in)    :: dt            !< time step (sec)
  type(forcing),                             intent(in)    :: fluxes        !< forcing fields
  real, dimension(:,:),                      pointer       :: p_surf_begin  !< surf pressure at start of this dynamic time step (Pa)
  real, dimension(:,:),                      pointer       :: p_surf_end    !< surf pressure at end   of this dynamic time step (Pa)
  real,                                      intent(in)    :: dt_since_flux !< elapsed time since fluxes were applied (sec)
  real,                                      intent(in)    :: dt_therm      !< thermodynamic time step (sec)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(inout) :: uh    !< zonal volume/mass transport (m3/s or kg/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: vh    !< merid volume/mass transport (m3/s or kg/s)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr          !< accumulatated zonal volume/mass transport since last tracer advection (m3 or kg)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr          !< accumulatated merid volume/mass transport since last tracer advection (m3 or kg)
  real, dimension(SZI_(G),SZJ_(G)),          intent(out)   :: eta_av        !< free surface height or column mass time averaged over time step (m or kg/m2)
  type(mom_dyn_split_rk2_cs),                pointer       :: CS            !< module control structure
  logical,                                   intent(in)    :: calc_dtbt     !< if true, recalculate barotropic time step
  type(varmix_cs),                           pointer       :: VarMix        !< specify the spatially varying viscosities
  type(meke_type),                           pointer       :: MEKE          !< related to mesoscale eddy kinetic energy param

  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 velocity in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp   ! Predicted meridional velocity 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_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(), &
    eta => null()

  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.
  !---For group halo pass
  logical :: showCallTree, sym

  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: cont_stencil
  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
  idt = 1.0 / dt

  sym=.false.;if (g%Domain%symmetric) sym=.true.  ! switch to include symmetric domain in checksums

  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("step_MOM_dyn_split_RK2(), MOM_dynamics_split_RK2.F90")

  !$OMP parallel do default(shared)
  do k = 1, nz
    do j=g%jsd,g%jed   ; do i=g%isdB,g%iedB ;  up(i,j,k) = 0.0 ; enddo ; enddo
    do j=g%jsdB,g%jedB ; do i=g%isd,g%ied   ;  vp(i,j,k) = 0.0 ; enddo ; enddo
    do j=g%jsd,g%jed   ; do i=g%isd,g%ied   ;  hp(i,j,k) = h(i,j,k) ; enddo ; enddo
  enddo

  ! Update CFL truncation value as function of time
  call updatecfltruncationvalue(time_local, cs%vertvisc_CSp)

  if (cs%debug) then
    call mom_state_chksum("Start predictor ", u, v, h, uh, vh, g, gv, symmetric=sym)
    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-1,je+1 ; do i=is-2,ie+1
      u_old_rad_obc(i,j,k) = u_av(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=js-2,je+1 ; do i=is-1,ie+1
      v_old_rad_obc(i,j,k) = v_av(i,j,k)
    enddo ; enddo ; enddo
  endif

  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

  !--- begin set up for group halo pass

  call cpu_clock_begin(id_clock_pass)
  cont_stencil = continuity_stencil(cs%continuity_CSp)
  !### Apart from circle_OBCs halo for eta could be 1, but halo>=3 is required
  !### to match circle_OBCs solutions. Why?
  call create_group_pass(cs%pass_eta, eta, g%Domain) !### , halo=1)
  call create_group_pass(cs%pass_visc_rem, cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                         to_all+scalar_pair, cgrid_ne, halo=max(1,cont_stencil))
  call create_group_pass(cs%pass_uvp, up, vp, g%Domain, halo=max(1,cont_stencil))
  call create_group_pass(cs%pass_hp_uv, hp, g%Domain, halo=2)
  call create_group_pass(cs%pass_hp_uv, u_av, v_av, g%Domain, halo=2)
  call create_group_pass(cs%pass_hp_uv, uh(:,:,:), vh(:,:,:), g%Domain, halo=2)

  call create_group_pass(cs%pass_uv, u, v, g%Domain, halo=max(2,cont_stencil))
  call create_group_pass(cs%pass_h, h, g%Domain, halo=max(2,cont_stencil))
  call create_group_pass(cs%pass_av_uvh, u_av, v_av, g%Domain, halo=2)
  call create_group_pass(cs%pass_av_uvh, uh(:,:,:), vh(:,:,:), g%Domain, halo=2)
  call cpu_clock_end(id_clock_pass)
  !--- end set up for group halo pass


! 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
    pa_to_eta = 1.0 / gv%H_to_Pa
    !$OMP parallel do default(shared)
    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 (showcalltree) call calltree_waypoint("done with PressureForce (step_MOM_dyn_split_RK2)")

  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 (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    call start_group_pass(cs%pass_eta, g%Domain)
    call cpu_clock_end(id_clock_pass)
  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)
  if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")

! u_bc_accel = CAu + PFu + diffu(u[n-1])
  call cpu_clock_begin(id_clock_btforce)
  !$OMP parallel do default(shared)
  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
    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 mom_accel_chksum("pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, cs%pbce, u_bc_accel, v_bc_accel, &
                          symmetric=sym)
    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

  call cpu_clock_begin(id_clock_vertvisc)
  !$OMP parallel do default(shared)
  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
    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)

  if (cs%debug) then
    call uvchksum("before vertvisc: up", up, vp, g%HI, haloshift=0, symmetric=sym)
  endif
  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)
  if (showcalltree) call calltree_waypoint("done with vertvisc_coef (step_MOM_dyn_split_RK2)")


  call cpu_clock_begin(id_clock_pass)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_eta, g%Domain)
    call start_group_pass(cs%pass_visc_rem, g%Domain)
  else
    call do_group_pass(cs%pass_eta, g%Domain)
    call do_group_pass(cs%pass_visc_rem, g%Domain)
  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 btcalc(h, g, gv, cs%barotropic_CSp, obc=cs%OBC)
  call 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 complete_group_pass(cs%pass_visc_rem, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

! u_accel_bt = layer accelerations due to barotropic solver
  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 btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
                  obc=cs%OBC)
    endif
    if (showcalltree) call calltree_waypoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2)")
  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 set_dtbt(g, gv, cs%barotropic_CSp, eta, cs%pbce)
  if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
  ! This is the predictor step call to btstep.
  call btstep(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)
  if (showcalltree) call calltree_leave("btstep()")
  call cpu_clock_end(id_clock_btstep)

! up = u + dt_pred*( u_bc_accel + u_accel_bt )
  dt_pred = dt * cs%be
  call cpu_clock_begin(id_clock_mom_update)

  !$OMP parallel do default(shared)
  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
    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, symmetric=sym)
    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, &
                  symmetric=sym, 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, symmetric=sym)
    call mom_state_chksum("Predictor 1 init", u_init, v_init, h, uh, vh, g, gv, haloshift=2, &
                          symmetric=sym)
    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)
  if (cs%debug) then
    call uvchksum("0 before vertvisc: [uv]p", up, vp, g%HI,haloshift=0, symmetric=sym)
  endif
  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 (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
  if (g%nonblocking_updates) then
    call cpu_clock_end(id_clock_vertvisc) ; call cpu_clock_begin(id_clock_pass)
    call start_group_pass(cs%pass_uvp, 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 do_group_pass(cs%pass_visc_rem, g%Domain)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_uvp, g%Domain)
  else
    call do_group_pass(cs%pass_uvp, 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)
  if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")


  if (associated(cs%OBC)) then
    call cpu_clock_begin(id_clock_pass)
    ! These should be done with a pass that excludes uh & vh.
    call do_group_pass(cs%pass_hp_uv, g%Domain)
    call cpu_clock_end(id_clock_pass)

    call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, dt_pred)
  endif

  call cpu_clock_begin(id_clock_pass)
  call do_group_pass(cs%pass_hp_uv, g%Domain)
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_av_uvh, g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

  ! h_av = (h + hp)/2
  !$OMP parallel do default(shared)
  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 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.
    !$OMP parallel do default(shared)
    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)
    if (showcalltree) call calltree_waypoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2)")
  endif

  if (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    call complete_group_pass(cs%pass_av_uvh, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  if (bt_cont_bt_thick) then
    call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
                obc=cs%OBC)
    if (showcalltree) call calltree_waypoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2)")
  endif

  if (cs%debug) then
    call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv, symmetric=sym)
    call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym)
    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)
  if (showcalltree) call calltree_waypoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")

! 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)
  if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")

! 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)
  !$OMP parallel do default(shared)
  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
    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 mom_accel_chksum("corr pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, cs%pbce, u_bc_accel, v_bc_accel, &
                          symmetric=sym)
    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%BT_use_layer_fluxes) then
    uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
  endif

  if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
  ! This is the corrector step call to btstep.
  call btstep(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, &
              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)
  do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo
  call cpu_clock_end(id_clock_btstep)
  if (showcalltree) call calltree_leave("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)
  !$OMP parallel do default(shared)
  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
    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, symmetric=sym)
    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, &
                  symmetric=sym, 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, &
                          symmetric=sym)
  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)
    call start_group_pass(cs%pass_uv, 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)
  if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")

! Later, h_av = (h_in + h_out)/2, but for now use h_av to store h_in.
  !$OMP parallel do default(shared)
  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 do_group_pass(cs%pass_visc_rem, g%Domain)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_uv, g%Domain)
  else
    call do_group_pass(cs%pass_uv, g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

  ! uh = u_av * h
  ! h  = h + dt * div . uh
  ! 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)
  call cpu_clock_begin(id_clock_pass)
  call do_group_pass(cs%pass_h, g%Domain)
  call cpu_clock_end(id_clock_pass)
  ! 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 (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")

  call cpu_clock_begin(id_clock_pass)
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_av_uvh, g%Domain)
  else
    call do_group_pass(cs%pass_av_uvh, 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.
  !$OMP parallel do default(shared)
  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 complete_group_pass(cs%pass_av_uvh, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif
  !$OMP parallel do default(shared)
  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
    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 time 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 time 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%debug) then
    call mom_state_chksum("Corrector ", u, v, h, uh, vh, g, gv, symmetric=sym)
    call uvchksum("Corrector avg [uv]", u_av, v_av, g%HI,haloshift=1, symmetric=sym)
    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

  if (showcalltree) call calltree_leave("step_MOM_dyn_split_RK2()")

end subroutine step_mom_dyn_split_rk2

!> 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.
subroutine register_restarts_dyn_split_rk2(HI, GV, param_file, CS, restart_CS, uh, vh)
  type(hor_index_type),          intent(in)    :: HI         !< Horizontal index structure
  type(verticalgrid_type),       intent(in)    :: GV         !< ocean vertical grid structure
  type(param_file_type),         intent(in)    :: param_file !< parameter file
  type(mom_dyn_split_rk2_cs),    pointer       :: CS         !< module control structure
  type(mom_restart_cs),          pointer       :: restart_CS !< restart control structure
  real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), target, intent(inout) :: uh !< zonal volume/mass transport (m3/s or kg/s)
  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), target, intent(inout) :: vh !< merid volume/mass transport (m3/s or kg/s)

  type(vardesc)      :: vd
  character(len=40)  :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
  character(len=48)  :: thickness_units, flux_units

  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 specific to this module would be allocated.
  if (associated(cs)) then
    call mom_error(warning, "register_restarts_dyn_split_RK2 called with an associated "// &
                             "control structure.")
    return
  endif
  allocate(cs)

  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

  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_barotropic_restarts(hi, gv, param_file, cs%barotropic_CSp, &
                                    restart_cs)

end subroutine register_restarts_dyn_split_rk2

!> This subroutine initializes all of the variables that are used by this
!! dynamic core, including diagnostics and the cpu clocks.
subroutine initialize_dyn_split_rk2(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           !< ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: GV          !< ocean vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u           !< zonal velocity (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v           !< merid velocity (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h           !< layer thickness (m or kg/m2)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(inout) :: uh  !< zonal volume/mass transport (m3 s-1 or kg s-1)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: vh  !< merid volume/mass transport (m3 s-1 or kg s-1)
  real, dimension(SZI_(G),SZJ_(G)),          intent(inout) :: eta         !< free surface height or column mass (m or kg m-2)
  type(time_type),                   target, intent(in)    :: Time        !< current model time
  type(param_file_type),                     intent(in)    :: param_file  !< parameter file for parsing
  type(diag_ctrl),                   target, intent(inout) :: diag        !< to control diagnostics
  type(mom_dyn_split_rk2_cs),                pointer       :: CS          !< module control structure
  type(mom_restart_cs),                      pointer       :: restart_CS  !< restart control structure
  real,                                      intent(in)    :: dt          !< time step (sec)
  type(accel_diag_ptrs),             target, intent(inout) :: Accel_diag  !< points to momentum equation terms for budget analysis
  type(cont_diag_ptrs),              target, intent(inout) :: Cont_diag   !< points to terms in continuity equation
  type(ocean_internal_state),                intent(inout) :: MIS         !< "MOM6 internal state" used to pass diagnostic pointers
  type(varmix_cs),                           pointer       :: VarMix      !< points to spatially variable viscosities
  type(meke_type),                           pointer       :: MEKE        !< points to mesoscale eddy kinetic energy fields
  type(ocean_obc_type),                      pointer       :: OBC         !< points to OBC related fields
  type(update_obc_cs),                       pointer       :: update_OBC_CSp !< points to OBC update related fields
  type(ale_cs),                              pointer       :: ALE_CSp     !< points to ALE control structure
  type(set_visc_cs),                         pointer       :: setVisc_CSp !< points to the set_visc control structure.
  type(vertvisc_type),                       intent(inout) :: visc        !< vertical viscosities, bottom drag, and related
  type(directories),                         intent(in)    :: dirs        !< contains 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).

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp
  character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  type(group_pass_type) :: pass_h_tmp, pass_av_h_uvh
  logical :: 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_split_RK2 called with an unassociated control structure.")
  if (cs%module_is_initialized) then
    call mom_error(warning, "initialize_dyn_split_RK2 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, "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, "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_split_RK2 called with setVisc_CSp unassociated.")
  cs%set_visc_CSp => setvisc_csp
  call updatecfltruncationvalue(time, cs%vertvisc_CSp, activate= &
            ((dirs%input_filename(1:1) == 'n') .and. &
             (len_trim(dirs%input_filename) == 1))   )

  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) * gv%m_to_H ; 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 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, obc=cs%OBC)
  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 create_group_pass(pass_h_tmp, h_tmp, g%Domain)
    call do_group_pass(pass_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

  call cpu_clock_begin(id_clock_pass_init)
  call create_group_pass(pass_av_h_uvh, cs%u_av, cs%v_av, g%Domain, halo=2)
  call create_group_pass(pass_av_h_uvh, cs%h_av, g%Domain, halo=2)
  call create_group_pass(pass_av_h_uvh, uh, vh, g%Domain, halo=2)
  call do_group_pass(pass_av_h_uvh, 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, y_cell_method='sum', v_extensive=.true.)
  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
      'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true.)

  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')

  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')

  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_split_rk2


!> Close the dyn_split_RK2 module
subroutine end_dyn_split_rk2(CS)
  type(mom_dyn_split_rk2_cs), pointer :: CS  !< module control structure

  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)

  call dealloc_bt_cont_type(cs%BT_cont)

  deallocate(cs)
end subroutine end_dyn_split_rk2


!> \namespace mom_dynamics_split_rk2
!!
!!  This file time steps the adiabatic dynamic core by splitting
!!  between baroclinic and barotropic modes. It uses 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). Additional
!!  issues related to exact tracer conservation and how to
!!  ensure consistency between the barotropic and layered estimates
!!  of the free surface height are described 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_split_RK2 actually does the time
!!  stepping, while register_restarts_dyn_split_RK2 sets the fields
!!  that are found in a full restart file with this scheme, and
!!  initialize_dyn_split_RK2 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.

end module mom_dynamics_split_rk2