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 : mech_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_averages
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, pass_var
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, is_new_run, mom_restart_cs
use mom_time_manager,      only : time_type, 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_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_open_boundary,         only : open_boundary_test_extern_h
use mom_pressureforce,         only : pressureforce, pressureforce_init, pressureforce_cs
use mom_set_visc,              only : set_viscous_ml, set_visc_cs
use mom_thickness_diffuse,     only : thickness_diffuse_cs
use mom_tidal_forcing,         only : tidal_forcing_init, tidal_forcing_cs
use mom_unit_scaling,          only : unit_scale_type
use mom_vert_friction,         only : vertvisc, vertvisc_coef, vertvisc_remnant
use mom_vert_friction,         only : 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
use mom_wave_interface, only: wave_parameters_cs
 
implicit none ; private
 
#include <MOM_memory.h>
 
!> MOM_dynamics_split_RK2 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) [L T-2 ~> m s-2]
    pfu, &    !< PFu = -dM/dx [L T-2 ~> m s-2]
    diffu     !< Zonal acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]
 
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
    cav, &    !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2]
    pfv, &    !< PFv = -dM/dy [L T-2 ~> m s-2]
    diffv     !< Meridional acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]
 
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: 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 worth of a barotropic acceleration
              !! that a layer experiences after viscosity is applied.
              !! Nondimensional between 0 (at the bottom) and 1 (far above).
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_accel_bt
              !< The zonal layer accelerations due to the difference between
              !! the barotropic accelerations and the baroclinic accelerations
              !! that were fed into the barotopic calculation [L T-2 ~> m s-2]
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: 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 worth of a barotropic acceleration
              !! that a layer experiences after viscosity is applied.
              !! Nondimensional between 0 (at the bottom) and 1 (far above).
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_accel_bt
              !< The meridional layer accelerations due to the difference between
              !! the barotropic accelerations and the baroclinic accelerations
              !! that were fed into the barotopic calculation [L T-2 ~> 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 Boussinesq
                                                                  !! mode) or column mass anomaly (in non-Boussinesq
                                                                  !! mode) [H ~> m or kg m-2]
  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 [L T-1 ~> 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 [L T-1 ~> m s-1]
  real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_)      :: h_av   !< arithmetic mean of two successive layer
                                                                  !! thicknesses [H ~> m or kg m-2]
  real allocable_, dimension(NIMEM_,NJMEM_)             :: eta_pf !< instantaneous SSH used in calculating PFu and
                                                                  !! PFv [H ~> m or kg m-2]
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_)        :: uhbt   !< average x-volume or mass flux determined by the
                                                                  !! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].
                                                                  !! uhbt is roughly equal to the vertical sum of uh.
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_)        :: vhbt   !< average y-volume or mass flux determined by the
                                                                  !! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].
                                                                  !! vhbt is 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 [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
 
  real, pointer, dimension(:,:) :: taux_bot => null() !< frictional x-bottom stress from the ocean
                                                      !! to the seafloor [R L Z T-2 ~> Pa]
  real, pointer, dimension(:,:) :: tauy_bot => null() !< frictional y-bottom stress from the ocean
                                                      !! to the seafloor [R L Z T-2 ~> 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 :: debug_obc !< If true, do debugging calls for open boundary conditions.
 
  logical :: module_is_initialized = .false. !< Record whether this mouled has been initialzed.
 
  !>@{ Diagnostic IDs
  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.
 
  ! The remainder of the structure points to child subroutines' control structures.
  !> A pointer to the horizontal viscosity control structure
  type(hor_visc_cs),      pointer :: hor_visc_csp      => null()
  !> A pointer to the continuity control structure
  type(continuity_cs),    pointer :: continuity_csp    => null()
  !> A pointer to the CoriolisAdv control structure
  type(coriolisadv_cs),   pointer :: coriolisadv_csp   => null()
  !> A pointer to the PressureForce control structure
  type(pressureforce_cs), pointer :: pressureforce_csp => null()
  !> A pointer to the barotropic stepping control structure
  type(barotropic_cs),    pointer :: barotropic_csp    => null()
  !> A pointer to a structure containing interface height diffusivities
  type(thickness_diffuse_cs), pointer :: thickness_diffuse_csp => null()
  !> A pointer to the vertical viscosity control structure
  type(vertvisc_cs),      pointer :: vertvisc_csp      => null()
  !> A pointer to the set_visc control structure
  type(set_visc_cs),      pointer :: set_visc_csp      => null()
  !> A pointer to the tidal forcing control structure
  type(tidal_forcing_cs), pointer :: tides_csp         => null()
  !> A pointer to the ALE control structure.
  type(ale_cs), pointer :: ale_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.
  !> A pointer to the update_OBC control structure
  type(update_obc_cs),    pointer :: update_obc_csp => null()
 
  type(group_pass_type) :: pass_eta  !< Structure for group halo pass
  type(group_pass_type) :: pass_visc_rem  !< Structure for group halo pass
  type(group_pass_type) :: pass_uvp  !< Structure for group halo pass
  type(group_pass_type) :: pass_hp_uv  !< Structure for group halo pass
  type(group_pass_type) :: pass_uv  !< Structure for group halo pass
  type(group_pass_type) :: pass_h  !< Structure for group halo pass
  type(group_pass_type) :: pass_av_uvh  !< Structure for group halo pass
 
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
 
!>@{ CPU time clock IDs
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, forces, p_surf_begin, p_surf_end, &
                 uh, vh, uhtr, vhtr, eta_av, &
                 G, GV, US, CS, calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, Waves)
  type(ocean_grid_type),             intent(inout) :: g            !< ocean grid structure
  type(verticalgrid_type),           intent(in)    :: gv           !< ocean vertical grid structure
  type(unit_scale_type),             intent(in)    :: us           !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                             target, intent(inout) :: u            !< zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                             target, intent(inout) :: v            !< merid velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                                     intent(inout) :: h            !< layer thickness [H ~> m or kg m-2]
  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 [T ~> s]
  type(mech_forcing),                intent(in)    :: forces       !< A structure with the driving mechanical forces
  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, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                             target, intent(inout) :: uh           !< zonal volume/mass transport
                                                                   !! [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                             target, intent(inout) :: vh           !< merid volume/mass transport
                                                                   !! [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                     intent(inout) :: uhtr         !< accumulatated zonal volume/mass transport
                                                                   !! since last tracer advection [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                     intent(inout) :: vhtr         !< accumulatated merid volume/mass transport
                                                                   !! since last tracer advection [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJ_(G)),  intent(out)   :: eta_av       !< free surface height or column mass time
                                                                   !! averaged over time step [H ~> m or kg m-2]
  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
  type(thickness_diffuse_cs),        pointer       :: thickness_diffuse_csp!< Pointer to a structure containing
                                                                    !! interface height diffusivities
  type(wave_parameters_cs), optional, pointer      :: waves        !< A pointer to a structure containing
                                                                   !! fields related to the surface wave conditions
 
  ! local variables
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up   ! Predicted zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp   ! Predicted meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: hp   ! Predicted thickness [H ~> m or kg m-2].
 
  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 [L T-2 ~> 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 [H L2 T-1 ~> 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 [H m2 s-1 ~> 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,
    ! [H ~> m or kg m-2].
 
  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 [L T-1 ~> 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_av, & ! The zonal velocity time-averaged over a time step [L T-1 ~> m s-1].
    v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1].
    h_av    ! The layer thickness time-averaged over a time step [H ~> m or kg m-2].
  real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
 
  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
 
  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, us, 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 => forces%p_surf
  endif
 
  if (associated(cs%OBC)) then
    if (cs%debug_OBC) call open_boundary_test_extern_h(g, gv, cs%OBC, h)
 
    do k=1,nz ; do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
      u_old_rad_obc(i,j,k) = u_av(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
      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 = &
    (allocated(cs%BT_cont%h_u) .and. allocated(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
 
  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 cpu_clock_begin(id_clock_pass)
  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_averages(dt, time_local, cs%diag)
  call cpu_clock_begin(id_clock_pres)
  call pressureforce(h, tv, cs%PFu, cs%PFv, g, gv, us, 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, us, tv, h, cs%update_OBC_CSp, time_local)
  endif; endif
  if (associated(cs%OBC) .and. cs%debug_OBC) &
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
 
  if (g%nonblocking_updates) &
    call start_group_pass(cs%pass_eta, g%Domain, clock=id_clock_pass)
 
! 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, us, 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, us, 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_averages(dt, time_local, cs%diag)
  call set_viscous_ml(u, v, h, tv, forces, visc, dt, g, gv, us, &
                      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, scale=us%L_T_to_m_s)
  endif
  call vertvisc_coef(up, vp, h, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, 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, .true., g, gv, cs%barotropic_CSp)
  call cpu_clock_end(id_clock_btcalc)
 
  if (g%nonblocking_updates) &
    call complete_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
 
! 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, us, 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
 
  call cpu_clock_begin(id_clock_btstep)
  if (calc_dtbt) call set_dtbt(g, gv, us, 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, forces, 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, us, 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(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(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, scale=us%L_T_to_m_s)
    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*us%L_to_m**2*us%s_to_T)
!   call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
    call mom_accel_chksum("Predictor accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
             cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, symmetric=sym)
    call mom_state_chksum("Predictor 1 init", u, v, h, uh, vh, g, gv, us, 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, scale=us%L_T_to_m_s)
  endif
  call vertvisc_coef(up, vp, h, forces, visc, dt_pred, g, gv, us, cs%vertvisc_CSp, &
                     cs%OBC)
  call vertvisc(up, vp, h, forces, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, g, &
                gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
  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 start_group_pass(cs%pass_uvp, g%Domain, clock=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, us, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)
 
  call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
  else
    call do_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
  endif
 
  ! 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, us, 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)")
 
  call do_group_pass(cs%pass_hp_uv, g%Domain, clock=id_clock_pass)
 
  if (associated(cs%OBC)) then
 
    if (cs%debug) &
      call uvchksum("Pre OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
 
    call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, us, dt_pred)
 
    if (cs%debug) &
      call uvchksum("Post OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
 
    ! These should be done with a pass that excludes uh & vh.
!   call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
  endif
 
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
  endif
 
  ! 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_averages(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, .false., 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, us, 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) &
    call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
 
  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, us, symmetric=sym)
    call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
    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, US)
    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, us, cs%hor_visc_CSp, &
                            obc=cs%OBC, bt=cs%barotropic_CSp)
  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, us, 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, us, 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, forces, 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, us, 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(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(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, scale=us%L_T_to_m_s)
    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*us%L_to_m**2*us%s_to_T)
  ! call MOM_state_chksum("Corrector 1", u, v, h, uh, vh, G, GV, US, haloshift=1)
    call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us, 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, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(u, v, h, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, us, &
                cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot,waves=waves)
  if (g%nonblocking_updates) then
    call cpu_clock_end(id_clock_vertvisc)
    call start_group_pass(cs%pass_uv, g%Domain, clock=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, us, 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 do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
  else
    call do_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
  endif
 
  ! 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, us, 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 do_group_pass(cs%pass_h, g%Domain, clock=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)")
 
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
  else
    call do_group_pass(cs%pass_av_uvh, g%domain, clock=id_clock_pass)
  endif
 
  if (associated(cs%OBC)) then
    call radiation_open_bdry_conds(cs%OBC, u, u_old_rad_obc, v, v_old_rad_obc, g, us, 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) &
    call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
 
  !$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, us, symmetric=sym)
    call uvchksum("Corrector avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
    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, US)
  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 [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
                         target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
 
  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_H
 
  thickness_units = get_thickness_units(gv)
  flux_units = get_flux_units(gv)
 
  if (gv%Boussinesq) then
    vd = var_desc("sfc",thickness_units,"Free surface Height",'h','1')
  else
    vd = var_desc("p_bot",thickness_units,"Bottom Pressure",'h','1')
  endif
  call register_restart_field(cs%eta, vd, .false., restart_cs)
 
  vd = var_desc("u2","m s-1","Auxiliary Zonal velocity",'u','L')
  call register_restart_field(cs%u_av, vd, .false., restart_cs)
 
  vd = var_desc("v2","m s-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","m s-2","Zonal horizontal viscous acceleration",'u','L')
  call register_restart_field(cs%diffu, vd, .false., restart_cs)
 
  vd = var_desc("diffv","m s-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, US, param_file, &
                      diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
                      VarMix, MEKE, thickness_diffuse_CSp,                  &
                      OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
                      visc, dirs, ntrunc, calc_dtbt)
  type(ocean_grid_type),            intent(inout) :: g          !< ocean grid structure
  type(verticalgrid_type),          intent(in)    :: gv         !< ocean vertical grid structure
  type(unit_scale_type),            intent(in)    :: us         !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                    intent(inout) :: u          !< zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                    intent(inout) :: v          !< merid velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                            target, intent(inout) :: uh    !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                            target, intent(inout) :: vh    !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta        !< free surface height or column mass [H ~> 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 [T ~> s]
  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(Barotropic_CS),              pointer       :: Barotropic_CSp !< Pointer to the control structure for
!                                                                !! the barotropic module
  type(thickness_diffuse_cs),       pointer       :: thickness_diffuse_csp !< Pointer to the control structure
                                                  !! used for the isopycnal height diffusive transport.
  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).
  logical,                          intent(out)   :: calc_dtbt  !< If true, recalculate the barotropic time step
 
  ! local variables
  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, eta_rest_name
  real :: h_rescale  ! A rescaling factor for thicknesses from the representation in
                     ! a restart file to the internal representation in this run.
  real :: vel_rescale  ! A rescaling factor for velocities from the representation in
                     ! a restart file to the internal representation in this run.
  real :: uh_rescale ! A rescaling factor for thickness transports from the representation in
                     ! a restart file to the internal representation in this run.
  real :: accel_rescale ! A rescaling factor for accelerations from the representation in
                     ! a restart file to the internal representation in this run.
  real :: h_convert
  type(group_pass_type) :: 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 "//&
                 "of a  2nd-order Runga-Kutta baroclinic time stepping "//&
                 "scheme (0.5) and a backward Euler scheme (1) that is "//&
                 "used for the Coriolis and inertial terms.  BE may be "//&
                 "from 0.5 to 1, but instability may occur near 0.5. "//&
                 "BE is also applicable if SPLIT is false and USE_RK2 "//&
                 "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 "//&
                 "controls the extent to which the treatment of gravity "//&
                 "waves is forward-backward (0) or simulated backward "//&
                 "Euler (1).  0 is almost always used. "//&
                 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
                 "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 "//&
                 "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 "//&
                 "adjustment due to the change in the barotropic velocity "//&
                 "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., debuggingparam=.true.)
  call get_param(param_file, mdl, "DEBUG_OBC", cs%debug_OBC, 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
 
  id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', &
                                     grain=clock_routine)
 
  call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp)
  call coriolisadv_init(time, g, gv, us, 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, us, param_file, diag, cs%PressureForce_CSp, &
                          cs%tides_CSp)
  call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke)
  call vertvisc_init(mis, time, g, gv, us, 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=is_new_run(restart_cs) )
 
  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
 
  eta_rest_name = "sfc" ; if (.not.gv%Boussinesq) eta_rest_name = "p_bot"
  if (.not. query_initialized(cs%eta, trim(eta_rest_name), 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) = -gv%Z_to_H * 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
  elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
    h_rescale = gv%m_to_H / gv%m_to_H_restart
    do j=js,je ; do i=is,ie ; cs%eta(i,j) = h_rescale * cs%eta(i,j) ; 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, us, param_file, diag, &
                       cs%barotropic_CSp, restart_cs, calc_dtbt, cs%BT_cont, &
                       cs%tides_CSp)
 
  if (.not. query_initialized(cs%diffu,"diffu",restart_cs) .or. &
      .not. query_initialized(cs%diffv,"diffv",restart_cs)) then
    call horizontal_viscosity(u, v, h, cs%diffu, cs%diffv, meke, varmix, &
                              g, gv, us, cs%hor_visc_CSp, &
                              obc=cs%OBC, bt=cs%barotropic_CSp)
  else
    if ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
         (us%m_to_L * us%s_to_T_restart**2 /= us%m_to_L_restart * us%s_to_T**2) ) then
      accel_rescale = (us%m_to_L * us%s_to_T_restart**2) / (us%m_to_L_restart * us%s_to_T**2)
      do k=1,nz ; do j=js,je ; do i=g%IscB,g%IecB
        cs%diffu(i,j,k) = accel_rescale * cs%diffu(i,j,k)
      enddo ; enddo ; enddo
      do k=1,nz ; do j=g%JscB,g%JecB ; do i=is,ie
        cs%diffv(i,j,k) = accel_rescale * cs%diffv(i,j,k)
      enddo ; enddo ; enddo
    endif
  endif
 
  if (.not. query_initialized(cs%u_av,"u2", restart_cs) .or. &
      .not. query_initialized(cs%u_av,"v2", restart_cs)) then
    do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb ; cs%u_av(i,j,k) = u(i,j,k) ; enddo ; enddo ; enddo
    do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied ; cs%v_av(i,j,k) = v(i,j,k) ; enddo ; enddo ; enddo
  elseif ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
         ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) ) then
    vel_rescale = (us%m_to_L * us%s_to_T_restart) /  (us%m_to_L_restart * us%s_to_T)
    do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb ; cs%u_av(i,j,k) = vel_rescale * cs%u_av(i,j,k) ; enddo ; enddo ; enddo
    do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied ; cs%v_av(i,j,k) = vel_rescale * cs%v_av(i,j,k) ; enddo ; enddo ; enddo
  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
    do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo
    call continuity(u, v, h, h_tmp, uh, vh, dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
    call pass_var(h_tmp, g%Domain, clock=id_clock_pass_init)
    do k=1,nz ; do j=jsd,jed ; do i=isd,ied
      cs%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
    enddo ; enddo ; enddo
  else
    if (.not. query_initialized(cs%h_av,"h2",restart_cs)) then
      cs%h_av(:,:,:) = h(:,:,:)
    elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
      h_rescale = gv%m_to_H / gv%m_to_H_restart
      do k=1,nz ; do j=js,je ; do i=is,ie ; cs%h_av(i,j,k) = h_rescale * cs%h_av(i,j,k) ; enddo ; enddo ; enddo
    endif
    if ( (gv%m_to_H_restart * us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
         ((gv%m_to_H * us%m_to_L**2 * us%s_to_T_restart) /= &
          (gv%m_to_H_restart * us%m_to_L_restart**2 * us%s_to_T)) ) then
      uh_rescale = (gv%m_to_H * us%m_to_L**2 * us%s_to_T_restart) / &
                   (gv%m_to_H_restart * us%m_to_L_restart**2 * us%s_to_T)
      do k=1,nz ; do j=js,je ; do i=g%IscB,g%IecB ; uh(i,j,k) = uh_rescale * uh(i,j,k) ; enddo ; enddo ; enddo
      do k=1,nz ; do j=g%JscB,g%JecB ; do i=is,ie ; vh(i,j,k) = uh_rescale * vh(i,j,k) ; enddo ; enddo ; enddo
    endif
  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)
  h_convert = gv%H_to_m ; if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
      'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extensive=.true., &
      conversion=h_convert*us%L_to_m**2*us%s_to_T)
  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
      'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true., &
      conversion=h_convert*us%L_to_m**2*us%s_to_T)
 
  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
      'Zonal Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
      'Meridional Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
      'Zonal Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
      'Meridional Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
 
  cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
      'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
      'Barotropic-step Averaged Meridional Velocity', 'm s-1', conversion=us%L_T_to_m_s)
 
  cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
    'Barotropic Anomaly Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
    'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
 
  id_clock_cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
  id_clock_continuity = cpu_clock_id('(Ocean continuity equation)',      grain=clock_module)
  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_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