MOM_barotropic.F90

Go to the documentation of this file.

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

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, April 1994 - January 2007                      *
!*                                                                     *
!*    This program contains the subroutines that time steps the        *
!*  linearized barotropic equations.  btstep is used to actually       *
!*  time step the barotropic equations, and contains most of the       *
!*  substance of this module.                                          *
!*                                                                     *
!*    btstep uses a forwards-backwards based scheme to time step       *
!*  the barotropic equations, returning the layers' accelerations due  *
!*  to the barotropic changes in the ocean state, the final free       *
!*  surface height (or column mass), and the volume (or mass) fluxes   *
!*  summed through the layers and averaged over the baroclinic time    *
!*  step.  As input, btstep takes the initial 3-D velocities, the      *
!*  inital free surface height, the 3-D accelerations of the layers,   *
!*  and the external forcing.  Everything in btstep is cast in terms   *
!*  of anomalies, so if everything is in balance, there is explicitly  *
!*  no acceleration due to btstep.                                     *
!*                                                                     *
!*    The spatial discretization of the continuity equation is second  *
!*  order accurate.  A flux conservative form is used to guarantee     *
!*  global conservation of volume.  The spatial discretization of the  *
!*  momentum equation is second order accurate.  The Coriolis force    *
!*  is written in a form which does not contribute to the energy       *
!*  tendency and which conserves linearized potential vorticity, f/D.  *
!*  These terms are exactly removed from the baroclinic momentum       *
!*  equations, so the linearization of vorticity advection will not    *
!*  degrade the overall solution.                                      *
!*                                                                     *
!*    btcalc calculates the fractional thickness of each layer at the  *
!*  velocity points, for later use in calculating the barotropic       *
!*  velocities and the averaged accelerations.  Harmonic mean          *
!*  thicknesses (i.e. 2*h_L*h_R/(h_L + h_R)) are used to avoid overly  *
!*  strong weighting of overly thin layers.  This may later be relaxed *
!*  to use thicknesses determined from the continuity equations.       *
!*                                                                     *
!*    bt_mass_source determines the real mass sources for the          *
!*  barotropic solver, along with the corrective pseudo-fluxes that    *
!*  keep the barotropic and baroclinic estimates of the free surface   *
!*  height close to each other.  Given the layer thicknesses and the   *
!*  free surface height that correspond to each other, it calculates   *
!*  a corrective mass source that is added to the barotropic continuity*
!*  equation, and optionally adjusts a slowly varying correction rate. *
!*  Newer algorithmic changes have deemphasized the need for this, but *
!*  it is still here to add net water sources to the barotropic solver.*
!*                                                                     *
!*    barotropic_init allocates and initializes any barotropic arrays  *
!*  that have not been read from a restart file, reads parameters from *
!*  the inputfile, and sets up diagnostic fields.                      *
!*                                                                     *
!*    barotropic_end deallocates anything allocated in barotropic_init *
!*  or register_barotropic_restarts.                                   *
!*                                                                     *
!*    register_barotropic_restarts is used to indicate any fields that *
!*  are private to the barotropic solver that need to be included in   *
!*  the restart files, and to ensure that they are read.               *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, CoriolisBu                            *
!*    j+1  > o > o >   At ^:  v_in, vbt, accel_layer_v, vbtav          *
!*    j    x ^ x ^ x   At >:  u_in, ubt, accel_layer_u, ubtav, amer    *
!*    j    > o > o >   At o:  eta, h, bathyT, pbce                     *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1                                                  *
!*           i  i+1                                                    *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_debugging, only : hchksum, uvchksum
use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
use mom_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field
use mom_diag_mediator, only : safe_alloc_ptr, diag_ctrl, enable_averaging
use mom_domains, only : min_across_pes, clone_mom_domain, pass_vector
use mom_domains, only : to_all, scalar_pair, agrid, corner, mom_domain_type
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_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_grid, only : ocean_grid_type
use mom_hor_index, only : hor_index_type
use mom_io, only : vardesc, var_desc
use mom_open_boundary, only : ocean_obc_type, obc_simple, obc_none
use mom_open_boundary, only : obc_direction_e, obc_direction_w
use mom_open_boundary, only : obc_direction_n, obc_direction_s, obc_segment_type
use mom_restart, only : register_restart_field, query_initialized, mom_restart_cs
use mom_tidal_forcing, only : tidal_forcing_sensitivity, tidal_forcing_cs
use mom_time_manager, only : time_type, set_time, operator(+), operator(-)
use mom_variables, only : bt_cont_type, alloc_bt_cont_type
use mom_verticalgrid, only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>
#ifdef STATIC_MEMORY_
#  ifndef BTHALO_
#    define BTHALO_ 0
#  endif
#  define WHALOI_ MAX(BTHALO_-NIHALO_,0)
#  define WHALOJ_ MAX(BTHALO_-NJHALO_,0)
#  define NIMEMW_   1-WHALOI_:NIMEM_+WHALOI_
#  define NJMEMW_   1-WHALOJ_:NJMEM_+WHALOJ_
#  define NIMEMBW_  -WHALOI_:NIMEM_+WHALOI_
#  define NJMEMBW_  -WHALOJ_:NJMEM_+WHALOJ_
#  define SZIW_(G)  NIMEMW_
#  define SZJW_(G)  NJMEMW_
#  define SZIBW_(G) NIMEMBW_
#  define SZJBW_(G) NJMEMBW_
#else
#  define NIMEMW_   :
#  define NJMEMW_   :
#  define NIMEMBW_  :
#  define NJMEMBW_  :
#  define SZIW_(G)  G%isdw:G%iedw
#  define SZJW_(G)  G%jsdw:G%jedw
#  define SZIBW_(G) G%isdw-1:G%iedw
#  define SZJBW_(G) G%jsdw-1:G%jedw
#endif

public btcalc, bt_mass_source, btstep, barotropic_init, barotropic_end
public register_barotropic_restarts, set_dtbt

type, private :: bt_obc_type
  real, dimension(:,:), pointer :: &
    cg_u => null(), &     ! The external wave speed at u-points, in m s-1.
    cg_v => null(), &     ! The external wave speed at u-points, in m s-1.
    h_u => null(), &      ! The total thickness at the u-points, in m or kg m-2.
    h_v => null(), &      ! The total thickness at the v-points, in m or kg m-2.
    uhbt => null(), &     ! The zonal and meridional barotropic thickness fluxes
    vhbt => null(), &     ! specified for open boundary conditions (if any),
                          ! in units of m3 s-1.
    ubt_outer => null(), & ! The zonal and meridional velocities just outside
    vbt_outer => null(), & ! the domain, as set by the open boundary conditions,
                           ! in units of m s-1.
    eta_outer_u => null(), & ! The surface height outside of the domain at a
    eta_outer_v => null()    ! u- or v- point with an open boundary condition,
                             ! in units of m or kg m-2.
  logical :: apply_u_obcs !< True if this PE has an open boundary at a u-point.
  logical :: apply_v_obcs !< True if this PE has an open boundary at a v-point.
  integer :: is_u_obc, ie_u_obc, js_u_obc, je_u_obc
  integer :: is_v_obc, ie_v_obc, js_v_obc, je_v_obc
  logical :: is_alloced = .false. !< True if BT_OBC is in use and has been allocated
end type bt_obc_type

type, public :: barotropic_cs ; private
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: frhatu
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: frhatv
      ! frhatu and frhatv are the fraction of the total column thickness
      ! interpolated to u or v grid points in each layer, nondimensional.
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
    idatu, &        ! Inverse of the basin depth at u grid points, in m-1.
    uhbt_ic, &      ! The barotropic solver's estimate of the zonal
                    ! transport as the initla condition for the next call
                    ! to btstep, in H m2 s-1.
    ubt_ic, &       ! The barotropic solver's estimate of the zonal velocity
                    ! that will be the initial condition for the next call
                    ! to btstep, in m s-1.
    ubtav           ! The barotropic zonal velocity averaged over the
                    ! baroclinic time step, m s-1.
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
    idatv, &        ! Inverse of the basin depth at v grid points, in m-1.
    vhbt_ic, &      ! The barotropic solver's estimate of the zonal
                    ! transport as the initla condition for the next call
                    ! to btstep, in H m2 s-1.
    vbt_ic, &       ! The barotropic solver's estimate of the zonal velocity
                    ! that will be the initial condition for the next call
                    ! to btstep, in m s-1.
    vbtav           ! The barotropic meridional velocity averaged over the
                    ! baroclinic time step, m s-1.
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    eta_source, &   ! The net mass source to be applied within the
                    ! barotropic solver, in H s-1.
    eta_cor, &      ! The difference between the free surface height from
                    ! the barotropic calculation and the sum of the layer
                    ! thicknesses. This difference is imposed as a forcing
                    ! term in the barotropic calculation over a baroclinic
                    ! timestep, in H (m or kg m-2).
    eta_cor_bound   ! A limit on the rate at which eta_cor can be applied
                    ! while avoiding instability, in units of H s-1. This
                    ! is only used if CS%bound_BT_corr is true.
  real allocable_, dimension(NIMEMW_,NJMEMW_) :: &
    ua_polarity, &  ! Test vector components for checking grid polarity.
    va_polarity, &  ! Test vector components for checking grid polarity.
    bathyt          !   A copy of bathyT (ocean bottom depth) with wide halos.
  real allocable_, dimension(NIMEMW_,NJMEMW_) :: &
    iareat          !   This is a copy of G%IareaT with wide halos, but will
                    ! still utilize the macro IareaT when referenced, m-2.
  real allocable_, dimension(NIMEMBW_,NJMEMW_) :: &
    d_u_cor, &      !   A simply averaged depth at u points, in m.
    dy_cu, &         !   A copy of G%dy_Cu with wide halos, in m.
    idxcu            !   A copy of G%IdxCu with wide halos, in m-1.
  real allocable_, dimension(NIMEMW_,NJMEMBW_) :: &
    d_v_cor, &      !   A simply averaged depth at v points, in m.
    dx_cv, &         !   A copy of G%dx_Cv with wide halos, in m.
    idycv            !   A copy of G%IdyCv with wide halos, in m-1.
  real allocable_, dimension(NIMEMBW_,NJMEMBW_) :: &
    q_d             ! f / D at PV points, in m-1 s-1.

  real, pointer, dimension(:,:,:) :: frhatu1 => null(), frhatv1 => null() ! Predictor values.

  type(bt_obc_type) :: bt_obc !< A structure with all of this module's fields
                              !! for applying open boundary conditions.

  real    :: rho0            !   The density used in the Boussinesq
                             ! approximation, in kg m-3.
  real    :: dtbt            ! The barotropic time step, in s.
  real    :: dtbt_fraction   !   The fraction of the maximum time-step that
                             ! should used.  The default is 0.98.
  real    :: dtbt_max        !   The maximum stable barotropic time step, in s.
  real    :: dt_bt_filter    !   The time-scale over which the barotropic mode
                             ! solutions are filtered, in s.  This can never
                             ! be taken to be longer than 2*dt.  The default, 0,
                             ! applies no filtering.
  integer :: nstep_last = 0  ! The number of barotropic timesteps per baroclinic
                             ! time step the last time btstep was called.
  real    :: bebt            ! A nondimensional number, from 0 to 1, that
                             ! determines the gravity wave time stepping scheme.
                             ! 0.0 gives a forward-backward scheme, while 1.0
                             ! give backward Euler. In practice, bebt should be
                             ! of order 0.2 or greater.
  logical :: split           ! If true, use the split time stepping scheme.
  real    :: eta_source_limit  ! The fraction of the initial depth of the ocean
                             ! that can be added or removed to the bartropic
                             ! solution within a thermodynamic time step.  By
                             ! default this is 0 (i.e., no correction). Nondim.
  logical :: bound_bt_corr   ! If true, the magnitude of the fake mass source
                             ! in the barotropic equation that drives the two
                             ! estimates of the free surface height toward each
                             ! other is bounded to avoid driving corrective
                             ! velocities that exceed MAXCFL_BT_CONT.
  logical :: gradual_bt_ics  ! If true, adjust the initial conditions for the
                             ! barotropic solver to the values from the layered
                             ! solution over a whole timestep instead of
                             ! instantly.  This is a decent approximation to the
                             ! inclusion of sum(u dh_dt) while also correcting
                             ! for truncation errors.
  logical :: sadourny        ! If true, the Coriolis terms are discretized
                             ! with Sadourny's energy conserving scheme,
                             ! otherwise the Arakawa & Hsu scheme is used.  If
                             ! the deformation radius is not resolved Sadourny's
                             ! scheme should probably be used.
  logical :: nonlinear_continuity ! If true, the barotropic continuity equation
                             ! uses the full ocean thickness for transport.
  integer :: nonlin_cont_update_period ! The number of barotropic time steps
                             ! between updates to the face area, or 0 only to
                             ! update at the start of a call to btstep.  The
                             ! default is 1.
  logical :: bt_project_velocity ! If true, step the barotropic velocity first
                             ! and project out the velocity tendancy by 1+BEBT
                             ! when calculating the transport.  The default
                             ! (false) is to use a predictor continuity step to
                             ! find the pressure field, and then do a corrector
                             ! continuity step using a weighted average of the
                             ! old and new velocities, with weights of (1-BEBT)
                             ! and BEBT.
  logical :: dynamic_psurf   ! If true, add a dynamic pressure due to a viscous
                             ! ice shelf, for instance.
  real    :: dmin_dyn_psurf  ! The minimum depth to use in limiting the size
                             ! of the dynamic surface pressure for stability,
                             ! in m.
  real    :: ice_strength_length  ! The length scale at which the damping rate
                             ! due to the ice strength should be the same as if
                             ! a Laplacian were applied, in m.
  real    :: const_dyn_psurf ! The constant that scales the dynamic surface
                             ! pressure, nondim.  Stable values are < ~1.0.
                             ! The default is 0.9.
  logical :: tides           ! If true, apply tidal momentum forcing.
  real    :: g_extra         ! A nondimensional factor by which gtot is enhanced.
  integer :: hvel_scheme     ! An integer indicating how the thicknesses at
                             ! velocity points are calculated. Valid values are
                             ! given by the parameters defined below:
                             !   HARMONIC, ARITHMETIC, HYBRID, and FROM_BT_CONT
  logical :: strong_drag     ! If true, use a stronger estimate of the retarding
                             ! effects of strong bottom drag.
  logical :: linearized_bt_pv  ! If true, the PV and interface thicknesses used
                             ! in the barotropic Coriolis calculation is time
                             ! invariant and linearized.
  logical :: use_wide_halos  ! If true, use wide halos and march in during the
                             ! barotropic time stepping for efficiency.
  logical :: clip_velocity   ! If true, limit any velocity components that are
                             ! are large enough for a CFL number to exceed
                             ! CFL_trunc.  This should only be used as a
                             ! desperate debugging measure.
  logical :: debug           ! If true, write verbose checksums for debugging purposes.
  logical :: debug_bt        ! If true, write verbose checksums for debugging purposes.
  real    :: maxvel          ! Velocity components greater than maxvel are
                             ! truncated to maxvel, in m s-1.
  real    :: cfl_trunc       ! If clip_velocity is true, velocity components will
                             ! be truncated when they are large enough that the
                             ! corresponding CFL number exceeds this value, nondim.
  real    :: maxcfl_bt_cont  ! The maximum permitted CFL number associated with the
                             ! barotropic accelerations from the summed velocities
                             ! times the time-derivatives of thicknesses.  The
                             ! default is 0.1, and there will probably be real
                             ! problems if this were set close to 1.
  logical :: bt_cont_bounds  ! If true, use the BT_cont_type variables to set
                             ! limits on the magnitude of the corrective mass
                             ! fluxes.
  logical :: visc_rem_u_uh0  ! If true, use the viscous remnants when estimating
                             ! the barotropic velocities that were used to
                             ! calculate uh0 and vh0.  False is probably the
                             ! better choice.
  logical :: adjust_bt_cont  ! If true, adjust the curve fit to the BT_cont type
                             ! that is used by the barotropic solver to match the
                             ! transport about which the flow is being linearized.
  type(time_type), pointer :: time ! A pointer to the ocean model's clock.
  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                             ! timing of diagnostic output.
  type(mom_domain_type), pointer :: bt_domain => null()
  type(hor_index_type), pointer :: debug_bt_hi ! debugging copy of horizontal index_type
  type(tidal_forcing_cs), pointer :: tides_csp => null()
  logical :: module_is_initialized = .false.

  integer :: isdw, iedw, jsdw, jedw ! The memory limits of the wide halo arrays.

  !--- for group halo pass
  type(group_pass_type) :: pass_q_dcor, pass_gtot
  type(group_pass_type) :: pass_tmp_uv, pass_eta_bt_rem
  type(group_pass_type) :: pass_force_hbt0_cor_ref, pass_dat_uv
  type(group_pass_type) :: pass_eta_ubt, pass_etaav, pass_ubt_cor
  type(group_pass_type) :: pass_ubta_uhbta, pass_e_anom

  integer :: id_pfu_bt = -1, id_pfv_bt = -1, id_coru_bt = -1, id_corv_bt = -1
  integer :: id_ubtforce = -1, id_vbtforce = -1, id_uaccel = -1, id_vaccel = -1
  integer :: id_visc_rem_u = -1, id_visc_rem_v = -1, id_eta_cor = -1
  integer :: id_ubt = -1, id_vbt = -1, id_eta_bt = -1, id_ubtav = -1, id_vbtav = -1
  integer :: id_ubt_st = -1, id_vbt_st = -1, id_eta_st = -1
  integer :: id_ubt_hifreq = -1, id_vbt_hifreq = -1, id_eta_hifreq = -1
  integer :: id_uhbt_hifreq = -1, id_vhbt_hifreq = -1, id_eta_pred_hifreq = -1
  integer :: id_gtotn = -1, id_gtots = -1, id_gtote = -1, id_gtotw = -1
  integer :: id_uhbt = -1, id_frhatu = -1, id_vhbt = -1, id_frhatv = -1
  integer :: id_frhatu1 = -1, id_frhatv1 = -1

  integer :: id_btc_fa_u_ee = -1, id_btc_fa_u_e0 = -1, id_btc_fa_u_w0 = -1, id_btc_fa_u_ww = -1
  integer :: id_btc_ubt_ee = -1, id_btc_ubt_ww = -1
  integer :: id_btc_fa_v_nn = -1, id_btc_fa_v_n0 = -1, id_btc_fa_v_s0 = -1, id_btc_fa_v_ss = -1
  integer :: id_btc_vbt_nn = -1, id_btc_vbt_ss = -1
  integer :: id_uhbt0 = -1, id_vhbt0 = -1

end type barotropic_cs

type, private :: local_bt_cont_u_type
  real :: fa_u_ee, fa_u_e0, fa_u_w0, fa_u_ww
  real :: ubt_ee, ubt_ww
  real :: uh_crve, uh_crvw
  real :: uh_ee, uh_ww
end type local_bt_cont_u_type
type, private :: local_bt_cont_v_type
  real :: fa_v_nn, fa_v_n0, fa_v_s0, fa_v_ss
  real :: vbt_nn, vbt_ss
  real :: vh_crvn, vh_crvs
  real :: vh_nn, vh_ss
end type local_bt_cont_v_type

type, private :: memory_size_type
  integer :: isdw, iedw, jsdw, jedw ! The memory limits of the wide halo arrays.
end type memory_size_type

integer :: id_clock_sync=-1, id_clock_calc=-1
integer :: id_clock_calc_pre=-1, id_clock_calc_post=-1
integer :: id_clock_pass_step=-1, id_clock_pass_pre=-1, id_clock_pass_post=-1

! Enumeration values for various schemes
integer, parameter :: harmonic        = 1
integer, parameter :: arithmetic      = 2
integer, parameter :: hybrid          = 3
integer, parameter :: from_bt_cont    = 4
integer, parameter :: hybrid_bt_cont  = 5
character*(20), parameter :: hybrid_string = "HYBRID"
character*(20), parameter :: harmonic_string = "HARMONIC"
character*(20), parameter :: arithmetic_string = "ARITHMETIC"
character*(20), parameter :: bt_cont_string = "FROM_BT_CONT"

contains

!> This subroutine time steps the barotropic equations explicitly.
!! For gravity waves, anything between a forwards-backwards scheme
!! and a simulated backwards Euler scheme is used, with bebt between
!! 0.0 and 1.0 determining the scheme.  In practice, bebt must be of
!! order 0.2 or greater.  A forwards-backwards treatment of the
!! Coriolis terms is always used.
subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
                  fluxes, pbce, eta_PF_in, U_Cor, V_Cor, &
                  accel_layer_u, accel_layer_v, eta_out, uhbtav, vhbtav, G, GV, CS, &
                  visc_rem_u, visc_rem_v, etaav, OBC, &
                  BT_cont, eta_PF_start, &
                  taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0)
  type(ocean_grid_type),                   intent(inout) :: G       !< The ocean's grid structure.
  type(verticalgrid_type),                   intent(in)  :: GV      !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)  :: U_in    !< The initial (3-D) zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)  :: V_in    !< The initial (3-D) meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G)),          intent(in)  :: eta_in  !< The initial barotropic free surface height
                                                         !! anomaly or column mass anomaly, in H (m or kg m-2).
  real,                                      intent(in)  :: dt      !< The time increment to integrate over.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)  :: bc_accel_u !< The zonal baroclinic accelerations, in m s-2.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)  :: bc_accel_v !< The meridional baroclinic accelerations,
                                                                       !! in m s-2.
  type(forcing),                             intent(in)  :: fluxes     !< A structure containing pointers to any
                                                         !! possible forcing fields.  Unused fields have NULL ptrs.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)  :: pbce       !< The baroclinic pressure anomaly in each layer
                                                         !! due to free surface height anomalies, in m2 H-1 s-2.
  real, dimension(SZI_(G),SZJ_(G)),          intent(in)  :: eta_PF_in  !< The 2-D eta field (either SSH anomaly or
                                                         !! column mass anomaly) that was used to calculate the input
                                                         !! pressure gradient accelerations (or its final value if
                                                         !! eta_PF_start is provided, in m or kg m-2.
                                                         !! Note: eta_in, pbce, and eta_PF_in must have up-to-date
                                                         !! values in the first point of their halos.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)  :: U_Cor      !< The (3-D) zonal-velocities used to
                                                         !! calculate the Coriolis terms in bc_accel_u, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)  :: V_Cor      !< Ditto for meridonal bc_accel_v.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: accel_layer_u !< The zonal acceleration of each layer due
                                                         !! to the barotropic calculation, in m s-2.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: accel_layer_v !< The meridional acceleration of each layer
                                                         !! due to the barotropic calculation, in m s-2.
  real, dimension(SZI_(G),SZJ_(G)),          intent(out) :: eta_out       !< The final barotropic free surface
                                                         !! height anomaly or column mass anomaly, in m or kg m-2.
  real, dimension(SZIB_(G),SZJ_(G)),         intent(out) :: uhbtav        !< the barotropic zonal volume or mass
                                                         !! fluxes averaged through the barotropic steps, in
                                                         !! m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJB_(G)),         intent(out) :: vhbtav        !< the barotropic meridional volume or mass
                                                         !! fluxes averaged through the barotropic steps, in
                                                         !! m3 s-1 or kg s-1.
  type(barotropic_cs),                       pointer     :: CS            !< The control structure returned by a
                                                         !! previous call to barotropic_init.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)  :: visc_rem_u    !< Both the fraction of the 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, in the zonal direction. Nondimensional
                                                         !! between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)  :: visc_rem_v    !< Ditto for meridional direction.
  real, dimension(SZI_(G),SZJ_(G)), intent(out), optional :: etaav        !< The free surface height or column mass
                                                         !! averaged over the barotropic integration, in m or kg m-2.
  type(ocean_obc_type),             pointer,     optional :: OBC          !< The open boundary condition structure.
  type(bt_cont_type),               pointer,     optional :: BT_cont      !< A structure with elements that describe
                                                         !! the effective open face areas as a function of barotropic
                                                         !! flow.
  real, dimension(:,:),             pointer,     optional :: eta_PF_start !< The eta field consistent with the pressure
                                                         !! gradient at the start of the barotropic stepping, in m or
                                                         !! kg m-2.
  real, dimension(:,:),             pointer,     optional :: taux_bot     !< The zonal bottom frictional stress from
                                                         !! ocean to the seafloor, in Pa.
  real, dimension(:,:),             pointer,     optional :: tauy_bot     !< The meridional bottom frictional stress
                                                         !! from ocean to the seafloor, in Pa.
  real, dimension(:,:,:),           pointer,     optional :: uh0, u_uh0
  real, dimension(:,:,:),           pointer,     optional :: vh0, v_vh0

  ! Local variables
  real :: ubt_Cor(szib_(g),szj_(g)) ! The barotropic velocities that had been
  real :: vbt_Cor(szi_(g),szjb_(g)) ! used to calculate the input Coriolis
                                    ! terms, in m s-1.
  real :: wt_u(szib_(g),szj_(g),szk_(g)) ! wt_u and wt_v are the
  real :: wt_v(szi_(g),szjb_(g),szk_(g)) ! normalized weights to
                ! be used in calculating barotropic velocities, possibly with
                ! sums less than one due to viscous losses.  Nondimensional.
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    av_rem_u, &   ! The weighted average of visc_rem_u, nondimensional.
    tmp_u         ! A temporary array at u points.
  real, dimension(SZI_(G),SZJB_(G)) :: &
    av_rem_v, &   ! The weighted average of visc_rem_v, nondimensional.
    tmp_v         ! A temporary array at v points.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    e_anom        ! The anomaly in the sea surface height or column mass
                  ! averaged between the beginning and end of the time step,
                  ! relative to eta_PF, with SAL effects included, in units
                  ! of H (m or kg m-2, the same as eta and h).

  ! These are always allocated with symmetric memory and wide halos.
  real :: q(szibw_(cs),szjbw_(cs))  ! A pseudo potential vorticity in s-1 m-1.
  real, dimension(SZIBW_(CS),SZJW_(CS)) :: &
    ubt, &        ! The zonal barotropic velocity in m s-1.
    bt_rem_u, &   ! The fraction of the barotropic zonal velocity that remains
                  ! after a time step, the remainder being lost to bottom drag.
                  ! bt_rem_u is a nondimensional number between 0 and 1.
    bt_force_u, & ! The vertical average of all of the u-accelerations that are
                  ! not explicitly included in the barotropic equation, m s-2.
    u_accel_bt, & ! The difference between the zonal acceleration from the
                  ! barotropic calculation and BT_force_u, in m s-2.
    uhbt, &       ! The zonal barotropic thickness fluxes, in H m2 s-1.
    uhbt0, &      ! The difference between the sum of the layer zonal thickness
                  ! fluxes and the barotropic thickness flux using the same
                  ! velocity, in H m2 s-1.
    ubt_old, &    ! The starting value of ubt in a barotropic step, in m s-1.
    ubt_first, &  ! The starting value of ubt in a series of barotropic steps, in m s-1.
    ubt_sum, &    ! The sum of ubt over the time steps, in m s-1.
    uhbt_sum, &   ! The sum of uhbt over the time steps, in H m2 s-1.
    ubt_wtd, &    ! A weighted sum used to find the filtered final ubt, in m s-1.
    ubt_trans, &  ! The latest value of ubt used for a transport, in m s-1.
    azon, bzon, & ! _zon & _mer are the values of the Coriolis force which
    czon, dzon, & ! are applied to the neighboring values of vbtav & ubtav,
    amer, bmer, & ! respectively to get the barotropic inertial rotation,
    cmer, dmer, & ! in units of s-1.
    cor_u, &      ! The zonal Coriolis acceleration, in m s-2.
    cor_ref_u, &  ! The zonal barotropic Coriolis acceleration due
                  ! to the reference velocities, in m s-2.
    pfu, &        ! The zonal pressure force acceleration, in m s-2.
    pfu_bt_sum, & ! The summed zonal barotropic pressure gradient force, in m s-2.
    coru_bt_sum, & ! The summed zonal barotropic Coriolis acceleration, in m s-2.
    dcor_u, &     ! A simply averaged depth at u points, in m.
    datu          ! Basin depth at u-velocity grid points times the y-grid
                  ! spacing, in H m.
  real, dimension(SZIW_(CS),SZJBW_(CS)) :: &
    vbt, &        ! The meridional barotropic velocity in m s-1.
    bt_rem_v, &   ! The fraction of the barotropic meridional velocity that
                  ! remains after a time step, the rest being lost to bottom
                  ! drag.  bt_rem_v is a nondimensional number between 0 and 1.
    bt_force_v, & ! The vertical average of all of the v-accelerations that are
                  ! not explicitly included in the barotropic equation, m s-2.
    v_accel_bt, & ! The difference between the meridional acceleration from the
                  ! barotropic calculation and BT_force_v, in m s-2.
    vhbt, &       ! The meridional barotropic thickness fluxes, in H m2 s-1.
    vhbt0, &      ! The difference between the sum of the layer meridional
                  ! thickness fluxes and the barotropic thickness flux using
                  ! the same velocities, in H m2 s-1.
    vbt_old, &    ! The starting value of vbt in a barotropic step, in m s-1.
    vbt_first, &  ! The starting value of ubt in a series of barotropic steps, in m s-1.
    vbt_sum, &    ! The sum of vbt over the time steps, in m s-1.
    vhbt_sum, &   ! The sum of vhbt over the time steps, in H m2 s-1.
    vbt_wtd, &    ! A weighted sum used to find the filtered final vbt, in m s-1.
    vbt_trans, &  ! The latest value of vbt used for a transport, in m s-1.
    cor_v, &      ! The meridional Coriolis acceleration, in m s-2.
    cor_ref_v, &  ! The meridional barotropic Coriolis acceleration due
                  ! to the reference velocities, in m s-2.
    pfv, &        ! The meridional pressure force acceleration, in m s-2.
    pfv_bt_sum, & ! The summed meridional barotropic pressure gradient force,
                  ! in m s-2.
    corv_bt_sum, & ! The summed meridional barotropic Coriolis acceleration,
                  ! in m s-2.
    dcor_v, &     ! A simply averaged depth at v points, in m.
    datv          ! Basin depth at v-velocity grid points times the x-grid
                  ! spacing, in H m.
  real, target, dimension(SZIW_(CS),SZJW_(CS)) :: &
    eta, &        ! The barotropic free surface height anomaly or column mass
                  ! anomaly, in H (m or kg m-2)
    eta_pred      ! A predictor value of eta, in H (m or kg m-2) like eta.
  real, pointer, dimension(:,:) :: &
    eta_PF_BT     ! A pointer to the eta array (either eta or eta_pred) that
                  ! determines the barotropic pressure force, in H (m or kg m-2)
  real, dimension(SZIW_(CS),SZJW_(CS)) :: &
    eta_sum, &    ! eta summed across the timesteps, in m or kg m-2.
    eta_wtd, &    ! A weighted estimate used to calculate eta_out, in m or kg m-2.
    eta_PF, &     ! A local copy of the 2-D eta field (either SSH anomaly or
                  ! column mass anomaly) that was used to calculate the input
                  ! pressure gradient accelerations, in m or kg m-2.
    eta_pf_1, &   ! The initial value of eta_PF, when interp_eta_PF is
                  ! true, in m or kg m-2.
    d_eta_pf, &   ! The change in eta_PF over the barotropic time stepping when
                  ! interp_eta_PF is true, in m or kg m-2.
    gtot_e, &     ! gtot_X is the effective total reduced gravity used to relate
    gtot_w, &     ! free surface height deviations to pressure forces (including
    gtot_n, &     ! GFS and baroclinic  contributions) in the barotropic momentum
    gtot_s, &     ! equations half a grid-point in the X-direction (X is N, S,
                  ! E, or W) from the thickness point. gtot_X has units of m2 H-1 s-2.
                  ! (See Hallberg, J Comp Phys 1997 for a discussion.)
    eta_src, &    ! The source of eta per barotropic timestep, in m or kg m-2.
    dyn_coef_eta, & ! The coefficient relating the changes in eta to the
                  ! dynamic surface pressure under rigid ice, in m2 s-2 H-1.
    p_surf_dyn    ! A dynamic surface pressure under rigid ice, in m2 s-2.
  type(local_bt_cont_u_type), dimension(SZIBW_(CS),SZJW_(CS)) :: &
    BTCL_u        ! A repackaged version of the u-point information in BT_cont.
  type(local_bt_cont_v_type), dimension(SZIW_(CS),SZJBW_(CS)) :: &
    BTCL_v        ! A repackaged version of the v-point information in BT_cont.
  ! End of wide-sized variables.

  real, dimension(SZIBW_(CS),SZJW_(CS)) :: &
    ubt_prev, uhbt_prev, ubt_sum_prev, uhbt_sum_prev, ubt_wtd_prev  ! for OBC
  real, dimension(SZIW_(CS),SZJBW_(CS)) :: &
    vbt_prev, vhbt_prev, vbt_sum_prev, vhbt_sum_prev, vbt_wtd_prev  ! for OBC

  real :: I_Rho0      ! The inverse of the mean density (Rho0), in m3 kg-1.
  real :: visc_rem    ! A work variable that may equal visc_rem_[uv].  Nondim.
  real :: vel_prev    ! The previous velocity in m s-1.
  real :: dtbt        ! The barotropic time step in s.
  real :: bebt        ! A copy of CS%bebt.
  real :: be_proj     ! The fractional amount by which velocities are projected
                      ! when project_velocity is true. For now be_proj is set
                      ! to equal bebt, as they have similar roles and meanings.
  real :: Idt         ! The inverse of dt, in s-1.
  real :: det_de      ! The partial derivative due to self-attraction and loading
                      ! of the reference geopotential with the sea surface height.
                      ! This is typically ~0.09 or less.
  real :: dgeo_de     ! The constant of proportionality between geopotential and
                      ! sea surface height.  It is a nondimensional number of
                      ! order 1.  For stability, this may be made larger
                      ! than physical problem would suggest.
  real :: Instep      ! The inverse of the number of barotropic time steps
                      ! to take.
  real :: wt_end      ! The weighting of the final value of eta_PF, ND.
  integer :: nstep    ! The number of barotropic time steps to take.
  type(time_type) :: &
    time_bt_start, &  ! The starting time of the barotropic steps.
    time_step_end, &  ! The end time of a barotropic step.
    time_end_in       ! The end time for diagnostics when this routine started.
  real :: time_int_in ! The diagnostics' time interval when this routine started.
  logical :: do_hifreq_output  ! If true, output occurs every barotropic step.
  logical :: use_BT_cont, do_ave, find_etaav, find_PF, find_Cor
  logical :: ice_is_rigid, nonblock_setup, interp_eta_PF
  logical :: project_velocity, add_uh0

  real :: dyn_coef_max ! The maximum stable value of dyn_coef_eta, in m2 s-2 H-1.
  real :: ice_strength = 0.0  ! The effective strength of the ice in m s-2.
  real :: Idt_max2    ! The squared inverse of the local maximum stable
                      ! barotropic time step, in s-2.
  real :: H_min_dyn   ! The minimum depth to use in limiting the size of the
                      ! dynamic surface pressure for stability, in H.
  real :: H_eff_dx2   ! The effective total thickness divided by the grid spacing
                      ! squared, in H m-2.
  real :: vel_tmp     ! A temporary velocity, in m s-1.
  real :: u_max_cor, v_max_cor ! The maximum corrective velocities, in m s-1.
  real :: Htot        ! The total thickness, in units of H.
  real :: eta_cor_max ! The maximum fluid that can be added as a correction to eta, in H.

  real, allocatable, dimension(:) :: wt_vel, wt_eta, wt_accel, wt_trans, wt_accel2
  real :: sum_wt_vel, sum_wt_eta, sum_wt_accel, sum_wt_trans
  real :: I_sum_wt_vel, I_sum_wt_eta, I_sum_wt_accel, I_sum_wt_trans
  real :: dt_filt     ! The half-width of the barotropic filter, in s.
  real :: trans_wt1, trans_wt2 ! weight used to compute ubt_trans and vbt_trans
  integer :: nfilter

  logical :: apply_OBCs, apply_OBC_flather, apply_OBC_open
  type(memory_size_type) :: MS
  character(len=200) :: mesg
  integer :: isv, iev, jsv, jev ! The valid array size at the end of a step.
  integer :: stencil  ! The stencil size of the algorithm, often 1 or 2.
  integer :: isvf, ievf, jsvf, jevf, num_cycles
  integer :: i, j, k, n
  integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  integer :: ioff, joff

  if (.not.associated(cs)) call mom_error(fatal, &
      "btstep: Module MOM_barotropic must be initialized before it is used.")
  if (.not.cs%split) return
  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
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  ms%isdw = cs%isdw ; ms%iedw = cs%iedw ; ms%jsdw = cs%jsdw ; ms%jedw = cs%jedw
  idt = 1.0 / dt

  use_bt_cont = .false.
  if (present(bt_cont)) use_bt_cont = (associated(bt_cont))

  interp_eta_pf = .false.
  if (present(eta_pf_start)) interp_eta_pf = (associated(eta_pf_start))

  project_velocity = cs%BT_project_velocity

  ! Figure out the fullest arrays that could be updated.
  stencil = 1
  if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
      (cs%Nonlin_cont_update_period > 0)) stencil = 2

  num_cycles = 1
  if (cs%use_wide_halos) &
    num_cycles = min((is-cs%isdw) / stencil, (js-cs%jsdw) / stencil)
  isvf = is - (num_cycles-1)*stencil ; ievf = ie + (num_cycles-1)*stencil
  jsvf = js - (num_cycles-1)*stencil ; jevf = je + (num_cycles-1)*stencil

  do_ave = query_averaging_enabled(cs%diag)
  find_etaav = present(etaav)
  find_pf = (do_ave .and. ((cs%id_PFu_bt > 0) .or. (cs%id_PFv_bt > 0)))
  find_cor = (do_ave .and. ((cs%id_Coru_bt > 0) .or. (cs%id_Corv_bt > 0)))

  add_uh0 = .false.
  if (present(uh0)) add_uh0 = associated(uh0)
  if (add_uh0 .and. .not.(present(vh0) .and. present(u_uh0) .and. &
                          present(v_vh0))) call mom_error(fatal, &
      "btstep: vh0, u_uh0, and v_vh0 must be present if uh0 is used.")
  if (add_uh0 .and. .not.(associated(vh0) .and. associated(u_uh0) .and. &
                          associated(v_vh0))) call mom_error(fatal, &
      "btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.")

  ! This can be changed to try to optimize the performance.
  nonblock_setup = g%nonblocking_updates

  if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)

  apply_obcs = .false. ; cs%BT_OBC%apply_u_OBCs = .false. ; cs%BT_OBC%apply_v_OBCs = .false.
  apply_obc_open = .false.
  apply_obc_flather = .false.
  if (present(obc)) then ; if (associated(obc)) then
    cs%BT_OBC%apply_u_OBCs = obc%open_u_BCs_exist_globally .or. obc%specified_u_BCs_exist_globally
    cs%BT_OBC%apply_v_OBCs = obc%open_v_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally
    apply_obc_flather = obc%Flather_u_BCs_exist_globally .or. obc%Flather_v_BCs_exist_globally
    apply_obc_open = obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally
    apply_obcs = obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
           apply_obc_flather .or. apply_obc_open
    if (.not.apply_obc_flather .and. obc%oblique_BCs_exist_globally) stencil = 2

    if (apply_obc_flather .and. .not.gv%Boussinesq) call mom_error(fatal, &
      "btstep: Flather open boundary conditions have not yet been "// &
      "implemented for a non-Boussinesq model.")
  endif ; endif

  nstep = ceiling(dt/cs%dtbt - 0.0001)
  if (is_root_pe() .and. (nstep /= cs%nstep_last)) then
    write(mesg,'("btstep is using a dynamic barotropic timestep of ", ES12.6, &
               & " seconds, max ", ES12.6, ".")') (dt/nstep), cs%dtbt_max
    call mom_mesg(mesg, 3)
  endif
  cs%nstep_last = nstep

  ! Set the actual barotropic time step.
  instep = 1.0 / real(nstep)
  dtbt = dt * instep
  bebt = cs%bebt
  be_proj = cs%bebt
  i_rho0 = 1.0/gv%Rho0

  !--- setup the weight when computing vbt_trans and ubt_trans
  if (project_velocity) then
    trans_wt1 = (1.0 + be_proj); trans_wt2 = -be_proj
  else
    trans_wt1 = bebt ;           trans_wt2 = (1.0-bebt)
  endif

  do_hifreq_output = .false.
  if ((cs%id_ubt_hifreq > 0) .or. (cs%id_vbt_hifreq > 0) .or. &
      (cs%id_eta_hifreq > 0) .or. (cs%id_eta_pred_hifreq > 0) .or. &
      (cs%id_uhbt_hifreq > 0) .or. (cs%id_vhbt_hifreq > 0)) then
    do_hifreq_output = query_averaging_enabled(cs%diag, time_int_in, time_end_in)
    if (do_hifreq_output) &
      time_bt_start = time_end_in - set_time(int(floor(dt+0.5)))
  endif

!--- begin setup for group halo update
  if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
  if (.not. cs%linearized_BT_PV) then
    call create_group_pass(cs%pass_q_DCor, q, cs%BT_Domain, to_all, position=corner)
    call create_group_pass(cs%pass_q_DCor, dcor_u, dcor_v, cs%BT_Domain, &
         to_all+scalar_pair)
  endif
  if ((isq > is-1) .or. (jsq > js-1)) &
    call create_group_pass(cs%pass_tmp_uv, tmp_u, tmp_v, g%Domain)
  call create_group_pass(cs%pass_gtot, gtot_e, gtot_n, cs%BT_Domain, &
       to_all+scalar_pair, agrid)
  call create_group_pass(cs%pass_gtot, gtot_w, gtot_s, cs%BT_Domain, &
       to_all+scalar_pair, agrid)

  if (cs%dynamic_psurf) &
    call create_group_pass(cs%pass_eta_bt_rem, dyn_coef_eta, cs%BT_Domain)
  if (interp_eta_pf) then
    call create_group_pass(cs%pass_eta_bt_rem, eta_pf_1, cs%BT_Domain)
    call create_group_pass(cs%pass_eta_bt_rem, d_eta_pf, cs%BT_Domain)
  else
    call create_group_pass(cs%pass_eta_bt_rem, eta_pf, cs%BT_Domain)
  endif
  call create_group_pass(cs%pass_eta_bt_rem, eta_src, cs%BT_Domain)
  ! The following halo update is not needed without wide halos.  RWH
  if (ievf > ie) call create_group_pass(cs%pass_eta_bt_rem, bt_rem_u, bt_rem_v, &
                      cs%BT_Domain, to_all+scalar_pair)
  ! The following halo update is not needed without wide halos.  RWH
  if (((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) .or. (isq <= is-1) .or. (jsq <= js-1)) &
    call create_group_pass(cs%pass_force_hbt0_Cor_ref, bt_force_u, bt_force_v, cs%BT_Domain)
  if (add_uh0) call create_group_pass(cs%pass_force_hbt0_Cor_ref, uhbt0, vhbt0, cs%BT_Domain)
  call create_group_pass(cs%pass_force_hbt0_Cor_ref, cor_ref_u, cor_ref_v, cs%BT_Domain)
  if (.not. use_bt_cont) then
    call create_group_pass(cs%pass_Dat_uv, datu, datv, cs%BT_Domain, to_all+scalar_pair)
  endif
  call create_group_pass(cs%pass_eta_ubt, eta, cs%BT_Domain)
  call create_group_pass(cs%pass_eta_ubt, ubt, vbt, cs%BT_Domain)

  call create_group_pass(cs%pass_ubt_Cor, ubt_cor, vbt_cor, g%Domain)
  ! These passes occur at the end of the routine, as data is being readied to
  ! share with the main part of the MOM6 code.
  if (find_etaav) then
    call create_group_pass(cs%pass_etaav, etaav, g%Domain)
  endif
  call create_group_pass(cs%pass_e_anom, e_anom, g%Domain)
  call create_group_pass(cs%pass_ubta_uhbta, cs%ubtav, cs%vbtav, g%Domain)
  call create_group_pass(cs%pass_ubta_uhbta, uhbtav, vhbtav, g%Domain)

  if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
!--- end setup for group halo update

!   Calculate the constant coefficients for the Coriolis force terms in the
! barotropic momentum equations.  This has to be done quite early to start
! the halo update that needs to be completed before the next calculations.
  if (cs%linearized_BT_PV) then
    !$OMP parallel do default(shared)
    do j=jsvf-2,jevf+1 ; do i=isvf-2,ievf+1
      q(i,j) = cs%q_D(i,j)
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=jsvf-1,jevf+1 ; do i=isvf-2,ievf+1
      dcor_u(i,j) = cs%D_u_Cor(i,j)
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=jsvf-2,jevf+1 ; do i=isvf-1,ievf+1
      dcor_v(i,j) = cs%D_v_Cor(i,j)
    enddo ; enddo
  else
    q(:,:) = 0.0 ; dcor_u(:,:) = 0.0 ; dcor_v(:,:) = 0.0
    !  This option has not yet been written properly.
    !  ### bathyT here should be replaced with bathyT+eta(Bous) or eta(non-Bous).
    !$OMP parallel do default(shared)
    do j=js,je ; do i=is-1,ie
      dcor_u(i,j) = 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=js-1,je ; do i=is,ie
      dcor_v(i,j) = 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j))
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=js-1,je ; do i=is-1,ie
      q(i,j) = 0.25 * g%CoriolisBu(i,j) * &
           ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
           ((g%areaT(i,j) * g%bathyT(i,j) + g%areaT(i+1,j+1) * g%bathyT(i+1,j+1)) + &
            (g%areaT(i+1,j) * g%bathyT(i+1,j) + g%areaT(i,j+1) * g%bathyT(i,j+1)))
    enddo ; enddo

    ! With very wide halos, q and D need to be calculated on the available data
    ! domain and then updated onto the full computational domain.
    ! These calculations can be done almost immediately, but the halo updates
    ! must be done before the [abcd]mer and [abcd]zon are calculated.
    if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
    if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
    if (nonblock_setup) then
      call start_group_pass(cs%pass_q_DCor, cs%BT_Domain)
    else
      call do_group_pass(cs%pass_q_DCor, cs%BT_Domain)
    endif
    if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
    if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)
  endif

  ! Zero out various wide-halo arrays.
  !$OMP parallel do default(shared)
  do j=cs%jsdw,cs%jedw ; do i=cs%isdw,cs%iedw
    gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
    gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
    eta(i,j) = 0.0
    eta_pf(i,j) = 0.0
    if (interp_eta_pf) then
      eta_pf_1(i,j) = 0.0 ; d_eta_pf(i,j) = 0.0
    endif
    p_surf_dyn(i,j) = 0.0
    if (cs%dynamic_psurf) dyn_coef_eta(i,j) = 0.0
  enddo ; enddo
  !   The halo regions of various arrays need to be initialized to
  ! non-NaNs in case the neighboring domains are not part of the ocean.
  ! Otherwise a halo update later on fills in the correct values.
  !$OMP parallel do default(shared)
  do j=cs%jsdw,cs%jedw ; do i=cs%isdw-1,cs%iedw
    cor_ref_u(i,j) = 0.0 ; bt_force_u(i,j) = 0.0 ; ubt(i,j) = 0.0
    datu(i,j) = 0.0 ; bt_rem_u(i,j) = 0.0 ; uhbt0(i,j) = 0.0
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=cs%jsdw-1,cs%jedw ; do i=cs%isdw,cs%iedw
    cor_ref_v(i,j) = 0.0 ; bt_force_v(i,j) = 0.0 ; vbt(i,j) = 0.0
    datv(i,j) = 0.0 ; bt_rem_v(i,j) = 0.0 ; vhbt0(i,j) = 0.0
  enddo ; enddo

  ! Copy input arrays into their wide-halo counterparts.
  if (interp_eta_pf) then
    !$OMP parallel do default(shared)
    do j=g%jsd,g%jed ; do i=g%isd,g%ied ! Was "do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1" but doing so breaks OBC. Not sure why?
      eta(i,j) = eta_in(i,j)
      eta_pf_1(i,j) = eta_pf_start(i,j)
      d_eta_pf(i,j) = eta_pf_in(i,j) - eta_pf_start(i,j)
    enddo ; enddo
  else
    !$OMP parallel do default(shared)
    do j=g%jsd,g%jed ; do i=g%isd,g%ied !: Was "do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1" but doing so breaks OBC. Not sure why?
      eta(i,j) = eta_in(i,j)
      eta_pf(i,j) = eta_pf_in(i,j)
    enddo ; enddo
  endif

  !$OMP parallel do default(shared) private(visc_rem)
  do k=1,nz ; do j=js,je ; do i=is-1,ie
    ! rem needs greater than visc_rem_u and 1-Instep/visc_rem_u.
    ! The 0.5 below is just for safety.
    if (visc_rem_u(i,j,k) <= 0.0) then ; visc_rem = 0.0
    elseif (visc_rem_u(i,j,k) >= 1.0) then ; visc_rem = 1.0
    elseif (visc_rem_u(i,j,k)**2 > visc_rem_u(i,j,k) - 0.5*instep) then
      visc_rem = visc_rem_u(i,j,k)
    else ; visc_rem = 1.0 - 0.5*instep/visc_rem_u(i,j,k) ; endif
    wt_u(i,j,k) = cs%frhatu(i,j,k) * visc_rem
  enddo ; enddo ; enddo
  !$OMP parallel do default(shared) private(visc_rem)
  do k=1,nz ; do j=js-1,je ; do i=is,ie
    ! rem needs greater than visc_rem_v and 1-Instep/visc_rem_v.
    if (visc_rem_v(i,j,k) <= 0.0) then ; visc_rem = 0.0
    elseif (visc_rem_v(i,j,k) >= 1.0) then ; visc_rem = 1.0
    elseif (visc_rem_v(i,j,k)**2 > visc_rem_v(i,j,k) - 0.5*instep) then
      visc_rem = visc_rem_v(i,j,k)
    else ; visc_rem = 1.0 - 0.5*instep/visc_rem_v(i,j,k) ; endif
    wt_v(i,j,k) = cs%frhatv(i,j,k) * visc_rem
  enddo ; enddo ; enddo

  !   Use u_Cor and v_Cor as the reference values for the Coriolis terms,
  ! including the viscous remnant.
  !$OMP parallel do default(shared)
  do j=js-1,je+1 ; do i=is-1,ie ; ubt_cor(i,j) = 0.0 ; enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js-1,je ; do i=is-1,ie+1 ; vbt_cor(i,j) = 0.0 ; enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js,je ; do k=1,nz ; do i=is-1,ie
    ubt_cor(i,j) = ubt_cor(i,j) + wt_u(i,j,k) * u_cor(i,j,k)
  enddo ; enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js-1,je ; do k=1,nz ; do i=is,ie
    vbt_cor(i,j) = vbt_cor(i,j) + wt_v(i,j,k) * v_cor(i,j,k)
  enddo ; enddo ; enddo

  ! The gtot arrays are the effective layer-weighted reduced gravities for
  ! accelerations across the various faces, with names for the relative
  ! locations of the faces to the pressure point.  They will have their halos
  ! updated later on.
  !$OMP parallel do default(shared)
  do j=js,je
    do k=1,nz ; do i=is-1,ie
      gtot_e(i,j)   = gtot_e(i,j)   + pbce(i,j,k)   * wt_u(i,j,k)
      gtot_w(i+1,j) = gtot_w(i+1,j) + pbce(i+1,j,k) * wt_u(i,j,k)
    enddo ; enddo
  enddo
  !$OMP parallel do default(shared)
  do j=js-1,je
     do k=1,nz ; do i=is,ie
      gtot_n(i,j)   = gtot_n(i,j)   + pbce(i,j,k)   * wt_v(i,j,k)
      gtot_s(i,j+1) = gtot_s(i,j+1) + pbce(i,j+1,k) * wt_v(i,j,k)
    enddo ; enddo
  enddo

  if (cs%tides) then
    call tidal_forcing_sensitivity(g, cs%tides_CSp, det_de)
    dgeo_de = 1.0 + det_de + cs%G_extra
  else
    dgeo_de = 1.0 + cs%G_extra
  endif

  if (nonblock_setup .and. .not.cs%linearized_BT_PV) then
    if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
    if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
    call complete_group_pass(cs%pass_q_DCor, cs%BT_Domain)
    if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
    if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)
  endif

  ! Calculate the open areas at the velocity points.
  ! The halo updates are needed before Datu is first used, either in set_up_BT_OBC or ubt_Cor.
  if (use_bt_cont) then
    call set_local_bt_cont_types(bt_cont, btcl_u, btcl_v, g, ms, cs%BT_Domain, 1+ievf-ie)
  else
    if (cs%Nonlinear_continuity) then
      call find_face_areas(datu, datv, g, gv, cs, ms, eta, 1)
    else
      call find_face_areas(datu, datv, g, gv, cs, ms, halo=1)
    endif
  endif

  ! Set up fields related to the open boundary conditions.
  if (apply_obcs) then
    call set_up_bt_obc(obc, eta, cs%BT_OBC, g, gv, ms, ievf-ie, use_bt_cont, &
                       datu, datv, btcl_u, btcl_v)
  endif

!   Here the vertical average accelerations due to the Coriolis, advective,
! pressure gradient and horizontal viscous terms in the layer momentum
! equations are calculated.  These will be used to determine the difference
! between the accelerations due to the average of the layer equations and the
! barotropic calculation.

  !$OMP parallel do default(shared)
  do j=js,je ; do i=is-1,ie
    ! ### IDatu here should be replaced with 1/D+eta(Bous) or 1/eta(non-Bous).
    ! ### although with BT_cont_types IDatu should be replaced by
    ! ###   CS%dy_Cu(I,j) / (d(uhbt)/du) (with appropriate bounds).
    bt_force_u(i,j) = fluxes%taux(i,j) * i_rho0*cs%IDatu(i,j)*visc_rem_u(i,j,1)
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js-1,je ; do i=is,ie
    ! ### IDatv here should be replaced with 1/D+eta(Bous) or 1/eta(non-Bous).
    ! ### although with BT_cont_types IDatv should be replaced by
    ! ###   CS%dx_Cv(I,j) / (d(vhbt)/dv) (with appropriate bounds).
    bt_force_v(i,j) = fluxes%tauy(i,j) * i_rho0*cs%IDatv(i,j)*visc_rem_v(i,j,1)
  enddo ; enddo
  if (present(taux_bot) .and. present(tauy_bot)) then
    if (associated(taux_bot) .and. associated(tauy_bot)) then
      !$OMP parallel do default(shared)
      do j=js,je ; do i=is-1,ie
        bt_force_u(i,j) = bt_force_u(i,j) - taux_bot(i,j) * i_rho0 * cs%IDatu(i,j)
      enddo ; enddo
      !$OMP parallel do default(shared)
      do j=js-1,je ; do i=is,ie
        bt_force_v(i,j) = bt_force_v(i,j) - tauy_bot(i,j) * i_rho0 * cs%IDatv(i,j)
      enddo ; enddo
    endif
  endif

  ! bc_accel_u & bc_accel_v are only available on the potentially
  ! non-symmetric computational domain.
  !$OMP parallel do default(shared)
  do j=js,je ; do k=1,nz ; do i=isq,ieq
    bt_force_u(i,j) = bt_force_u(i,j) + wt_u(i,j,k) * bc_accel_u(i,j,k)
  enddo ; enddo ; enddo
  !$OMP parallel do default(shared)
  do j=jsq,jeq ; do k=1,nz ; do i=is,ie
    bt_force_v(i,j) = bt_force_v(i,j) + wt_v(i,j,k) * bc_accel_v(i,j,k)
  enddo ; enddo ; enddo

  ! Determine the difference between the sum of the layer fluxes and the
  ! barotropic fluxes found from the same input velocities.
  if (add_uh0) then
    !$OMP parallel do default(shared)
    do j=js,je ; do i=is-1,ie ; uhbt(i,j) = 0.0 ; ubt(i,j) = 0.0 ; enddo ; enddo
    !$OMP parallel do default(shared)
    do j=js-1,je ; do i=is,ie ; vhbt(i,j) = 0.0 ; vbt(i,j) = 0.0 ; enddo ; enddo
    if (cs%visc_rem_u_uh0) then
      !$OMP parallel do default(shared)
      do j=js,je ; do k=1,nz ; do i=is-1,ie
        uhbt(i,j) = uhbt(i,j) + uh0(i,j,k)
        ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_uh0(i,j,k)
      enddo ; enddo ; enddo
      !$OMP parallel do default(shared)
      do j=js-1,je ; do k=1,nz ; do i=is,ie
        vhbt(i,j) = vhbt(i,j) + vh0(i,j,k)
        vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_vh0(i,j,k)
      enddo ; enddo ; enddo
    else
      !$OMP parallel do default(shared)
      do j=js,je ; do k=1,nz ; do i=is-1,ie
        uhbt(i,j) = uhbt(i,j) + uh0(i,j,k)
        ubt(i,j) = ubt(i,j) + cs%frhatu(i,j,k) * u_uh0(i,j,k)
      enddo ; enddo ; enddo
      !$OMP parallel do default(shared)
      do j=js-1,je ; do k=1,nz ; do i=is,ie
        vhbt(i,j) = vhbt(i,j) + vh0(i,j,k)
        vbt(i,j) = vbt(i,j) + cs%frhatv(i,j,k) * v_vh0(i,j,k)
      enddo ; enddo ; enddo
    endif
    if (use_bt_cont) then
      if (cs%adjust_BT_cont) then
        ! Use the additional input transports to broaden the fits
        ! over which the bt_cont_type applies.

        ! Fill in the halo data for ubt, vbt, uhbt, and vhbt.
        if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
        if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
        call pass_vector(ubt, vbt, cs%BT_Domain, complete=.false., halo=1+ievf-ie)
        call pass_vector(uhbt, vhbt, cs%BT_Domain, complete=.true., halo=1+ievf-ie)
        if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
        if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)

        call adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, btcl_u, btcl_v, &
                                        g, ms, 1+ievf-ie)
      endif
      !$OMP parallel do default(shared)
      do j=js,je ; do i=is-1,ie
        uhbt0(i,j) = uhbt(i,j) - find_uhbt(ubt(i,j),btcl_u(i,j))
      enddo ; enddo
      !$OMP parallel do default(shared)
      do j=js-1,je ; do i=is,ie
        vhbt0(i,j) = vhbt(i,j) - find_vhbt(vbt(i,j),btcl_v(i,j))
      enddo ; enddo
    else
      !$OMP parallel do default(shared)
      do j=js,je ; do i=is-1,ie
        uhbt0(i,j) = uhbt(i,j) - datu(i,j)*ubt(i,j)
      enddo ; enddo
      !$OMP parallel do default(shared)
      do j=js-1,je ; do i=is,ie
        vhbt0(i,j) = vhbt(i,j) - datv(i,j)*vbt(i,j)
      enddo ; enddo
    endif
  endif

! Calculate the initial barotropic velocities from the layer's velocities.
  !$OMP parallel do default(shared)
  do j=jsvf-1,jevf+1 ; do i=isvf-2,ievf+1
    ubt(i,j) = 0.0 ; uhbt(i,j) = 0.0 ; u_accel_bt(i,j) = 0.0
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=jsvf-2,jevf+1 ; do i=isvf-1,ievf+1
    vbt(i,j) = 0.0 ; vhbt(i,j) = 0.0 ; v_accel_bt(i,j) = 0.0
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js,je ; do k=1,nz ; do i=is-1,ie
    ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_in(i,j,k)
  enddo ; enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js-1,je ; do k=1,nz ; do i=is,ie
    vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_in(i,j,k)
  enddo ; enddo ;  enddo

  if (apply_obcs) then
    ubt_first(:,:) = ubt(:,:) ; vbt_first(:,:) = vbt(:,:)
  endif

  if (cs%gradual_BT_ICs) then
    !$OMP parallel do default(shared)
    do j=js,je ; do i=is-1,ie
      bt_force_u(i,j) = bt_force_u(i,j) + (ubt(i,j) - cs%ubt_IC(i,j)) * idt
      ubt(i,j) = cs%ubt_IC(i,j)
    enddo ; enddo
    !$OMP parallel do default(shared)
    do j=js-1,je ; do i=is,ie
      bt_force_v(i,j) = bt_force_v(i,j) + (vbt(i,j) - cs%vbt_IC(i,j)) * idt
      vbt(i,j) = cs%vbt_IC(i,j)
    enddo ; enddo
  endif

  if ((isq > is-1) .or. (jsq > js-1)) then
    ! Non-symmetric memory is being used, so the edge values need to be
    ! filled in with a halo update of a non-symmetric array.
    if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
    if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
    tmp_u(:,:) = 0.0 ; tmp_v(:,:) = 0.0
    do j=js,je ; do i=isq,ieq ; tmp_u(i,j) = bt_force_u(i,j) ; enddo ; enddo
    do j=jsq,jeq ; do i=is,ie ; tmp_v(i,j) = bt_force_v(i,j) ; enddo ; enddo
    if (nonblock_setup) then
      call start_group_pass(cs%pass_tmp_uv, g%Domain)
    else
      call do_group_pass(cs%pass_tmp_uv, g%Domain)
      do j=jsd,jed ; do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ; enddo ; enddo
      do j=jsdb,jedb ; do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ; enddo ; enddo
    endif
    if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
    if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)
  endif

  if (nonblock_setup) then
    if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
    if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
    call start_group_pass(cs%pass_gtot, cs%BT_Domain)
    call start_group_pass(cs%pass_ubt_Cor, g%Domain)
    if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
    if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)
  endif

  ! Determine the weighted Coriolis parameters for the neighboring velocities.
  !$OMP parallel do default(shared)
  do j=jsvf-1,jevf ; do i=isvf-1,ievf+1
    if (cs%Sadourny) then
      amer(i-1,j) = dcor_u(i-1,j) * q(i-1,j)
      bmer(i,j) = dcor_u(i,j) * q(i,j)
      cmer(i,j+1) = dcor_u(i,j+1) * q(i,j)
      dmer(i-1,j+1) = dcor_u(i-1,j+1) * q(i-1,j)
    else
      amer(i-1,j) = dcor_u(i-1,j) * &
                    ((q(i,j) + q(i-1,j-1)) + q(i-1,j)) / 3.0
      bmer(i,j) = dcor_u(i,j) * &
                  (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
      cmer(i,j+1) = dcor_u(i,j+1) * &
                    (q(i,j) + (q(i-1,j) + q(i,j+1))) / 3.0
      dmer(i-1,j+1) = dcor_u(i-1,j+1) * &
                      ((q(i,j) + q(i-1,j+1)) + q(i-1,j)) / 3.0
    endif
  enddo ; enddo

  !$OMP parallel do default(shared)
  do j=jsvf-1,jevf+1 ; do i=isvf-1,ievf
    if (cs%Sadourny) then
      azon(i,j) = dcor_v(i+1,j) * q(i,j)
      bzon(i,j) = dcor_v(i,j) * q(i,j)
      czon(i,j) = dcor_v(i,j-1) * q(i,j-1)
      dzon(i,j) = dcor_v(i+1,j-1) * q(i,j-1)
    else
      azon(i,j) = dcor_v(i+1,j) * &
                  (q(i,j) + (q(i+1,j) + q(i,j-1))) / 3.0
      bzon(i,j) = dcor_v(i,j) * &
                  (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
      czon(i,j) = dcor_v(i,j-1) * &
                  ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) / 3.0
      dzon(i,j) = dcor_v(i+1,j-1) * &
                  ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) / 3.0
    endif
  enddo ; enddo

! Complete the previously initiated message passing.
  if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
  if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
  if (nonblock_setup) then
    if ((isq > is-1) .or. (jsq > js-1)) then
      call complete_group_pass(cs%pass_tmp_uv, g%Domain)
      do j=jsd,jed ; do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ; enddo ; enddo
      do j=jsdb,jedb ; do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ; enddo ; enddo
    endif
    call complete_group_pass(cs%pass_gtot, cs%BT_Domain)
    call complete_group_pass(cs%pass_ubt_Cor, g%Domain)
  else
    call do_group_pass(cs%pass_gtot, cs%BT_Domain)
    call do_group_pass(cs%pass_ubt_Cor, g%Domain)
  endif
  ! The various elements of gtot are positive definite but directional, so use
  ! the polarity arrays to sort out when the directions have shifted.
  do j=jsvf-1,jevf+1 ; do i=isvf-1,ievf+1
    if (cs%ua_polarity(i,j) < 0.0) call swap(gtot_e(i,j), gtot_w(i,j))
    if (cs%va_polarity(i,j) < 0.0) call swap(gtot_n(i,j), gtot_s(i,j))
  enddo ; enddo

  !$OMP parallel do default(shared)
  do j=js,je ; do i=is-1,ie
    cor_ref_u(i,j) =  &
        ((azon(i,j) * vbt_cor(i+1,j) + czon(i,j) * vbt_cor(i  ,j-1)) + &
         (bzon(i,j) * vbt_cor(i  ,j) + dzon(i,j) * vbt_cor(i+1,j-1)))
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js-1,je ; do i=is,ie
    cor_ref_v(i,j) = -1.0 * &
        ((amer(i-1,j) * ubt_cor(i-1,j) + cmer(i  ,j+1) * ubt_cor(i  ,j+1)) + &
         (bmer(i  ,j) * ubt_cor(i  ,j) + dmer(i-1,j+1) * ubt_cor(i-1,j+1)))
  enddo ; enddo

  ! Now start new halo updates.
  if (nonblock_setup) then
    if (.not.use_bt_cont) &
      call start_group_pass(cs%pass_Dat_uv, cs%BT_Domain)

    ! The following halo update is not needed without wide halos.  RWH
    call start_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
  endif
  if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
  if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)
!$OMP parallel default(none) shared(is,ie,js,je,nz,av_rem_u,av_rem_v,CS,visc_rem_u, &
!$OMP                               visc_rem_v,bt_rem_u,G,GV,nstep,bt_rem_v,Instep, &
!$OMP                               find_etaav,jsvf,jevf,isvf,ievf,eta_sum,eta_wtd, &
!$OMP                               ubt_sum,uhbt_sum,PFu_bt_sum,Coru_bt_sum,ubt_wtd,&
!$OMP                               ubt_trans,vbt_sum,vhbt_sum,PFv_bt_sum,          &
!$OMP                               Corv_bt_sum,vbt_wtd,vbt_trans,eta_src,dt,dtbt,  &
!$OMP                               use_BT_Cont,BTCL_u,uhbt0,BTCL_v,vhbt0,eta,Idt)  &
!$OMP                       private(u_max_cor,v_max_cor,eta_cor_max,Htot)
!$OMP do
  do j=js-1,je+1 ; do i=is-1,ie ; av_rem_u(i,j) = 0.0 ; enddo ; enddo
!$OMP do
  do j=js-1,je ; do i=is-1,ie+1 ; av_rem_v(i,j) = 0.0 ; enddo ; enddo
!$OMP do
  do j=js,je ; do k=1,nz ; do i=is-1,ie
    av_rem_u(i,j) = av_rem_u(i,j) + cs%frhatu(i,j,k) * visc_rem_u(i,j,k)
  enddo ; enddo ; enddo
!$OMP do
  do j=js-1,je ; do k=1,nz ; do i=is,ie
    av_rem_v(i,j) = av_rem_v(i,j) + cs%frhatv(i,j,k) * visc_rem_v(i,j,k)
  enddo ; enddo ; enddo
  if (cs%strong_drag) then
!$OMP do
    do j=js,je ; do i=is-1,ie
      bt_rem_u(i,j) = g%mask2dCu(i,j) * &
         ((nstep * av_rem_u(i,j)) / (1.0 + (nstep-1)*av_rem_u(i,j)))
    enddo ; enddo
!$OMP do
    do j=js-1,je ; do i=is,ie
      bt_rem_v(i,j) = g%mask2dCv(i,j) * &
         ((nstep * av_rem_v(i,j)) / (1.0 + (nstep-1)*av_rem_v(i,j)))
    enddo ; enddo
  else
!$OMP do
    do j=js,je ; do i=is-1,ie
      bt_rem_u(i,j) = 0.0
      if (g%mask2dCu(i,j) * av_rem_u(i,j) > 0.0) &
        bt_rem_u(i,j) = g%mask2dCu(i,j) * (av_rem_u(i,j)**instep)
    enddo ; enddo
!$OMP do
    do j=js-1,je ; do i=is,ie
      bt_rem_v(i,j) = 0.0
      if (g%mask2dCv(i,j) * av_rem_v(i,j) > 0.0) &
        bt_rem_v(i,j) = g%mask2dCv(i,j) * (av_rem_v(i,j)**instep)
    enddo ; enddo
  endif
  ! Zero out the arrays for various time-averaged quantities.
  if (find_etaav) then
!$OMP do
    do j=jsvf-1,jevf+1 ; do i=isvf-1,ievf+1
      eta_sum(i,j) = 0.0 ; eta_wtd(i,j) = 0.0
    enddo ; enddo
  else
!$OMP do
    do j=jsvf-1,jevf+1 ; do i=isvf-1,ievf+1
      eta_wtd(i,j) = 0.0
    enddo ; enddo
   endif
!$OMP do
  do j=jsvf-1,jevf+1 ; do i=isvf-1,ievf
    ubt_sum(i,j) = 0.0 ; uhbt_sum(i,j) = 0.0
    pfu_bt_sum(i,j) = 0.0 ; coru_bt_sum(i,j) = 0.0
    ubt_wtd(i,j) = 0.0 ; ubt_trans(i,j) = 0.0
  enddo ; enddo
!$OMP do
  do j=jsvf-1,jevf ; do i=isvf-1,ievf+1
    vbt_sum(i,j) = 0.0 ; vhbt_sum(i,j) = 0.0
    pfv_bt_sum(i,j) = 0.0 ; corv_bt_sum(i,j) = 0.0
    vbt_wtd(i,j) = 0.0 ; vbt_trans(i,j) = 0.0
  enddo ; enddo

  ! Set the mass source, after first initializing the halos to 0.
!$OMP do
  do j=jsvf-1,jevf+1; do i=isvf-1,ievf+1 ; eta_src(i,j) = 0.0 ; enddo ; enddo
  if (cs%bound_BT_corr) then ; if (use_bt_cont .and. cs%BT_cont_bounds) then
    do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
      if (cs%eta_cor(i,j) > 0.0) then
        !   Limit the source (outward) correction to be a fraction the mass that
        ! can be transported out of the cell by velocities with a CFL number of
        ! CFL_cor.
        u_max_cor = g%dxT(i,j) * (cs%maxCFL_BT_cont*idt)
        v_max_cor = g%dyT(i,j) * (cs%maxCFL_BT_cont*idt)
        eta_cor_max = dt * (cs%IareaT(i,j) * &
                   (((find_uhbt(u_max_cor,btcl_u(i,j)) + uhbt0(i,j)) - &
                     (find_uhbt(-u_max_cor,btcl_u(i-1,j)) + uhbt0(i-1,j))) + &
                    ((find_vhbt(v_max_cor,btcl_v(i,j)) + vhbt0(i,j)) - &
                     (find_vhbt(-v_max_cor,btcl_v(i,j-1)) + vhbt0(i,j-1))) ) - &
                   cs%eta_source(i,j))
        cs%eta_cor(i,j) = min(cs%eta_cor(i,j), max(0.0, eta_cor_max))
      else
        ! Limit the sink (inward) correction to the amount of mass that is already
        ! inside the cell, plus any mass added by eta_source.
        htot = eta(i,j)
        if (gv%Boussinesq) htot = cs%bathyT(i,j)*gv%m_to_H + eta(i,j)

        cs%eta_cor(i,j) = max(cs%eta_cor(i,j), -max(0.0,htot + dt*cs%eta_source(i,j)))
      endif
    endif ; enddo ; enddo
  else ; do j=js,je ; do i=is,ie
    if (abs(cs%eta_cor(i,j)) > dt*cs%eta_cor_bound(i,j)) &
      cs%eta_cor(i,j) = sign(dt*cs%eta_cor_bound(i,j),cs%eta_cor(i,j))
  enddo ; enddo ; endif ; endif
!$OMP do
  do j=js,je ; do i=is,ie
    eta_src(i,j) = g%mask2dT(i,j) * (instep * cs%eta_cor(i,j) + dtbt * cs%eta_source(i,j))
  enddo ; enddo
!$OMP end parallel

  if (cs%dynamic_psurf) then
    ice_is_rigid = (associated(fluxes%rigidity_ice_u) .and. &
                    associated(fluxes%rigidity_ice_v))
    h_min_dyn = gv%m_to_H * cs%Dmin_dyn_psurf
    if (ice_is_rigid .and. use_bt_cont) &
      call bt_cont_to_face_areas(bt_cont, datu, datv, g, ms, 0, .true.)
    if (ice_is_rigid) then
!$OMP parallel do default(none) shared(is,ie,js,je,dgeo_de,bebt,G,GV,gtot_E,Datu, &
!$OMP                                  gtot_W,gtot_N,gtot_S,Datv,H_min_dyn,     &
!$OMP                                  CS,dtbt,fluxes,dyn_coef_eta )            &
!$OMP                          private(Idt_max2,H_eff_dx2,dyn_coef_max,ice_strength)
      do j=js,je ; do i=is,ie
      ! First determine the maximum stable value for dyn_coef_eta.

      !   This estimate of the maximum stable time step is pretty accurate for
      ! gravity waves, but it is a conservative estimate since it ignores the
      ! stabilizing effect of the bottom drag.
      idt_max2 = 0.5 * (dgeo_de * (1.0 + 2.0*bebt)) * (g%IareaT(i,j) * &
            ((gtot_e(i,j) * (datu(i,j)*g%IdxCu(i,j)) + &
              gtot_w(i,j) * (datu(i-1,j)*g%IdxCu(i-1,j))) + &
             (gtot_n(i,j) * (datv(i,j)*g%IdyCv(i,j)) + &
              gtot_s(i,j) * (datv(i,j-1)*g%IdyCv(i,j-1)))) + &
            ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
             (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)))
      h_eff_dx2 = max(h_min_dyn * (g%IdxT(i,j)**2 + g%IdyT(i,j)**2), &
                      g%IareaT(i,j) * &
                        ((datu(i,j)*g%IdxCu(i,j) + datu(i-1,j)*g%IdxCu(i-1,j)) + &
                         (datv(i,j)*g%IdyCv(i,j) + datv(i,j-1)*g%IdyCv(i,j-1)) ) )
      dyn_coef_max = cs%const_dyn_psurf * max(0.0, 1.0 - dtbt**2 * idt_max2) / &
                     (dtbt**2 * h_eff_dx2)

      ! ice_strength has units of m s-2. rigidity_ice_[uv] has units of m3 s-1.
      ice_strength = ((fluxes%rigidity_ice_u(i,j) + fluxes%rigidity_ice_u(i-1,j)) + &
                      (fluxes%rigidity_ice_v(i,j) + fluxes%rigidity_ice_v(i,j-1))) / &
                      (cs%ice_strength_length**2 * dtbt)

      ! Units of dyn_coef: m2 s-2 H-1
      dyn_coef_eta(i,j) = min(dyn_coef_max, ice_strength * gv%H_to_m)
    enddo ; enddo ; endif
  endif

  if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
  if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
  if (nonblock_setup) then
    call start_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
    ! The following halo update is not needed without wide halos.  RWH
  else
    call do_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
    if (.not.use_bt_cont) &
      call do_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
    call do_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
  endif
  if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
  if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)

  ! Complete all of the outstanding halo updates.
  if (nonblock_setup) then
    if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
    if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)

    if (.not.use_bt_cont) &  !### IS THIS OK HERE?
      call complete_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
    call complete_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
    call complete_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)

    if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
    if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)
  endif

  if (cs%debug) then
    call uvchksum("BT [uv]hbt", uhbt, vhbt, cs%debug_BT_HI, haloshift=0, &
                  scale=gv%H_to_m)
    call uvchksum("BT Initial [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=0)
    call hchksum(eta, "BT Initial eta", cs%debug_BT_HI, haloshift=0, scale=gv%H_to_m)
    call uvchksum("BT BT_force_[uv]", bt_force_u, bt_force_v, &
                  cs%debug_BT_HI, haloshift=0)
    if (interp_eta_pf) then
      call hchksum(eta_pf_1, "BT eta_PF_1",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
      call hchksum(d_eta_pf, "BT d_eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
    else
      call hchksum(eta_pf, "BT eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
      call hchksum(eta_pf_in, "BT eta_PF_in",g%HI,haloshift=0, scale=gv%H_to_m)
    endif
    call uvchksum("BT Cor_ref_[uv]", cor_ref_u, cor_ref_v, cs%debug_BT_HI, haloshift=0)
    call uvchksum("BT [uv]hbt0", uhbt0, vhbt0, cs%debug_BT_HI, &
                  haloshift=0, scale=gv%H_to_m)
    if (.not. use_bt_cont) then
      call uvchksum("BT Dat[uv]", datu, datv, cs%debug_BT_HI, haloshift=1, &
                    scale=gv%H_to_m)
    endif
    call uvchksum("BT wt_[uv]", wt_u, wt_v, g%HI, 0, .true., .true.)
    call uvchksum("BT frhat[uv]", cs%frhatu, cs%frhatv, g%HI, 0, .true., .true.)
    call uvchksum("BT bc_accel_[uv]", bc_accel_u, bc_accel_v, &
                  g%HI, haloshift=0)
    call uvchksum("BT IDat[uv]", cs%IDatu, cs%IDatv, g%HI, haloshift=0)
    call uvchksum("BT visc_rem_[uv]", visc_rem_u, visc_rem_v, &
                  g%HI, haloshift=1)
  endif

  if (query_averaging_enabled(cs%diag)) then
    if (cs%id_eta_st > 0) call post_data(cs%id_eta_st, eta(isd:ied,jsd:jed), cs%diag)
    if (cs%id_ubt_st > 0) call post_data(cs%id_ubt_st, ubt(isdb:iedb,jsd:jed), cs%diag)
    if (cs%id_vbt_st > 0) call post_data(cs%id_vbt_st, vbt(isd:ied,jsdb:jedb), cs%diag)
  endif

  if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
  if (id_clock_calc > 0) call cpu_clock_begin(id_clock_calc)

  if (project_velocity) then ; eta_pf_bt => eta ; else ; eta_pf_bt => eta_pred ; endif

  if (cs%dt_bt_filter >= 0.0) then
    dt_filt = 0.5 * max(0.0, min(cs%dt_bt_filter, 2.0*dt))
  else
    dt_filt = 0.5 * max(0.0, dt * min(-cs%dt_bt_filter, 2.0))
  endif
  nfilter = ceiling(dt_filt / dtbt)

  if (nstep+nfilter==0 ) call mom_error(fatal, &
      "btstep: number of barotropic step (nstep+nfilter) is 0")

  ! Set up the normalized weights for the filtered velocity.
  sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
  allocate(wt_vel(nstep+nfilter)) ; allocate(wt_eta(nstep+nfilter))
  allocate(wt_trans(nstep+nfilter+1)) ; allocate(wt_accel(nstep+nfilter+1))
  allocate(wt_accel2(nstep+nfilter+1))
  do n=1,nstep+nfilter
    ! Modify this to use a different filter...
    if ( (n==nstep) .or. (dt_filt - abs(n-nstep)*dtbt >= 0.0)) then
      wt_vel(n) = 1.0  ; wt_eta(n) = 1.0
    elseif (dtbt + dt_filt - abs(n-nstep)*dtbt > 0.0) then
      wt_vel(n) = 1.0 + (dt_filt / dtbt) - abs(n-nstep) ; wt_eta(n) = wt_vel(n)
    else
      wt_vel(n) = 0.0  ; wt_eta(n) = 0.0
    endif
!###    if (n < nstep-nfilter) then ; wt_vel(n) = 0.0 ; else ; wt_vel(n) = 1.0 ; endif
!###    if (n < nstep-nfilter) then ; wt_eta(n) = 0.0 ; else ; wt_eta(n) = 1.0 ; endif

    ! The rest should not be changed.
    sum_wt_vel = sum_wt_vel + wt_vel(n) ; sum_wt_eta = sum_wt_eta + wt_eta(n)
  enddo
  wt_trans(nstep+nfilter+1) = 0.0 ; wt_accel(nstep+nfilter+1) = 0.0
  do n=nstep+nfilter,1,-1
    wt_trans(n) = wt_trans(n+1) + wt_eta(n)
    wt_accel(n) = wt_accel(n+1) + wt_vel(n)
    sum_wt_accel = sum_wt_accel + wt_accel(n) ; sum_wt_trans = sum_wt_trans + wt_trans(n)
  enddo
  ! Normalize the weights.
  i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_accel = 1.0 / sum_wt_accel
  i_sum_wt_eta = 1.0 / sum_wt_eta ; i_sum_wt_trans = 1.0 / sum_wt_trans
  do n=1,nstep+nfilter
    wt_vel(n) = wt_vel(n) * i_sum_wt_vel
    wt_accel2(n) = wt_accel(n)
    wt_accel(n) = wt_accel(n) * i_sum_wt_accel
    wt_eta(n) = wt_eta(n) * i_sum_wt_eta
!    wt_trans(n) = wt_trans(n) * I_sum_wt_trans
  enddo


  sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0

  ! The following loop contains all of the time steps.
  isv=is ; iev=ie ; jsv=js ; jev=je
  do n=1,nstep+nfilter

    sum_wt_vel = sum_wt_vel + wt_vel(n)
    sum_wt_eta = sum_wt_eta + wt_eta(n)
    sum_wt_accel = sum_wt_accel + wt_accel2(n)
    sum_wt_trans = sum_wt_trans + wt_trans(n)

    if (cs%clip_velocity) then
      do j=jsv,jev ; do i=isv-1,iev
        if ((ubt(i,j) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then
          ! Add some error reporting later.
          ubt(i,j) = (-0.95*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
        elseif ((ubt(i,j) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
          ! Add some error reporting later.
          ubt(i,j) = (0.95*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
        endif
      enddo ; enddo
      do j=jsv-1,jev ; do i=isv,iev
        if ((vbt(i,j) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then
          ! Add some error reporting later.
          vbt(i,j) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
        elseif ((vbt(i,j) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
          ! Add some error reporting later.
          vbt(i,j) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
        endif
      enddo ; enddo
    endif

    if ((iev - stencil < ie) .or. (jev - stencil < je)) then
      if (id_clock_calc > 0) call cpu_clock_end(id_clock_calc)
      if (id_clock_pass_step > 0) call cpu_clock_begin(id_clock_pass_step)
      call do_group_pass(cs%pass_eta_ubt, cs%BT_Domain)
      isv = isvf ; iev = ievf ; jsv = jsvf ; jev = jevf
      if (id_clock_pass_step > 0) call cpu_clock_end(id_clock_pass_step)
      if (id_clock_calc > 0) call cpu_clock_begin(id_clock_calc)
    else
      isv = isv+stencil ; iev = iev-stencil
      jsv = jsv+stencil ; jev = jev-stencil
    endif

    if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
        (cs%Nonlin_cont_update_period > 0)) then
      if ((n>1) .and. (mod(n-1,cs%Nonlin_cont_update_period) == 0)) &
        call find_face_areas(datu, datv, g, gv, cs, ms, eta, 1+iev-ie)
    endif

!GOMP parallel default(none) shared(CS,isv,iev,jsv,jev,project_velocity,use_BT_cont, &
!GOMP                               uhbt,vhbt,ubt,BTCL_u,uhbt0,vbt,BTCL_v,vhbt0,     &
!GOMP                               eta_pred,eta,eta_src,dtbt,Datu,Datv,p_surf_dyn,  &
!GOMP                               dyn_coef_eta,find_etaav,is,ie,js,je,eta_sum,     &
!GOMP                               wt_accel2,n,eta_PF_BT,interp_eta_PF,wt_end,      &
!GOMP                               Instep,eta_PF,eta_PF_1,d_eta_PF,                 &
!GOMP                               apply_OBC_flather,ubt_old,vbt_old )
    if (cs%dynamic_psurf .or. .not.project_velocity) then
      if (use_bt_cont) then
!GOMP do
        do j=jsv-1,jev+1 ; do i=isv-2,iev+1
          uhbt(i,j) = find_uhbt(ubt(i,j),btcl_u(i,j)) + uhbt0(i,j)
        enddo ; enddo
!GOMP do
        do j=jsv-2,jev+1 ; do i=isv-1,iev+1
          vhbt(i,j) = find_vhbt(vbt(i,j),btcl_v(i,j)) + vhbt0(i,j)
        enddo ; enddo
!GOMP do
        do j=jsv-1,jev+1 ; do i=isv-1,iev+1
          eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
                     ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
        enddo ; enddo
      else
!GOMP do
        do j=jsv-1,jev+1 ; do i=isv-1,iev+1
          eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
              (((datu(i-1,j)*ubt(i-1,j) + uhbt0(i-1,j)) - &
                (datu(i,j)*ubt(i,j) + uhbt0(i,j))) + &
               ((datv(i,j-1)*vbt(i,j-1) + vhbt0(i,j-1)) - &
                (datv(i,j)*vbt(i,j) + vhbt0(i,j))))
        enddo ; enddo
      endif

      if (cs%dynamic_psurf) then
!GOMP do
        do j=jsv-1,jev+1 ; do i=isv-1,iev+1
          p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j))
        enddo ; enddo
      endif
    endif

    ! Recall that just outside the do n loop, there is code like...
    !  eta_PF_BT => eta_pred ; if (project_velocity) eta_PF_BT => eta

    if (find_etaav) then
!GOMP do
      do j=js,je ; do i=is,ie
        eta_sum(i,j) = eta_sum(i,j) + wt_accel2(n) * eta_pf_bt(i,j)
      enddo ; enddo
    endif

    if (interp_eta_pf) then
      wt_end = n*instep  ! This could be (n-0.5)*Instep.
!GOMP do
      do j=jsv-1,jev+1 ; do i=isv-1,iev+1
        eta_pf(i,j) = eta_pf_1(i,j) + wt_end*d_eta_pf(i,j)
      enddo ; enddo
    endif

    if (apply_obc_flather .or. apply_obc_open) then
!GOMP do
      do j=jsv,jev ; do i=isv-2,iev+1
        ubt_old(i,j) = ubt(i,j)
      enddo ; enddo
!GOMP do
      do j=jsv-2,jev+1 ; do i=isv,iev
        vbt_old(i,j) = vbt(i,j)
      enddo ; enddo
    endif
!GOMP end parallel

    if (apply_obcs) then
      if (mod(n+g%first_direction,2)==1) then
        ioff = 1; joff = 0
      else
        ioff = 0; joff = 1
      endif

      if (cs%BT_OBC%apply_u_OBCs) then  ! save the old value of ubt and uhbt
!GOMP parallel do default(none) shared(isv,iev,jsv,jev,ioff,joff,ubt_prev,ubt,uhbt_prev,  &
!GOMP                                  uhbt,ubt_sum_prev,ubt_sum,uhbt_sum_prev, &
!GOMP                                  uhbt_sum,ubt_wtd_prev,ubt_wtd)
        do j=jsv-joff,jev+joff ; do i=isv-1,iev
          ubt_prev(i,j) = ubt(i,j); uhbt_prev(i,j) = uhbt(i,j)
          ubt_sum_prev(i,j)=ubt_sum(i,j); uhbt_sum_prev(i,j)=uhbt_sum(i,j) ; ubt_wtd_prev(i,j)=ubt_wtd(i,j)
        enddo ; enddo
      endif

      if (cs%BT_OBC%apply_v_OBCs) then  ! save the old value of vbt and vhbt
!GOMP parallel do default(none) shared(isv,iev,jsv,jev,ioff,joff,vbt_prev,vbt,vhbt_prev,  &
!GOMP                                  vhbt,vbt_sum_prev,vbt_sum,vhbt_sum_prev, &
!GOMP                                  vhbt_sum,vbt_wtd_prev,vbt_wtd)
        do j=jsv-1,jev ; do i=isv-ioff,iev+ioff
          vbt_prev(i,j) = vbt(i,j); vhbt_prev(i,j) = vhbt(i,j)
          vbt_sum_prev(i,j)=vbt_sum(i,j); vhbt_sum_prev(i,j)=vhbt_sum(i,j) ; vbt_wtd_prev(i,j) = vbt_wtd(i,j)
        enddo ; enddo
      endif
    endif

!GOMP parallel default(none) shared(isv,iev,jsv,jev,G,amer,ubt,cmer,bmer,dmer,eta_PF_BT, &
!GOMP                               eta_PF,gtot_N,gtot_S,dgeo_de,CS,p_surf_dyn,n,        &
!GOMP                               v_accel_bt,wt_accel,wt_accel2,vbt,bt_rem_v,          &
!GOMP                               BT_force_v,vhbt,Cor_ref_v,dtbt,trans_wt1,trans_wt2,  &
!GOMP                               use_BT_cont,BTCL_v,vhbt0,Datv,wt_vel,azon,bzon,czon, &
!GOMP                               dzon,Cor_ref_u,gtot_E,gtot_W,u_accel_bt,bt_rem_u,    &
!GOMP                               BT_force_u,uhbt,BTCL_u,uhbt0,Datu,Cor_u,Cor_v,       &
!GOMP                               PFu,PFv,ubt_trans,vbt_trans,apply_v_OBCs,            &
!GOMP                               vbt_prev,vhbt_prev,apply_u_OBCs,ubt_prev,uhbt_prev ) &
!GOMP                       private(vel_prev)
    if (mod(n+g%first_direction,2)==1) then
      ! On odd-steps, update v first.
!GOMP do
      do j=jsv-1,jev ; do i=isv-1,iev+1
        cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
               (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
        pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
                     (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
                   dgeo_de * cs%IdyCv(i,j)
      enddo ; enddo
      if (cs%dynamic_psurf) then
!GOMP do
        do j=jsv-1,jev ; do i=isv-1,iev+1
          pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
        enddo ; enddo
      endif

      if (cs%BT_OBC%apply_v_OBCs) then  ! zero out PF across boundary
!GOMP do
        do j=jsv-1,jev ; do i=isv-1,iev+1 ; if (obc%segnum_v(i,j) /= obc_none) then
          pfv(i,j) = 0.0
        endif ; enddo ; enddo
      endif
!GOMP do
      do j=jsv-1,jev ; do i=isv-1,iev+1
        v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
        vel_prev = vbt(i,j)
        vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
                    dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
        vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
      enddo ; enddo

      if (use_bt_cont) then
!GOMP do
        do j=jsv-1,jev ; do i=isv-1,iev+1
          vhbt(i,j) = find_vhbt(vbt_trans(i,j),btcl_v(i,j)) + vhbt0(i,j)
        enddo ; enddo
      else
!GOMP do
        do j=jsv-1,jev ; do i=isv-1,iev+1
          vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
        enddo ; enddo
      endif
      if (cs%BT_OBC%apply_v_OBCs) then  ! copy back the value for v-points on the boundary.
!GOMP do
        do j=jsv-1,jev ; do i=isv-1,iev+1 ; if (obc%segnum_v(i,j) /= obc_none) then
          vbt(i,j) = vbt_prev(i,j); vhbt(i,j) = vhbt_prev(i,j)
        endif ; enddo ; enddo
      endif
      ! Now update the zonal velocity.
!GOMP do
      do j=jsv,jev ; do i=isv-1,iev
        cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
                      (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
                     cor_ref_u(i,j)
        pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
                     (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
                    dgeo_de * cs%IdxCu(i,j)
      enddo ; enddo

      if (cs%dynamic_psurf) then
!GOMP do
        do j=jsv,jev ; do i=isv-1,iev
          pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
        enddo ; enddo
      endif

      if (cs%BT_OBC%apply_u_OBCs) then  ! zero out pressure force across boundary
!GOMP do
        do j=jsv,jev ; do i=isv-1,iev ; if (obc%segnum_u(i,j) /= obc_none) then
          pfu(i,j) = 0.0
        endif ; enddo ; enddo
      endif
!GOMP do
      do j=jsv,jev ; do i=isv-1,iev
        u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
        vel_prev = ubt(i,j)
        ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
             dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
        ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
      enddo ; enddo

      if (use_bt_cont) then
!GOMP do
        do j=jsv,jev ; do i=isv-1,iev
          uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j)) + uhbt0(i,j)
        enddo ; enddo
      else
!GOMP do
        do j=jsv,jev ; do i=isv-1,iev
          uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
        enddo ; enddo
      endif
     if (cs%BT_OBC%apply_u_OBCs) then  ! copy back the value for u-points on the boundary.
!GOMP do
        do j=jsv,jev ; do i=isv-1,iev ; if (obc%segnum_u(i,j) /= obc_none) then
          ubt(i,j) = ubt_prev(i,j); uhbt(i,j) = uhbt_prev(i,j)
        endif ; enddo ; enddo
      endif
    else
      ! On even steps, update u first.
!GOMP do
      do j=jsv-1,jev+1 ; do i=isv-1,iev
        cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
                      (bzon(i,j) * vbt(i,j) +  dzon(i,j) * vbt(i+1,j-1))) - &
                     cor_ref_u(i,j)
        pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
                     (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
                     dgeo_de * cs%IdxCu(i,j)
      enddo ; enddo

      if (cs%dynamic_psurf) then
!GOMP do
        do j=jsv-1,jev+1 ; do i=isv-1,iev
          pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
        enddo ; enddo
      endif

      if (cs%BT_OBC%apply_u_OBCs) then  ! zero out pressure force across boundary
!GOMP do
        do j=jsv,jev ; do i=isv-1,iev ; if (obc%segnum_u(i,j) /= obc_none) then
          pfu(i,j) = 0.0
        endif ; enddo ; enddo
      endif

!GOMP do
      do j=jsv-1,jev+1 ; do i=isv-1,iev
        u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
        vel_prev = ubt(i,j)
        ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
             dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
        ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
      enddo ; enddo

      if (use_bt_cont) then
!GOMP do
        do j=jsv-1,jev+1 ; do i=isv-1,iev
          uhbt(i,j) = find_uhbt(ubt_trans(i,j),btcl_u(i,j)) + uhbt0(i,j)
        enddo ; enddo
      else
!GOMP do
        do j=jsv-1,jev+1 ; do i=isv-1,iev
          uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
        enddo ; enddo
      endif
      if (cs%BT_OBC%apply_u_OBCs) then  ! copy back the value for u-points on the boundary.
!GOMP do
        do j=jsv-1,jev+1 ; do i=isv-1,iev ; if (obc%segnum_u(i,j) /= obc_none) then
          ubt(i,j) = ubt_prev(i,j); uhbt(i,j) = uhbt_prev(i,j)
        endif ; enddo ; enddo
      endif

      ! Now update the meridional velocity.
!GOMP do
      do j=jsv-1,jev ; do i=isv,iev
        cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
                (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
        pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
                     (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
                    dgeo_de * cs%IdyCv(i,j)
      enddo ; enddo

      if (cs%dynamic_psurf) then
!GOMP do
        do j=jsv-1,jev ; do i=isv,iev
          pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
        enddo ; enddo
      endif

      if (cs%BT_OBC%apply_v_OBCs) then  ! zero out PF across boundary
!GOMP do
        do j=jsv-1,jev ; do i=isv-1,iev+1 ; if (obc%segnum_v(i,j) /= obc_none) then
          pfv(i,j) = 0.0
        endif ; enddo ; enddo
      endif

!GOMP do
      do j=jsv-1,jev ; do i=isv,iev
        v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
        vel_prev = vbt(i,j)
        vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
             dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
        vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
      enddo ; enddo
      if (use_bt_cont) then
!GOMP do
        do j=jsv-1,jev ; do i=isv,iev
          vhbt(i,j) = find_vhbt(vbt_trans(i,j),btcl_v(i,j)) + vhbt0(i,j)
        enddo ; enddo
      else
!GOMP do
        do j=jsv-1,jev ; do i=isv,iev
          vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
        enddo ; enddo
      endif
      if (cs%BT_OBC%apply_v_OBCs) then  ! copy back the value for v-points on the boundary.
!GOMP do
        do j=jsv-1,jev ; do i=isv,iev ; if (obc%segnum_v(i,j) /= obc_none) then
          vbt(i,j) = vbt_prev(i,j); vhbt(i,j) = vhbt_prev(i,j)
        endif ; enddo ; enddo
      endif
    endif
!GOMP end parallel

!GOMP parallel default(none) shared(is,ie,js,je,find_PF,PFu_bt_sum,wt_accel2, &
!GOMP                               PFu,PFv_bt_sum,PFv,find_Cor,Coru_bt_sum,  &
!GOMP                               Cor_u,Corv_bt_sum,Cor_v,ubt_sum,wt_trans, &
!GOMP                               ubt_trans,uhbt_sum,uhbt,ubt_wtd,wt_vel,   &
!GOMP                               ubt,vbt_sum,vbt_trans,vhbt_sum,vhbt,      &
!GOMP                               vbt_wtd,vbt,n )
    if (find_pf) then
!GOMP do
      do j=js,je ; do i=is-1,ie
        pfu_bt_sum(i,j)  = pfu_bt_sum(i,j) + wt_accel2(n) * pfu(i,j)
      enddo ; enddo
!GOMP do
      do j=js-1,je ; do i=is,ie
        pfv_bt_sum(i,j)  = pfv_bt_sum(i,j) + wt_accel2(n) * pfv(i,j)
      enddo ; enddo
    endif
    if (find_cor) then
!GOMP do
      do j=js,je ; do i=is-1,ie
        coru_bt_sum(i,j) = coru_bt_sum(i,j) + wt_accel2(n) * cor_u(i,j)
      enddo ; enddo
!GOMP do
      do j=js-1,je ; do i=is,ie
        corv_bt_sum(i,j) = corv_bt_sum(i,j) + wt_accel2(n) * cor_v(i,j)
      enddo ; enddo
    endif

!GOMP do
    do j=js,je ; do i=is-1,ie
      ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
      uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
      ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
    enddo ; enddo
!GOMP do
    do j=js-1,je ; do i=is,ie
      vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
      vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
      vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
    enddo ; enddo
!GOMP end parallel

    if (apply_obcs) then
      if (cs%BT_OBC%apply_u_OBCs) then  ! copy back the value for u-points on the boundary.
!GOMP parallel do default(none) shared(is,ie,js,je,ubt_sum_prev,ubt_sum,uhbt_sum_prev,&
!GOMP                                  uhbt_sum,ubt_wtd_prev,ubt_wtd)
        do j=js,je ; do i=is-1,ie
          if (obc%segnum_u(i,j) /= obc_none) then
            ubt_sum(i,j)=ubt_sum_prev(i,j); uhbt_sum(i,j)=uhbt_sum_prev(i,j) ; ubt_wtd(i,j)=ubt_wtd_prev(i,j)
          endif
        enddo ; enddo
      endif

      if (cs%BT_OBC%apply_v_OBCs) then  ! copy back the value for v-points on the boundary.
!GOMP parallel do default(none) shared(is,ie,js,je,vbt_sum_prev,vbt_sum,vhbt_sum_prev, &
!GOMP                                  vhbt_sum,vbt_wtd_prev,vbt_wtd)
        do j=js-1,je ; do i=is,ie
          if (obc%segnum_v(i,j) /= obc_none) then
            vbt_sum(i,j)=vbt_sum_prev(i,j); vhbt_sum(i,j)=vhbt_sum_prev(i,j) ; vbt_wtd(i,j)=vbt_wtd_prev(i,j)
          endif
        enddo ; enddo
      endif

      call apply_velocity_obcs(obc, ubt, vbt, uhbt, vhbt, &
           ubt_trans, vbt_trans, eta, ubt_old, vbt_old, cs%BT_OBC, &
           g, ms, iev-ie, dtbt, bebt, use_bt_cont, datu, datv, btcl_u, btcl_v, &
           uhbt0, vhbt0)
      if (cs%BT_OBC%apply_u_OBCs) then ; do j=js,je ; do i=is-1,ie
        if (obc%segnum_u(i,j) /= obc_none) then
          ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
          uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
          ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
        endif
      enddo ; enddo ; endif
      if (cs%BT_OBC%apply_v_OBCs) then ; do j=js-1,je ; do i=is,ie
        if (obc%segnum_v(i,j) /= obc_none) then
          vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
          vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
          vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
        endif
      enddo ; enddo ; endif
    endif

    if (cs%debug_bt) then
      call uvchksum("BT [uv]hbt just after OBC", uhbt, vhbt, &
                    cs%debug_BT_HI, haloshift=iev-ie, scale=gv%H_to_m)
    endif

!$OMP parallel do default(none) shared(isv,iev,jsv,jev,n,eta,eta_src,dtbt,CS,uhbt,vhbt,eta_wtd,wt_eta)
    do j=jsv,jev ; do i=isv,iev
      eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
                 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
      eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
      ! Should there be a concern if eta drops below 0 or G%bathyT?
    enddo ; enddo
    if (apply_obcs) call apply_eta_obcs(obc, eta, ubt_old, vbt_old, cs%BT_OBC, &
                                        g, ms, iev-ie, dtbt)

    if (do_hifreq_output) then
      time_step_end = time_bt_start + set_time(int(floor(n*dtbt+0.5)))
      call enable_averaging(dtbt, time_step_end, cs%diag)
      if (cs%id_ubt_hifreq > 0) call post_data(cs%id_ubt_hifreq, ubt(isdb:iedb,jsd:jed), cs%diag)
      if (cs%id_vbt_hifreq > 0) call post_data(cs%id_vbt_hifreq, vbt(isd:ied,jsdb:jedb), cs%diag)
      if (cs%id_eta_hifreq > 0) call post_data(cs%id_eta_hifreq, eta(isd:ied,jsd:jed), cs%diag)
      if (cs%id_uhbt_hifreq > 0) call post_data(cs%id_uhbt_hifreq, uhbt(isdb:iedb,jsd:jed), cs%diag)
      if (cs%id_vhbt_hifreq > 0) call post_data(cs%id_vhbt_hifreq, vhbt(isd:ied,jsdb:jedb), cs%diag)
      if (cs%id_eta_pred_hifreq > 0) call post_data(cs%id_eta_pred_hifreq, eta_pf_bt(isd:ied,jsd:jed), cs%diag)
    endif

    if (cs%debug_bt) then
      write(mesg,'("BT step ",I4)') n
      call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, &
                    cs%debug_BT_HI, haloshift=iev-ie)
      call hchksum(eta, trim(mesg)//" eta",cs%debug_BT_HI,haloshift=iev-ie, scale=gv%H_to_m)
    endif

  enddo ! end of do n=1,ntimestep
  if (id_clock_calc > 0) call cpu_clock_end(id_clock_calc)
  if (id_clock_calc_post > 0) call cpu_clock_begin(id_clock_calc_post)

  ! Reset the time information in the diag type.
  if (do_hifreq_output) call enable_averaging(time_int_in, time_end_in, cs%diag)

  i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_eta = 1.0 / sum_wt_eta
  i_sum_wt_accel = 1.0 / sum_wt_accel ; i_sum_wt_trans = 1.0 / sum_wt_trans

  if (find_etaav) then ; do j=js,je ; do i=is,ie
    etaav(i,j) = eta_sum(i,j) * i_sum_wt_accel
  enddo ; enddo ; endif
  do j=js-1,je+1 ; do i=is-1,ie+1 ; e_anom(i,j) = 0.0 ; enddo ; enddo
  if (interp_eta_pf) then
    do j=js,je ; do i=is,ie
      e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - &
                               (eta_pf_1(i,j) + 0.5*d_eta_pf(i,j)))
    enddo ; enddo
  else
    do j=js,je ; do i=is,ie
      e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - eta_pf(i,j))
    enddo ; enddo
  endif
  if (apply_obcs) then
    !!! Not safe for wide halos...
    if (cs%BT_OBC%apply_u_OBCs) then  ! copy back the value for u-points on the boundary.
!GOMP parallel do default(none) shared(is,ie,js,je,ubt_sum_prev,ubt_sum,uhbt_sum_prev,&
!GOMP                                  uhbt_sum,ubt_wtd_prev,ubt_wtd)
      do j=js,je ; do i=is-1,ie
        if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
          e_anom(i+1,j) = e_anom(i,j)
        elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
          e_anom(i,j) = e_anom(i+1,j)
        endif
      enddo ; enddo
    endif

    if (cs%BT_OBC%apply_v_OBCs) then  ! copy back the value for v-points on the boundary.
!GOMP parallel do default(none) shared(is,ie,js,je,vbt_sum_prev,vbt_sum,vhbt_sum_prev, &
!GOMP                                  vhbt_sum,vbt_wtd_prev,vbt_wtd)
      do j=js-1,je ; do i=is,ie
        if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
          e_anom(i,j+1) = e_anom(i,j)
        elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
          e_anom(i,j) = e_anom(i,j+1)
        endif
      enddo ; enddo
    endif
  endif

  ! It is possible that eta_out and eta_in are the same.
  do j=js,je ; do i=is,ie
    eta_out(i,j) = eta_wtd(i,j) * i_sum_wt_eta
  enddo ; enddo

  if (id_clock_calc_post > 0) call cpu_clock_end(id_clock_calc_post)
  if (id_clock_pass_post > 0) call cpu_clock_begin(id_clock_pass_post)
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_e_anom, g%Domain)
  else
    if (find_etaav) call do_group_pass(cs%pass_etaav, g%Domain)
    call do_group_pass(cs%pass_e_anom, g%Domain)
  endif
  if (id_clock_pass_post > 0) call cpu_clock_end(id_clock_pass_post)
  if (id_clock_calc_post > 0) call cpu_clock_begin(id_clock_calc_post)

  do j=js,je ; do i=is-1,ie
    cs%ubtav(i,j) = ubt_sum(i,j) * i_sum_wt_trans
    uhbtav(i,j) = uhbt_sum(i,j) * i_sum_wt_trans
 !###   u_accel_bt(I,j) = u_accel_bt(I,j) * I_sum_wt_accel
    ubt_wtd(i,j) = ubt_wtd(i,j) * i_sum_wt_vel
  enddo ; enddo

  do j=js-1,je ; do i=is,ie
    cs%vbtav(i,j) = vbt_sum(i,j) * i_sum_wt_trans
    vhbtav(i,j) = vhbt_sum(i,j) * i_sum_wt_trans
 !###   v_accel_bt(i,J) = v_accel_bt(i,J)  * I_sum_wt_accel
    vbt_wtd(i,j) = vbt_wtd(i,j) * i_sum_wt_vel
  enddo ; enddo

  if (id_clock_calc_post > 0) call cpu_clock_end(id_clock_calc_post)
  if (id_clock_pass_post > 0) call cpu_clock_begin(id_clock_pass_post)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_e_anom, g%Domain)
    if (find_etaav) call start_group_pass(cs%pass_etaav, g%Domain)
    call start_group_pass(cs%pass_ubta_uhbta, g%Domain)
  else
    call do_group_pass(cs%pass_ubta_uhbta, g%Domain)
  endif
  if (id_clock_pass_post > 0) call cpu_clock_end(id_clock_pass_post)
  if (id_clock_calc_post > 0) call cpu_clock_begin(id_clock_calc_post)

! Now calculate each layer's accelerations.
  if (apply_obcs) then
    !!! Not safe for wide halos...
    if (cs%BT_OBC%apply_u_OBCs) then ; do j=js,je ; do i=is-1,ie
      if (obc%segnum_u(i,j) /= obc_none) then
        u_accel_bt(i,j) = (ubt_wtd(i,j) - ubt_first(i,j)) / dt
      endif
    enddo ; enddo ; endif
    if (cs%BT_OBC%apply_v_OBCs) then ; do j=js-1,je ; do i=is,ie
      if (obc%segnum_v(i,j) /= obc_none) then
        v_accel_bt(i,j) = (vbt_wtd(i,j) - vbt_first(i,j)) / dt
      endif
    enddo ; enddo ; endif
  endif
!$OMP parallel do default(none) shared(is,ie,js,je,nz,accel_layer_u,u_accel_bt,pbce,gtot_W, &
!$OMP                                  e_anom,gtot_E,CS,accel_layer_v,v_accel_bt,      &
!$OMP                                  gtot_S,gtot_N)
  do k=1,nz
    do j=js,je ; do i=is-1,ie
      accel_layer_u(i,j,k) = u_accel_bt(i,j) - &
           ((pbce(i+1,j,k) - gtot_w(i+1,j)) * e_anom(i+1,j) - &
            (pbce(i,j,k) - gtot_e(i,j)) * e_anom(i,j)) * cs%IdxCu(i,j)
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      accel_layer_v(i,j,k) = v_accel_bt(i,j) - &
           ((pbce(i,j+1,k) - gtot_s(i,j+1))*e_anom(i,j+1) - &
            (pbce(i,j,k) - gtot_n(i,j))*e_anom(i,j)) * cs%IdyCv(i,j)
    enddo ; enddo
  enddo

  if (id_clock_calc_post > 0) call cpu_clock_end(id_clock_calc_post)

  ! Calculate diagnostic quantities.
  if (query_averaging_enabled(cs%diag)) then

    do j=js,je ; do i=is-1,ie ; cs%ubt_IC(i,j) = ubt_wtd(i,j) ; enddo ; enddo
    do j=js-1,je ; do i=is,ie ; cs%vbt_IC(i,j) = vbt_wtd(i,j) ; enddo ; enddo
    if (use_bt_cont) then
      do j=js,je ; do i=is-1,ie
        cs%uhbt_IC(i,j) = find_uhbt(ubt_wtd(i,j), btcl_u(i,j)) + uhbt0(i,j)
      enddo ; enddo
      do j=js-1,je ; do i=is,ie
        cs%vhbt_IC(i,j) = find_vhbt(vbt_wtd(i,j), btcl_v(i,j)) + vhbt0(i,j)
      enddo ; enddo
    else
      do j=js,je ; do i=is-1,ie
        cs%uhbt_IC(i,j) = ubt_wtd(i,j) * datu(i,j) + uhbt0(i,j)
      enddo ; enddo
      do j=js-1,je ; do i=is,ie
        cs%vhbt_IC(i,j) = vbt_wtd(i,j) * datv(i,j) + vhbt0(i,j)
      enddo ; enddo
    endif

!  Offer various barotropic terms for averaging.
    if (cs%id_PFu_bt > 0) then
      do j=js,je ; do i=is-1,ie
        pfu_bt_sum(i,j) = pfu_bt_sum(i,j) * i_sum_wt_accel
      enddo ; enddo
      call post_data(cs%id_PFu_bt, pfu_bt_sum(isdb:iedb,jsd:jed), cs%diag)
    endif
    if (cs%id_PFv_bt > 0) then
      do j=js-1,je ; do i=is,ie
        pfv_bt_sum(i,j) = pfv_bt_sum(i,j) * i_sum_wt_accel
      enddo ; enddo
      call post_data(cs%id_PFv_bt, pfv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
    endif
    if (cs%id_Coru_bt > 0) then
      do j=js,je ; do i=is-1,ie
        coru_bt_sum(i,j) = coru_bt_sum(i,j) * i_sum_wt_accel
      enddo ; enddo
      call post_data(cs%id_Coru_bt, coru_bt_sum(isdb:iedb,jsd:jed), cs%diag)
    endif
    if (cs%id_Corv_bt > 0) then
      do j=js-1,je ; do i=is,ie
        corv_bt_sum(i,j) = corv_bt_sum(i,j) * i_sum_wt_accel
      enddo ; enddo
      call post_data(cs%id_Corv_bt, corv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
    endif
    if (cs%id_ubtforce > 0) call post_data(cs%id_ubtforce, bt_force_u(isdb:iedb,jsd:jed), cs%diag)
    if (cs%id_vbtforce > 0) call post_data(cs%id_vbtforce, bt_force_v(isd:ied,jsdb:jedb), cs%diag)
    if (cs%id_uaccel > 0) call post_data(cs%id_uaccel, u_accel_bt(isdb:iedb,jsd:jed), cs%diag)
    if (cs%id_vaccel > 0) call post_data(cs%id_vaccel, v_accel_bt(isd:ied,jsdb:jedb), cs%diag)

    if (cs%id_eta_cor > 0) call post_data(cs%id_eta_cor, cs%eta_cor, cs%diag)
    if (cs%id_eta_bt > 0) call post_data(cs%id_eta_bt, eta_out, cs%diag)
    if (cs%id_gtotn > 0) call post_data(cs%id_gtotn, gtot_n(isd:ied,jsd:jed), cs%diag)
    if (cs%id_gtots > 0) call post_data(cs%id_gtots, gtot_s(isd:ied,jsd:jed), cs%diag)
    if (cs%id_gtote > 0) call post_data(cs%id_gtote, gtot_e(isd:ied,jsd:jed), cs%diag)
    if (cs%id_gtotw > 0) call post_data(cs%id_gtotw, gtot_w(isd:ied,jsd:jed), cs%diag)
    if (cs%id_ubt > 0) call post_data(cs%id_ubt, ubt_wtd(isdb:iedb,jsd:jed), cs%diag)
    if (cs%id_vbt > 0) call post_data(cs%id_vbt, vbt_wtd(isd:ied,jsdb:jedb), cs%diag)
    if (cs%id_ubtav > 0) call post_data(cs%id_ubtav, cs%ubtav, cs%diag)
    if (cs%id_vbtav > 0) call post_data(cs%id_vbtav, cs%vbtav, cs%diag)
    if (cs%id_visc_rem_u > 0) call post_data(cs%id_visc_rem_u, visc_rem_u, cs%diag)
    if (cs%id_visc_rem_v > 0) call post_data(cs%id_visc_rem_v, visc_rem_v, cs%diag)

    if (cs%id_frhatu > 0) call post_data(cs%id_frhatu, cs%frhatu, cs%diag)
    if (cs%id_uhbt > 0) call post_data(cs%id_uhbt, uhbtav, cs%diag)
    if (cs%id_frhatv > 0) call post_data(cs%id_frhatv, cs%frhatv, cs%diag)
    if (cs%id_vhbt > 0) call post_data(cs%id_vhbt, vhbtav, cs%diag)
    if (cs%id_uhbt0 > 0) call post_data(cs%id_uhbt0, uhbt0(isdb:iedb,jsd:jed), cs%diag)
    if (cs%id_vhbt0 > 0) call post_data(cs%id_vhbt0, vhbt0(isd:ied,jsdb:jedb), cs%diag)

    if (cs%id_frhatu1 > 0) call post_data(cs%id_frhatu1, cs%frhatu1, cs%diag)
    if (cs%id_frhatv1 > 0) call post_data(cs%id_frhatv1, cs%frhatv1, cs%diag)

    if (use_bt_cont) then
      if (cs%id_BTC_FA_u_EE > 0) call post_data(cs%id_BTC_FA_u_EE, bt_cont%FA_u_EE, cs%diag)
      if (cs%id_BTC_FA_u_E0 > 0) call post_data(cs%id_BTC_FA_u_E0, bt_cont%FA_u_E0, cs%diag)
      if (cs%id_BTC_FA_u_W0 > 0) call post_data(cs%id_BTC_FA_u_W0, bt_cont%FA_u_W0, cs%diag)
      if (cs%id_BTC_FA_u_WW > 0) call post_data(cs%id_BTC_FA_u_WW, bt_cont%FA_u_WW, cs%diag)
      if (cs%id_BTC_uBT_EE > 0) call post_data(cs%id_BTC_uBT_EE, bt_cont%uBT_EE, cs%diag)
      if (cs%id_BTC_uBT_WW > 0) call post_data(cs%id_BTC_uBT_WW, bt_cont%uBT_WW, cs%diag)
      if (cs%id_BTC_FA_v_NN > 0) call post_data(cs%id_BTC_FA_v_NN, bt_cont%FA_v_NN, cs%diag)
      if (cs%id_BTC_FA_v_N0 > 0) call post_data(cs%id_BTC_FA_v_N0, bt_cont%FA_v_N0, cs%diag)
      if (cs%id_BTC_FA_v_S0 > 0) call post_data(cs%id_BTC_FA_v_S0, bt_cont%FA_v_S0, cs%diag)
      if (cs%id_BTC_FA_v_SS > 0) call post_data(cs%id_BTC_FA_v_SS, bt_cont%FA_v_SS, cs%diag)
      if (cs%id_BTC_vBT_NN > 0) call post_data(cs%id_BTC_vBT_NN, bt_cont%vBT_NN, cs%diag)
      if (cs%id_BTC_vBT_SS > 0) call post_data(cs%id_BTC_vBT_SS, bt_cont%vBT_SS, cs%diag)
    endif
  else
    if (cs%id_frhatu1 > 0) cs%frhatu1(:,:,:) = cs%frhatu(:,:,:)
    if (cs%id_frhatv1 > 0) cs%frhatv1(:,:,:) = cs%frhatv(:,:,:)
  endif

  if (g%nonblocking_updates) then
    if (find_etaav) call complete_group_pass(cs%pass_etaav, g%Domain)
    call complete_group_pass(cs%pass_ubta_uhbta, g%Domain)
  endif

end subroutine btstep

!> This subroutine automatically determines an optimal value for dtbt based
!! on some state of the ocean.
subroutine set_dtbt(G, GV, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
  type(ocean_grid_type),                 intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),               intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(barotropic_cs),                   pointer       :: CS   !< Barotropic control structure.
  real, dimension(SZI_(G),SZJ_(G)),         intent(in), optional :: eta   !< The barotropic free surface height
                                                               !! anomaly or column mass anomaly, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in), optional :: pbce  !< The baroclinic pressure anomaly in each
                                                                          !! layer due to free surface height
                                                                          !! anomalies, in m2 H-1 s-2.
  type(bt_cont_type),                    pointer,    optional :: BT_cont  !< A structure with elements that describe
                                                                          !! the effective open face areas as a
                                                                          !! function of barotropic flow.
  real,                                  intent(in), optional :: gtot_est !< An estimate of the total gravitational
                                                                          !! acceleration, in m s-2.
  real,                                  intent(in), optional :: SSH_add  !< An additional contribution to SSH to
                                                                          !! provide a margin of error when
                                                                          !! calculating the external wave speed, in m.

  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: &
    gtot_E, &     ! gtot_X is the effective total reduced gravity used to relate
    gtot_W, &     ! free surface height deviations to pressure forces (including
    gtot_N, &     ! GFS and baroclinic  contributions) in the barotropic momentum
    gtot_S        ! equations half a grid-point in the X-direction (X is N, S,
                  ! E, or W) from the thickness point. gtot_X has units of m2 H-1 s-2.
                  ! (See Hallberg, J Comp Phys 1997 for a discussion.)
  real, dimension(SZIBS_(G),SZJ_(G)) :: &
    Datu          ! Basin depth at u-velocity grid points times the y-grid
                  ! spacing, in m2.
  real, dimension(SZI_(G),SZJBS_(G)) :: &
    Datv          ! Basin depth at v-velocity grid points times the x-grid
                  ! spacing, in m2.
  real :: det_de  ! The partial derivative due to self-attraction and loading
                  ! of the reference geopotential with the sea surface height.
                  ! This is typically ~0.09 or less.
  real :: dgeo_de ! The constant of proportionality between geopotential and
                  ! sea surface height.  It is a nondimensional number of
                  ! order 1.  For stability, this may be made larger
                  ! than physical problem would suggest.
  real :: add_SSH ! An additional contribution to SSH to provide a margin of error
                  ! when calculating the external wave speed, in m.
  real :: min_max_dt2, Idt_max2, dtbt_max
  logical :: use_BT_cont
  type(memory_size_type) :: MS

  character(len=200) :: mesg
  integer :: i, j, k, is, ie, js, je, nz

  if (.not.associated(cs)) call mom_error(fatal, &
      "set_dtbt: Module MOM_barotropic must be initialized before it is used.")
  if (.not.cs%split) return
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  ms%isdw = g%isd ; ms%iedw = g%ied ; ms%jsdw = g%jsd ; ms%jedw = g%jed

  if (.not.(present(pbce) .or. present(gtot_est))) call mom_error(fatal, &
      "set_dtbt: Either pbce or gtot_est must be present.")

  add_ssh = 0.0 ; if (present(ssh_add)) add_ssh = ssh_add

  use_bt_cont = .false.
  if (present(bt_cont)) use_bt_cont = (associated(bt_cont))

  if (use_bt_cont) then
    call bt_cont_to_face_areas(bt_cont, datu, datv, g, ms, 0, .true.)
  elseif (cs%Nonlinear_continuity .and. present(eta)) then
    call find_face_areas(datu, datv, g, gv, cs, ms, eta=eta, halo=0)
  else
    call find_face_areas(datu, datv, g, gv, cs, ms, halo=0, add_max=add_ssh)
  endif

  det_de = 0.0
  if (cs%tides) call tidal_forcing_sensitivity(g, cs%tides_CSp, det_de)
  dgeo_de = 1.0 + max(0.0, det_de + cs%G_extra)
  if (present(pbce)) then
    do j=js,je ; do i=is,ie
      gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
      gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
    enddo ; enddo
    do k=1,nz ; do j=js,je ; do i=is,ie
      gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * cs%frhatu(i,j,k)
      gtot_w(i,j) = gtot_w(i,j) + pbce(i,j,k) * cs%frhatu(i-1,j,k)
      gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * cs%frhatv(i,j,k)
      gtot_s(i,j) = gtot_s(i,j) + pbce(i,j,k) * cs%frhatv(i,j-1,k)
    enddo ; enddo ; enddo
  else
    do j=js,je ; do i=is,ie
      gtot_e(i,j) = gtot_est * gv%H_to_m ; gtot_w(i,j) = gtot_est * gv%H_to_m
      gtot_n(i,j) = gtot_est * gv%H_to_m ; gtot_s(i,j) = gtot_est * gv%H_to_m
    enddo ; enddo
  endif

  min_max_dt2 = 1.0e38  ! A huge number.
  do j=js,je ; do i=is,ie
    !   This is pretty accurate for gravity waves, but it is a conservative
    ! estimate since it ignores the stabilizing effect of the bottom drag.
    idt_max2 = 0.5 * (1.0 + 2.0*cs%bebt) * (g%IareaT(i,j) * &
      ((gtot_e(i,j)*datu(i,j)*g%IdxCu(i,j) + gtot_w(i,j)*datu(i-1,j)*g%IdxCu(i-1,j)) + &
       (gtot_n(i,j)*datv(i,j)*g%IdyCv(i,j) + gtot_s(i,j)*datv(i,j-1)*g%IdyCv(i,j-1))) + &
      ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
       (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)))
    if (idt_max2 * min_max_dt2 > 1.0) min_max_dt2 = 1.0 / idt_max2
  enddo ; enddo
  dtbt_max = sqrt(min_max_dt2 / dgeo_de)
  if (id_clock_sync > 0) call cpu_clock_begin(id_clock_sync)
  call min_across_pes(dtbt_max)
  if (id_clock_sync > 0) call cpu_clock_end(id_clock_sync)

  cs%dtbt = cs%dtbt_fraction * dtbt_max
  cs%dtbt_max = dtbt_max
end subroutine set_dtbt

!> The following 4 subroutines apply the open boundary conditions.
!! This subroutine applies the open boundary conditions on barotropic
!! velocities and mass transports, as developed by Mehmet Ilicak.
subroutine apply_velocity_obcs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, &
                               eta, ubt_old, vbt_old, BT_OBC, &
                               G, MS, halo, dtbt, bebt, use_BT_cont, Datu, Datv, &
                               BTCL_u, BTCL_v, uhbt0, vhbt0)
  type(ocean_obc_type),                  pointer       :: OBC     !< An associated pointer to an OBC type.
  type(ocean_grid_type),                 intent(inout) :: G       !< The ocean's grid structure.
  type(memory_size_type),                intent(in)    :: MS      !< A type that describes the memory sizes of
                                                                  !! the argument arrays.
  real, dimension(SZIBW_(MS),SZJW_(MS)), intent(inout) :: ubt     !< the zonal barotropic velocity, in m s-1.
  real, dimension(SZIBW_(MS),SZJW_(MS)), intent(inout) :: uhbt    !< the zonal barotropic transport, in H m2 s-1.
  real, dimension(SZIBW_(MS),SZJW_(MS)), intent(inout) :: ubt_trans !< the zonal barotropic velocity used in
                                                                  !! transport, m s-1.
  real, dimension(SZIW_(MS),SZJBW_(MS)), intent(inout) :: vbt     !< the meridional barotropic velocity, in m s-1.
  real, dimension(SZIW_(MS),SZJBW_(MS)), intent(inout) :: vhbt    !< the meridional barotropic transport, in H m2 s-1.
  real, dimension(SZIW_(MS),SZJBW_(MS)), intent(inout) :: vbt_trans !< the meridional BT velocity used in transports,
                                                                  !! m s-1.
  real, dimension(SZIW_(MS),SZJW_(MS)),  intent(in)    :: eta     !< The barotropic free surface height anomaly or
                                                                  !! column mass anomaly, in m or kg m-2.
  real, dimension(SZIBW_(MS),SZJW_(MS)), intent(in)    :: ubt_old !< The starting value of ubt in a barotropic step,
                                                                  !! m s-1.
  real, dimension(SZIW_(MS),SZJBW_(MS)), intent(in)    :: vbt_old !< The starting value of vbt in a barotropic step,
                                                                  !! m s-1.
  type(bt_obc_type),                     intent(in)    :: BT_OBC  !< A structure with the private barotropic arrays
                                                                  !! related to the open boundary conditions,
                                                                  !! set by set_up_BT_OBC.
  integer,                               intent(in)    :: halo    !< The extra halo size to use here.
  real,                                  intent(in)    :: dtbt    !< The time step, in s.
  real,                                  intent(in)    :: bebt    !< The fractional weighting of the future velocity
                                                                  !! in determining the transport.
  logical,                               intent(in)    :: use_BT_cont !< If true, use the BT_cont_types to calculate
                                                                  !! transports.
  real, dimension(SZIBW_(MS),SZJW_(MS)), intent(in)    :: Datu    !< A fixed estimate of the face areas at u points.
  real, dimension(SZIW_(MS),SZJBW_(MS)), intent(in)    :: Datv    !< A fixed estimate of the face areas at v points.
  type(local_bt_cont_u_type), dimension(SZIBW_(MS),SZJW_(MS)), intent(in) :: BTCL_u !< Structure of information used
                                                                  !! for a dynamic estimate of the face areas at
                                                                  !! u-points.
  type(local_bt_cont_v_type), dimension(SZIW_(MS),SZJBW_(MS)), intent(in) :: BTCL_v !< Structure of information used
                                                                  !! for a dynamic estimate of the face areas at
                                                                  !! v-points.
  real, dimension(SZIBW_(MS),SZJW_(MS)), intent(in)    :: uhbt0
  real, dimension(SZIW_(MS),SZJBW_(MS)), intent(in)    :: vhbt0

  ! Local variables
  real :: vel_prev    ! The previous velocity in m s-1.
  real :: vel_trans   ! The combination of the previous and current velocity
                      ! that does the mass transport, in m s-1.
  real :: H_u         ! The total thickness at the u-point, in m or kg m-2.
  real :: H_v         ! The total thickness at the v-point, in m or kg m-2.
  real :: cfl         ! The CFL number at the point in question, ND.
  real :: u_inlet
  real :: v_inlet
  real :: h_in
  real :: cff, Cx, Cy, tau
  real :: dhdt, dhdx, dhdy
  integer :: i, j, is, ie, js, je
  real, dimension(SZIB_(G),SZJB_(G)) :: grad
  real, parameter :: eps = 1.0e-20
  real :: rx_max, ry_max ! coefficients for radiation
  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
  rx_max = obc%rx_max ; ry_max = obc%rx_max

  if (bt_obc%apply_u_OBCs) then
    do j=js,je ; do i=is-1,ie ; if (obc%segnum_u(i,j) /= obc_none) then
      if (obc%segment(obc%segnum_u(i,j))%specified) then
        uhbt(i,j) = bt_obc%uhbt(i,j)
        ubt(i,j) = bt_obc%ubt_outer(i,j)
        vel_trans = ubt(i,j)
      elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
        if (obc%segment(obc%segnum_u(i,j))%Flather) then
          cfl = dtbt * bt_obc%Cg_u(i,j) * g%IdxCu(i,j)           ! CFL
          u_inlet = cfl*ubt_old(i-1,j) + (1.0-cfl)*ubt_old(i,j)  ! Valid for cfl<1
        !  h_in = 2.0*cfl*eta(i,j) + (1.0-2.0*cfl)*eta(i+1,j)    ! external
          h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j))      ! internal
          h_u = bt_obc%H_u(i,j)
          vel_prev = ubt(i,j)
          ubt(i,j) = 0.5*((u_inlet + bt_obc%ubt_outer(i,j)) + &
              (bt_obc%Cg_u(i,j)/h_u) * (h_in-bt_obc%eta_outer_u(i,j)))
          vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
        elseif (obc%segment(obc%segnum_u(i,j))%oblique) then
          grad(i,j) = (ubt_old(i,j+1) - ubt_old(i,j)) * g%mask2dBu(i,j)
          grad(i,j-1) = (ubt_old(i,j) - ubt_old(i,j-1)) * g%mask2dBu(i,j-1)
          grad(i-1,j) = (ubt_old(i-1,j+1) - ubt_old(i-1,j)) * g%mask2dBu(i-1,j)
          grad(i-1,j-1) = (ubt_old(i-1,j) - ubt_old(i-1,j-1)) * g%mask2dBu(i-1,j-1)
          dhdt = ubt_old(i-1,j)-ubt(i-1,j) !old-new
          dhdx = ubt(i-1,j)-ubt(i-2,j) !in new time backward sasha for I-1
!         if (OBC%segment(OBC%segnum_u(I,j))%oblique) then
            if (dhdt*(grad(i-1,j) + grad(i-1,j-1)) > 0.0) then
              dhdy = grad(i-1,j-1)
            elseif (dhdt*(grad(i-1,j) + grad(i-1,j-1)) == 0.0) then
              dhdy = 0.0
            else
              dhdy = grad(i-1,j)
            endif
!         endif
          if (dhdt*dhdx < 0.0) dhdt = 0.0
          if (dhdx == 0.0) dhdx=eps  ! avoid segv
          cx = min(dhdt/dhdx,rx_max) ! default to normal flow only
!         Cy = 0
          cff = max(dhdx*dhdx, eps)
!         if (OBC%segment(OBC%segnum_u(I,j))%oblique) then
            cff = max(dhdx*dhdx + dhdy*dhdy, eps)
            if (dhdy==0.) dhdy=eps ! avoid segv
            cy = min(cff, max(dhdt/dhdy, -cff))
!         endif
          ubt(i,j) = ((cff*ubt_old(i,j) + cx*ubt(i-1,j)) - &
              (max(cy,0.0)*grad(i,j-1) + min(cy,0.0)*grad(i,j))) / (cff + cx)
          vel_trans = ubt(i,j)
        elseif (obc%segment(obc%segnum_u(i,j))%gradient) then
          ubt(i,j) = ubt(i-1,j)
          vel_trans = ubt(i,j)
        endif
      elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
        if (obc%segment(obc%segnum_u(i,j))%Flather) then
          cfl = dtbt * bt_obc%Cg_u(i,j) * g%IdxCu(i,j)           ! CFL
          u_inlet = cfl*ubt_old(i+1,j) + (1.0-cfl)*ubt_old(i,j)  ! Valid for cfl<1
!         h_in = 2.0*cfl*eta(i+1,j) + (1.0-2.0*cfl)*eta(i,j)     ! external
          h_in = eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j))  ! external

          h_u = bt_obc%H_u(i,j)
          vel_prev = ubt(i,j)
          ubt(i,j) = 0.5*((u_inlet+bt_obc%ubt_outer(i,j)) + &
              (bt_obc%Cg_u(i,j)/h_u) * (bt_obc%eta_outer_u(i,j)-h_in))

          vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
        elseif (obc%segment(obc%segnum_u(i,j))%oblique) then
          grad(i,j) = (ubt_old(i,j+1) - ubt_old(i,j)) * g%mask2dBu(i,j)
          grad(i,j-1) = (ubt_old(i,j) - ubt_old(i,j-1)) * g%mask2dBu(i,j-1)
          grad(i+1,j) = (ubt_old(i+1,j+1) - ubt_old(i+1,j)) * g%mask2dBu(i+1,j)
          grad(i+1,j-1) = (ubt_old(i+1,j) - ubt_old(i+1,j-1)) * g%mask2dBu(i+1,j-1)
          dhdt = ubt_old(i+1,j)-ubt(i+1,j) !old-new
          dhdx = ubt(i+1,j)-ubt(i+2,j) !in new time backward sasha for I+1
!         if (OBC%segment(OBC%segnum_u(I,j))%oblique) then
            if (dhdt*(grad(i+1,j) + grad(i+1,j-1)) > 0.0) then
              dhdy = grad(i+1,j-1)
            elseif (dhdt*(grad(i+1,j) + grad(i+1,j-1)) == 0.0) then
              dhdy = 0.0
            else
              dhdy = grad(i+1,j)
            endif
!         endif
          if (dhdt*dhdx < 0.0) dhdt = 0.0
          if (dhdx == 0.0) dhdx=eps  ! avoid segv
          cx = min(dhdt/dhdx,rx_max) ! default to normal flow only
!         Cy = 0
          cff = max(dhdx*dhdx, eps)
!         if (OBC%segment(OBC%segnum_u(I,j))%oblique) then
            cff = max(dhdx*dhdx + dhdy*dhdy, eps)
            if (dhdy==0.) dhdy=eps ! avoid segv
            cy = min(cff,max(dhdt/dhdy,-cff))
!         endif
          ubt(i,j) = ((cff*ubt_old(i,j) + cx*ubt(i+1,j)) - &
              (max(cy,0.0)*grad(i,j-1) + min(cy,0.0)*grad(i,j))) / (cff + cx)
!         vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(I,j)
          vel_trans = ubt(i,j)
        elseif (obc%segment(obc%segnum_u(i,j))%gradient) then
          ubt(i,j) = ubt(i+1,j)
          vel_trans = ubt(i,j)
        endif
      endif

      if (.not. obc%segment(obc%segnum_u(i,j))%specified) then
        if (use_bt_cont) then
          uhbt(i,j) = find_uhbt(vel_trans,btcl_u(i,j)) + uhbt0(i,j)
        else
          uhbt(i,j) = datu(i,j)*vel_trans + uhbt0(i,j)
        endif
      endif

      ubt_trans(i,j) = vel_trans
    endif ; enddo ; enddo
  endif

  if (bt_obc%apply_v_OBCs) then
    do j=js-1,je ; do i=is,ie ; if (obc%segnum_v(i,j) /= obc_none) then
      if (obc%segment(obc%segnum_v(i,j))%specified) then
        vhbt(i,j) = bt_obc%vhbt(i,j)
        vbt(i,j) = bt_obc%vbt_outer(i,j)
        vel_trans = vbt(i,j)
      elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
        if (obc%segment(obc%segnum_v(i,j))%Flather) then
          cfl = dtbt * bt_obc%Cg_v(i,j) * g%IdyCv(i,j)            ! CFL
          v_inlet = cfl*vbt_old(i,j-1) + (1.0-cfl)*vbt_old(i,j)  ! Valid for cfl<1
        !  h_in = 2.0*cfl*eta(i,j) + (1.0-2.0*cfl)*eta(i,j+1)    ! external
          h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1))      ! internal

          h_v = bt_obc%H_v(i,j)
          vel_prev = vbt(i,j)
          vbt(i,j) = 0.5*((v_inlet+bt_obc%vbt_outer(i,j)) + &
              (bt_obc%Cg_v(i,j)/h_v) * (h_in-bt_obc%eta_outer_v(i,j)))

          vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
        elseif (obc%segment(obc%segnum_v(i,j))%oblique) then
          grad(i,j) = (vbt_old(i+1,j) - vbt_old(i,j)) * g%mask2dBu(i,j)
          grad(i-1,j) = (vbt_old(i,j) - vbt_old(i-1,j)) * g%mask2dBu(i-1,j)
          grad(i,j-1) = (vbt_old(i+1,j-1) - vbt_old(i,j-1)) * g%mask2dBu(i,j-1)
          grad(i-1,j-1) = (vbt_old(i,j-1) - vbt_old(i-1,j-1)) * g%mask2dBu(i-1,j-1)
          dhdt = vbt_old(i,j-1)-vbt(i,j-1) !old-new
          dhdy = vbt(i,j-1)-vbt(i,j-2) !in new time backward sasha for J-1
!         if (OBC%segment(OBC%segnum_v(i,J))%oblique) then
            if (dhdt*(grad(i,j-1) + grad(i-1,j-1)) > 0.0) then
              dhdx = grad(i-1,j-1)
            elseif (dhdt*(grad(i,j-1) + grad(i-1,j-1)) == 0.0) then
              dhdx = 0.0
            else
              dhdx = grad(i,j-1)
            endif
!         endif
          if (dhdt*dhdy < 0.0) dhdt = 0.0
          if (dhdy == 0.0) dhdy=eps  ! avoid segv
          cy = min(dhdt/dhdy,rx_max) ! default to normal flow only
!         Cx = 0
          cff = max(dhdy*dhdy, eps)
!         if (OBC%segment(OBC%segnum_v(i,J))%oblique) then
            cff = max(dhdx*dhdx + dhdy*dhdy, eps)
            if (dhdx==0.) dhdx=eps ! avoid segv
            cx = min(cff,max(dhdt/dhdx,-cff))
!         endif
          vbt(i,j) = ((cff*vbt_old(i,j) + cy*vbt(i,j-1)) - &
            (max(cx,0.0)*grad(i-1,j) + min(cx,0.0)*grad(i,j))) / (cff + cy)
!         vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,J)
          vel_trans = vbt(i,j)
        elseif (obc%segment(obc%segnum_v(i,j))%gradient) then
          vbt(i,j) = vbt(i,j-1)
          vel_trans = vbt(i,j)
        endif
      elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
        if (obc%segment(obc%segnum_v(i,j))%Flather) then
          cfl = dtbt * bt_obc%Cg_v(i,j) * g%IdyCv(i,j)            ! CFL
          v_inlet = cfl*vbt_old(i,j+1) + (1.0-cfl)*vbt_old(i,j)  ! Valid for cfl <1
        !  h_in = 2.0*cfl*eta(i,j+1) + (1.0-2.0*cfl)*eta(i,j)    ! external
          h_in = eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2))  ! internal

          h_v = bt_obc%H_v(i,j)
          vel_prev = vbt(i,j)
          vbt(i,j) = 0.5*((v_inlet+bt_obc%vbt_outer(i,j)) + &
              (bt_obc%Cg_v(i,j)/h_v) * (bt_obc%eta_outer_v(i,j)-h_in))

          vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
        elseif (obc%segment(obc%segnum_v(i,j))%oblique) then
          grad(i,j) = (vbt_old(i+1,j) - vbt_old(i,j)) * g%mask2dBu(i,j)
          grad(i-1,j) = (vbt_old(i,j) - vbt_old(i-1,j)) * g%mask2dBu(i-1,j)
          grad(i,j+1) = (vbt_old(i+1,j+1) - vbt_old(i,j+1)) * g%mask2dBu(i,j+1)
          grad(i-1,j+1) = (vbt_old(i,j+1) - vbt_old(i-1,j+1)) * g%mask2dBu(i-1,j+1)
          dhdt = vbt_old(i,j+1)-vbt(i,j+1) !old-new
          dhdy = vbt(i,j+1)-vbt(i,j+2) !in new time backward sasha for J+1
!         if (OBC%segment(OBC%segnum_v(i,J))%oblique) then
            if (dhdt*(grad(i,j+1) + grad(i-1,j+1)) > 0.0) then
              dhdx = grad(i-1,j+1)
            elseif (dhdt*(grad(i,j+1) + grad(i-1,j+1)) == 0.0) then
              dhdx = 0.0
            else
              dhdx = grad(i,j+1)
            endif
!         endif
          if (dhdt*dhdy < 0.0) dhdt = 0.0
          if (dhdy == 0.0) dhdy=eps  ! avoid segv
          cy = min(dhdt/dhdy,rx_max) ! default to normal flow only
!         Cx = 0
          cff = max(dhdy*dhdy, eps)
!         if (OBC%segment(OBC%segnum_v(i,J))%oblique) then
            cff = max(dhdx*dhdx + dhdy*dhdy, eps)
            if (dhdx==0.) dhdx=eps ! avoid segv
            cx = min(cff,max(dhdt/dhdx,-cff))
!         endif
          vbt(i,j) = ((cff*vbt_old(i,j) + cy*vbt(i,j+1)) - &
            (max(cx,0.0)*grad(i-1,j) + min(cx,0.0)*grad(i,j))) / (cff + cy)
!         vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,J)
          vel_trans = vbt(i,j)
        elseif (obc%segment(obc%segnum_v(i,j))%gradient) then
          vbt(i,j) = vbt(i,j+1)
          vel_trans = vbt(i,j)
        endif
      endif

      if (.not. obc%segment(obc%segnum_v(i,j))%specified) then
        if (use_bt_cont) then
           vhbt(i,j) = find_vhbt(vel_trans,btcl_v(i,j)) + vhbt0(i,j)
        else
           vhbt(i,j) = vel_trans*datv(i,j) + vhbt0(i,j)
        endif
      endif

      vbt_trans(i,j) = vel_trans
    endif ; enddo ; enddo
  endif

end subroutine apply_velocity_obcs

!> This subroutine applies the open boundary conditions on the free surface
!! height, as coded by Mehmet Ilicak.
subroutine apply_eta_obcs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt)
  type(ocean_obc_type),                  pointer       :: OBC  !< An associated pointer to an OBC type.
  type(memory_size_type),                intent(in)    :: MS   !< A type that describes the memory sizes of
                                                               !! the argument arrays.
  real, dimension(SZIW_(MS),SZJW_(MS)),  intent(inout) :: eta  !< The barotropic free surface height anomaly
                                                               !! or column mass anomaly, in m or kg m-2.
  real, dimension(SZIBW_(MS),SZJW_(MS)), intent(in)    :: ubt  !< the zonal barotropic velocity, in m s-1.
  real, dimension(SZIW_(MS),SZJBW_(MS)), intent(in)    :: vbt  !< the meridional barotropic velocity, in m s-1.
  type(bt_obc_type),                     intent(in)    :: BT_OBC !< A structure with the private barotropic arrays
                                                               !! related to the open boundary conditions,
                                                               !! set by set_up_BT_OBC.
  type(ocean_grid_type),                 intent(inout) :: G    !< The ocean's grid structure.
  integer,                               intent(in)    :: halo !< The extra halo size to use here.
  real,                                  intent(in)    :: dtbt !< The time step, in s.


  real :: H_u         ! The total thickness at the u-point, in m or kg m-2.
  real :: H_v         ! The total thickness at the v-point, in m or kg m-2.
  real :: cfl         ! The CFL number at the point in question, ND.
  real :: u_inlet
  real :: v_inlet
  real :: h_in
  integer :: i, j, is, ie, js, je
  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo

  if (obc%open_u_BCs_exist_globally .and. bt_obc%apply_u_OBCS) then
    do j=js,je ; do i=is-1,ie ; if (obc%segnum_u(i,j) /= obc_none) then
      if (obc%segment(obc%segnum_u(i,j))%Flather) then
        if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
          cfl = dtbt * bt_obc%Cg_u(i,j) * g%IdxCu(i,j)           ! CFL
          u_inlet = cfl*ubt(i-1,j) + (1.0-cfl)*ubt(i,j)          ! Valid for cfl <1
!          h_in = 2.0*cfl*eta(i,j) + (1.0-2.0*cfl)*eta(i+1,j)    ! external
          h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j))      ! internal

          h_u = bt_obc%H_u(i,j)
          eta(i+1,j) = 2.0 * 0.5*((bt_obc%eta_outer_u(i,j)+h_in) + &
              (h_u/bt_obc%Cg_u(i,j))*(u_inlet-bt_obc%ubt_outer(i,j))) - eta(i,j)
        elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
          cfl = dtbt*bt_obc%Cg_u(i,j)*g%IdxCu(i,j)               ! CFL
          u_inlet = cfl*ubt(i+1,j) + (1.0-cfl)*ubt(i,j)          ! Valid for cfl <1
!          h_in = 2.0*cfl*eta(i+1,j) + (1.0-2.0*cfl)*eta(i,j)    ! external
          h_in = eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j))  ! internal

          h_u = bt_obc%H_u(i,j)
          eta(i,j) = 2.0 * 0.5*((bt_obc%eta_outer_u(i,j)+h_in) + &
              (h_u/bt_obc%Cg_u(i,j))*(bt_obc%ubt_outer(i,j)-u_inlet)) - eta(i+1,j)
        endif
      elseif (obc%segment(obc%segnum_u(i,j))%radiation) then
        ! Chapman implicit from ROMS
        if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
          cfl = dtbt * bt_obc%Cg_u(i,j) * g%IdxCu(i,j)           ! CFL
          eta(i+1,j) = 1.0/(1 + cfl) * (eta(i,j) + cfl*eta(i-1,j))
        elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
          cfl = dtbt*bt_obc%Cg_u(i,j)*g%IdxCu(i,j)               ! CFL
          eta(i,j) = 1.0/(1 + cfl) * (eta(i+1,j) + cfl*eta(i+2,j))
        endif
      elseif (obc%segment(obc%segnum_u(i,j))%gradient) then
        if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
          eta(i+1,j) = eta(i,j)
        elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
          eta(i,j) = eta(i+1,j)
        endif
      endif
    endif ; enddo ; enddo
  endif

  if (obc%open_v_BCs_exist_globally .and. bt_obc%apply_v_OBCs) then
    do j=js-1,je ; do i=is,ie ; if (obc%segnum_v(i,j) /= obc_none) then
      if (obc%segment(obc%segnum_v(i,j))%Flather) then
        if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
          cfl = dtbt*bt_obc%Cg_v(i,j)*g%IdyCv(i,j)               ! CFL
          v_inlet = cfl*vbt(i,j-1) + (1.0-cfl)*vbt(i,j)          ! Valid for cfl <1
!          h_in = 2.0*cfl*eta(i,j) + (1.0-2.0*cfl)*eta(i,j+1)    ! external
          h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1))      ! internal

          h_v = bt_obc%H_v(i,j)
          eta(i,j+1) = 2.0 * 0.5*((bt_obc%eta_outer_v(i,j)+h_in) + &
              (h_v/bt_obc%Cg_v(i,j))*(v_inlet-bt_obc%vbt_outer(i,j))) - eta(i,j)
        elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
          cfl = dtbt*bt_obc%Cg_v(i,j)*g%IdyCv(i,j)               ! CFL
          v_inlet = cfl*vbt(i,j+1) + (1.0-cfl)*vbt(i,j)          ! Valid for cfl <1
!          h_in = 2.0*cfl*eta(i,j+1) + (1.0-2.0*cfl)*eta(i,j)    ! external
          h_in = eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2))  ! internal

          h_v = bt_obc%H_v(i,j)
          eta(i,j) = 2.0 * 0.5*((bt_obc%eta_outer_v(i,j)+h_in) + &
              (h_v/bt_obc%Cg_v(i,j))*(bt_obc%vbt_outer(i,j)-v_inlet)) - eta(i,j+1)
        endif
      elseif (obc%segment(obc%segnum_v(i,j))%radiation) then
        ! Chapman implicit from ROMS
        if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
          cfl = dtbt*bt_obc%Cg_v(i,j)*g%IdyCv(i,j)               ! CFL
          eta(i,j+1) = 1.0/(1 + cfl) * (eta(i,j) + cfl*eta(i,j-1))
        elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
          cfl = dtbt*bt_obc%Cg_v(i,j)*g%IdyCv(i,j)               ! CFL
          eta(i,j) = 1.0/(1 + cfl) * (eta(i,j+1) + cfl*eta(i,j+2))
        endif
      elseif (obc%segment(obc%segnum_v(i,j))%gradient) then
        if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
          eta(i,j+1) = eta(i,j)
        elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
          eta(i,j) = eta(i,j+1)
        endif
      endif
    endif ; enddo ; enddo
  endif

end subroutine apply_eta_obcs

!> This subroutine sets up the private structure used to apply the open
!! boundary conditions, as developed by Mehmet Ilicak.
subroutine set_up_bt_obc(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, Datv, BTCL_u, BTCL_v)
  type(ocean_obc_type),                  pointer       :: OBC    !< An associated pointer to an OBC type.
  type(memory_size_type),                intent(in)    :: MS     !< A type that describes the memory sizes of the
                                                                 !! argument arrays.
  real, dimension(SZIW_(MS),SZJW_(MS)),  intent(in)    :: eta    !< The barotropic free surface height anomaly or
                                                                 !! column mass anomaly, in m or kg m-2.
  type(bt_obc_type),                     intent(inout) :: BT_OBC !< A structure with the private barotropic arrays
                                                                 !! related to the open boundary conditions,
                                                                 !! set by set_up_BT_OBC.
  type(ocean_grid_type),                 intent(inout) :: G      !< The ocean's grid structure.
  type(verticalgrid_type),               intent(in)    :: GV     !< The ocean's vertical grid structure.
  integer,                               intent(in)    :: halo   !< The extra halo size to use here.
  logical,                               intent(in)    :: use_BT_cont !< If true, use the BT_cont_types to calculate
                                                                 !! transports.
  real, dimension(SZIBW_(MS),SZJW_(MS)), intent(in)    :: Datu   !< A fixed estimate of the face areas at u points.
  real, dimension(SZIW_(MS),SZJBW_(MS)), intent(in)    :: Datv   !< A fixed estimate of the face areas at u points.
  type(local_bt_cont_u_type), dimension(SZIBW_(MS),SZJW_(MS)), intent(in) :: BTCL_u !< Structure of information used
                                                                 !! for a dynamic estimate of the face areas at
                                                                 !! u-points.
  type(local_bt_cont_v_type), dimension(SZIW_(MS),SZJBW_(MS)), intent(in) :: BTCL_v !< Structure of information used
                                                                 !! for a dynamic estimate of the face areas at
                                                                 !! v-points.

  ! Local variables
  integer :: i, j, k, is, ie, js, je, n, nz, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  integer :: isdw, iedw, jsdw, jedw
  logical :: OBC_used
  type(obc_segment_type), pointer  :: segment !< Open boundary segment
  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  isdw = ms%isdw ; iedw = ms%iedw ; jsdw = ms%jsdw ; jedw = ms%jedw

  if ((isdw < isd) .or. (jsdw < jsd)) then
    call mom_error(fatal, "set_up_BT_OBC: Open boundary conditions are not "//&
                           "yet fully implemented with wide barotropic halos.")
  endif

  if (.not. bt_obc%is_alloced) then
    allocate(bt_obc%Cg_u(isdw-1:iedw,jsdw:jedw))        ; bt_obc%Cg_u(:,:) = 0.0
    allocate(bt_obc%H_u(isdw-1:iedw,jsdw:jedw))         ; bt_obc%H_u(:,:) = 0.0
    allocate(bt_obc%uhbt(isdw-1:iedw,jsdw:jedw))        ; bt_obc%uhbt(:,:) = 0.0
    allocate(bt_obc%ubt_outer(isdw-1:iedw,jsdw:jedw))   ; bt_obc%ubt_outer(:,:) = 0.0
    allocate(bt_obc%eta_outer_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%eta_outer_u(:,:) = 0.0

    allocate(bt_obc%Cg_v(isdw:iedw,jsdw-1:jedw))        ; bt_obc%Cg_v(:,:) = 0.0
    allocate(bt_obc%H_v(isdw:iedw,jsdw-1:jedw))         ; bt_obc%H_v(:,:) = 0.0
    allocate(bt_obc%vhbt(isdw:iedw,jsdw-1:jedw))        ; bt_obc%vhbt(:,:) = 0.0
    allocate(bt_obc%vbt_outer(isdw:iedw,jsdw-1:jedw))   ; bt_obc%vbt_outer(:,:) = 0.0
    allocate(bt_obc%eta_outer_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%eta_outer_v(:,:)=0.0
    bt_obc%is_alloced = .true.
  endif

  if (bt_obc%apply_u_OBCs) then
    if (obc%specified_u_BCs_exist_globally) then
      do n = 1, obc%number_of_segments
        segment => obc%segment(n)
        if (segment%is_E_or_W .and. segment%specified) then
          do j=segment%HI%jsd,segment%HI%jed ; do i=segment%HI%IsdB,segment%HI%IedB
            bt_obc%uhbt(i,j) = 0.
          enddo ; enddo
          do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed ; do i=segment%HI%IsdB,segment%HI%IedB
            bt_obc%uhbt(i,j) = bt_obc%uhbt(i,j) + segment%normal_trans(i,j,k)
          enddo ; enddo ; enddo
        endif
      enddo
    endif
    do j=js,je ; do i=is-1,ie ; if (obc%segnum_u(i,j) /= obc_none) then
      ! Can this go in segment loop above? Is loop above wrong for wide halos??
      if (obc%segment(obc%segnum_u(i,j))%specified) then
        if (use_bt_cont) then
          bt_obc%ubt_outer(i,j) = uhbt_to_ubt(bt_obc%uhbt(i,j),btcl_u(i,j))
        else
          if (datu(i,j) > 0.0) bt_obc%ubt_outer(i,j) = bt_obc%uhbt(i,j) / datu(i,j)
        endif
      else            ! This is assuming Flather as only other option
        bt_obc%Cg_u(i,j) = sqrt(gv%g_prime(1)*(0.5* &
                                (g%bathyT(i,j) + g%bathyT(i+1,j))))
        if (gv%Boussinesq) then
          if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
            bt_obc%H_u(i,j) = g%bathyT(i,j)*gv%m_to_H + eta(i,j)
          elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
            bt_obc%H_u(i,j) = g%bathyT(i+1,j)*gv%m_to_H + eta(i+1,j)
          endif
        else
          if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
            bt_obc%H_u(i,j) = eta(i,j)
          elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
            bt_obc%H_u(i,j) = eta(i+1,j)
          endif
        endif
      endif
    endif ; enddo ; enddo
    if (obc%Flather_u_BCs_exist_globally) then
     do n = 1, obc%number_of_segments
        segment => obc%segment(n)
        if (segment%is_E_or_W .and. segment%Flather) then
          do j=segment%HI%jsd,segment%HI%jed ; do i=segment%HI%IsdB,segment%HI%IedB
            bt_obc%ubt_outer(i,j) = segment%normal_vel_bt(i,j)
            bt_obc%eta_outer_u(i,j) = segment%eta(i,j)
          enddo ; enddo
        endif
      enddo
    endif
  endif

  if (bt_obc%apply_v_OBCs) then
    if (obc%specified_v_BCs_exist_globally) then
      do n = 1, obc%number_of_segments
        segment => obc%segment(n)
        if (segment%is_N_or_S .and. segment%specified) then
          do j=segment%HI%JsdB,segment%HI%JedB ; do i=segment%HI%isd,segment%HI%ied
            bt_obc%vhbt(i,j) = 0.
          enddo ; enddo
          do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB ; do i=segment%HI%isd,segment%HI%ied
            bt_obc%vhbt(i,j) = bt_obc%vhbt(i,j) + segment%normal_trans(i,j,k)
          enddo ; enddo ; enddo
        endif
      enddo
    endif
    do j=js-1,je ; do i=is,ie ; if (obc%segnum_v(i,j) /= obc_none) then
      ! Can this go in segment loop above? Is loop above wrong for wide halos??
      if (obc%segment(obc%segnum_v(i,j))%specified) then
        if (use_bt_cont) then
          bt_obc%vbt_outer(i,j) = vhbt_to_vbt(bt_obc%vhbt(i,j),btcl_v(i,j))
        else
          if (datv(i,j) > 0.0) bt_obc%vbt_outer(i,j) = bt_obc%vhbt(i,j) / datv(i,j)
        endif
      else              ! This is assuming Flather as only other option
        bt_obc%Cg_v(i,j) = sqrt(gv%g_prime(1)*(0.5* &
                                (g%bathyT(i,j) + g%bathyT(i,j+1))))
        if (gv%Boussinesq) then
          if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
            bt_obc%H_v(i,j) = g%bathyT(i,j)*gv%m_to_H + eta(i,j)
          elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
            bt_obc%H_v(i,j) = g%bathyT(i,j+1)*gv%m_to_H + eta(i,j+1)
          endif
        else
          if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
            bt_obc%H_v(i,j) = eta(i,j)
          elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
            bt_obc%H_v(i,j) = eta(i,j+1)
          endif
        endif
      endif
    endif ; enddo ; enddo
    if (obc%Flather_v_BCs_exist_globally) then
     do n = 1, obc%number_of_segments
        segment => obc%segment(n)
        if (segment%is_N_or_S .and. segment%Flather) then
          do j=segment%HI%JsdB,segment%HI%JedB ; do i=segment%HI%isd,segment%HI%ied
            bt_obc%vbt_outer(i,j) = segment%normal_vel_bt(i,j)
            bt_obc%eta_outer_v(i,j) = segment%eta(i,j)
          enddo ; enddo
        endif
      enddo
    endif
  endif

end subroutine set_up_bt_obc

!> Clean up the BT_OBC memory.
subroutine destroy_bt_obc(BT_OBC)
  type(bt_obc_type), intent(inout) :: BT_OBC !< A structure with the private barotropic arrays
                                             !! related to the open boundary conditions,
                                             !! set by set_up_BT_OBC.

  if (bt_obc%is_alloced) then
    deallocate(bt_obc%Cg_u)
    deallocate(bt_obc%H_u)
    deallocate(bt_obc%uhbt)
    deallocate(bt_obc%ubt_outer)
    deallocate(bt_obc%eta_outer_u)

    deallocate(bt_obc%Cg_v)
    deallocate(bt_obc%H_v)
    deallocate(bt_obc%vhbt)
    deallocate(bt_obc%vbt_outer)
    deallocate(bt_obc%eta_outer_v)
    bt_obc%is_alloced = .false.
  endif
end subroutine destroy_bt_obc

!> btcalc calculates the barotropic velocities from the full velocity and
!! thickness fields, determines the fraction of the total water column in each
!! layer at velocity points, and determines a corrective fictitious mass source
!! that will drive the barotropic estimate of the free surface height toward the
!! baroclinic estimate.
subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
  type(ocean_grid_type),                  intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),                intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)  :: h    !< Layer thicknesses, in H (usually m or kg m-2).
  type(barotropic_cs),                    pointer       :: CS   !< The control structure returned by a previous
                                                                !! call to barotropic_init.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: h_u !< The specified thicknesses at u-points,
                                                                !! in m or kg m-2.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in), optional :: h_v !< The specified thicknesses at v-points,
                                                                !! in m or kg m-2.
  logical,                                intent(in), optional :: may_use_default !< An optional logical argument
                                                                !! to indicate that the default velocity point
                                                                !! thickesses may be used for this particular
                                                                !! calculation, even though the setting of
                                                                !! CS%hvel_scheme would usually require that h_u
                                                                !! and h_v be passed in.
  type(ocean_obc_type),                   pointer, optional :: OBC !< Open boundary control structure.

  ! Local variables
! All of these variables are in the same units as h - usually m or kg m-2.
  real :: hatutot(szib_(g))    ! The sum of the layer thicknesses
  real :: hatvtot(szi_(g))     ! interpolated to the u & v grid points.
  real :: Ihatutot(szib_(g))   ! Ihatutot and Ihatvtot are the inverses
  real :: Ihatvtot(szi_(g))    ! of hatutot and hatvtot, both in H-1.
  real :: h_arith              ! The arithmetic mean thickness, in H.
  real :: h_harm               ! The harmonic mean thicknesses, in H.
  real :: h_neglect            ! A thickness that is so small it is usually lost
                               ! in roundoff and can be neglected, in H.
  real :: wt_arith             ! The nondimensional weight for the arithmetic
                               ! mean thickness.  The harmonic mean uses
                               ! a weight of (1 - wt_arith).
  real :: Rh                   ! A ratio of summed thicknesses, nondim.
  real :: e_u(szib_(g),szk_(g)+1) !   The interface heights at u-velocity and
  real :: e_v(szi_(g),szk_(g)+1)  ! v-velocity points in H.
  real :: D_shallow_u(szi_(g)) ! The shallower of the adjacent depths in H.
  real :: D_shallow_v(szib_(g))! The shallower of the adjacent depths in H.
  real :: htot                 ! The sum of the layer thicknesses, in H.
  real :: Ihtot                ! The inverse of htot, in H-1.

  logical :: use_default, test_dflt, apply_OBCs
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, i, j, k
  integer :: iss, ies, n

!    This section interpolates thicknesses onto u & v grid points with the
! second order accurate estimate h = 2*(h+ * h-)/(h+ + h-).
  if (.not.associated(cs)) call mom_error(fatal, &
      "btcalc: Module MOM_barotropic must be initialized before it is used.")
  if (.not.cs%split) return

  use_default = .false.
  test_dflt = .false. ; if (present(may_use_default)) test_dflt = may_use_default

  if (test_dflt) then
    if (.not.((present(h_u) .and. present(h_v)) .or. &
              (cs%hvel_scheme == harmonic) .or. (cs%hvel_scheme == hybrid) .or.&
              (cs%hvel_scheme == arithmetic))) use_default = .true.
  else
    if (.not.((present(h_u) .and. present(h_v)) .or. &
              (cs%hvel_scheme == harmonic) .or. (cs%hvel_scheme == hybrid) .or.&
              (cs%hvel_scheme == arithmetic))) call mom_error(fatal, &
        "btcalc: Inconsistent settings of optional arguments and hvel_scheme.")
  endif

  apply_obcs = .false.
  if (present(obc)) then ; if (associated(obc)) then ; if (obc%OBC_pe) then
    ! Some open boundary condition points might be in this processor's symmetric
    ! computational domain.
    apply_obcs = (obc%number_of_segments > 0)
  endif ; endif ; endif

  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
  h_neglect = gv%H_subroundoff

  !   This estimates the fractional thickness of each layer at the velocity
  ! points, using a harmonic mean estimate.
!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_u,CS,h_neglect,h,use_default,G,GV) &
!$OMP                          private(hatutot,Ihatutot,e_u,D_shallow_u,h_arith,h_harm,wt_arith)

  do j=js,je
    if (present(h_u)) then
      do i=is-1,ie ; hatutot(i) = h_u(i,j,1) ; enddo
      do k=2,nz ; do i=is-1,ie
        hatutot(i) = hatutot(i) + h_u(i,j,k)
      enddo ; enddo
      do i=is-1,ie ; ihatutot(i) = g%mask2dCu(i,j) / (hatutot(i) + h_neglect) ; enddo
      do k=1,nz ; do i=is-1,ie
        cs%frhatu(i,j,k) = h_u(i,j,k) * ihatutot(i)
      enddo ; enddo
    else
      if (cs%hvel_scheme == arithmetic) then
        do i=is-1,ie
          cs%frhatu(i,j,1) = 0.5 * (h(i+1,j,1) + h(i,j,1))
          hatutot(i) = cs%frhatu(i,j,1)
        enddo
        do k=2,nz ; do i=is-1,ie
          cs%frhatu(i,j,k) = 0.5 * (h(i+1,j,k) + h(i,j,k))
          hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
        enddo ; enddo
      elseif (cs%hvel_scheme == hybrid .or. use_default) then
        do i=is-1,ie
          e_u(i,nz+1) = -0.5 * gv%m_to_H * (g%bathyT(i+1,j) + g%bathyT(i,j))
          d_shallow_u(i) = -gv%m_to_H * min(g%bathyT(i+1,j), g%bathyT(i,j))
          hatutot(i) = 0.0
        enddo
        do k=nz,1,-1 ; do i=is-1,ie
          e_u(i,k) = e_u(i,k+1) + 0.5 * (h(i+1,j,k) + h(i,j,k))
          h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
          if (e_u(i,k+1) >= d_shallow_u(i)) then
            cs%frhatu(i,j,k) = h_arith
          else
            h_harm = (h(i+1,j,k) * h(i,j,k)) / (h_arith + h_neglect)
            if (e_u(i,k) <= d_shallow_u(i)) then
              cs%frhatu(i,j,k) = h_harm
            else
              wt_arith = (e_u(i,k) - d_shallow_u(i)) / (h_arith + h_neglect)
              cs%frhatu(i,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
            endif
          endif
          hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
        enddo ; enddo
      elseif (cs%hvel_scheme == harmonic) then
        do i=is-1,ie
          cs%frhatu(i,j,1) = 2.0*(h(i+1,j,1) * h(i,j,1)) / &
                             ((h(i+1,j,1) + h(i,j,1)) + h_neglect)
          hatutot(i) = cs%frhatu(i,j,1)
        enddo
        do k=2,nz ; do i=is-1,ie
          cs%frhatu(i,j,k) = 2.0*(h(i+1,j,k) * h(i,j,k)) / &
                             ((h(i+1,j,k) + h(i,j,k)) + h_neglect)
          hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
        enddo ; enddo
      endif
      do i=is-1,ie ; ihatutot(i) = g%mask2dCu(i,j) / (hatutot(i) + h_neglect) ; enddo
      do k=1,nz ; do i=is-1,ie
        cs%frhatu(i,j,k) = cs%frhatu(i,j,k) * ihatutot(i)
      enddo ; enddo
    endif
  enddo

!$OMP parallel do default(none) shared(is,ie,js,je,nz,CS,G,GV,h_v,h_neglect,h,use_default) &
!$OMP                          private(hatvtot,Ihatvtot,e_v,D_shallow_v,h_arith,h_harm,wt_arith)
  do j=js-1,je
    if (present(h_v)) then
      do i=is,ie ; hatvtot(i) = h_v(i,j,1) ; enddo
      do k=2,nz ; do i=is,ie
        hatvtot(i) = hatvtot(i) + h_v(i,j,k)
      enddo ; enddo
      do i=is,ie ; ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect) ; enddo
      do k=1,nz ; do i=is,ie
        cs%frhatv(i,j,k) = h_v(i,j,k) * ihatvtot(i)
      enddo ; enddo
    else
      if (cs%hvel_scheme == arithmetic) then
        do i=is,ie
          cs%frhatv(i,j,1) = 0.5 * (h(i,j+1,1) + h(i,j,1))
          hatvtot(i) = cs%frhatv(i,j,1)
        enddo
        do k=2,nz ; do i=is,ie
          cs%frhatv(i,j,k) = 0.5 * (h(i,j+1,k) + h(i,j,k))
          hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
        enddo ; enddo
      elseif (cs%hvel_scheme == hybrid .or. use_default) then
        do i=is,ie
          e_v(i,nz+1) = -0.5 * gv%m_to_H * (g%bathyT(i,j+1) + g%bathyT(i,j))
          d_shallow_v(i) = -gv%m_to_H * min(g%bathyT(i,j+1), g%bathyT(i,j))
          hatvtot(i) = 0.0
        enddo
        do k=nz,1,-1 ; do i=is,ie
          e_v(i,k) = e_v(i,k+1) + 0.5 * (h(i,j+1,k) + h(i,j,k))
          h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
          if (e_v(i,k+1) >= d_shallow_v(i)) then
            cs%frhatv(i,j,k) = h_arith
          else
            h_harm = (h(i,j+1,k) * h(i,j,k)) / (h_arith + h_neglect)
            if (e_v(i,k) <= d_shallow_v(i)) then
              cs%frhatv(i,j,k) = h_harm
            else
              wt_arith = (e_v(i,k) - d_shallow_v(i)) / (h_arith + h_neglect)
              cs%frhatv(i,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
            endif
          endif
          hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
        enddo ; enddo
      elseif (cs%hvel_scheme == harmonic) then
        do i=is,ie
          cs%frhatv(i,j,1) = 2.0*(h(i,j+1,1) * h(i,j,1)) / &
                             ((h(i,j+1,1) + h(i,j,1)) + h_neglect)
          hatvtot(i) = cs%frhatv(i,j,1)
        enddo
        do k=2,nz ; do i=is,ie
          cs%frhatv(i,j,k) = 2.0*(h(i,j+1,k) * h(i,j,k)) / &
                             ((h(i,j+1,k) + h(i,j,k)) + h_neglect)
          hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
        enddo ; enddo
      endif
      do i=is,ie ; ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect) ; enddo
      do k=1,nz ; do i=is,ie
        cs%frhatv(i,j,k) = cs%frhatv(i,j,k) * ihatvtot(i)
      enddo ; enddo
    endif
  enddo

  if (apply_obcs) then ; do n=1,obc%number_of_segments ! Test for segment type?
    if (.not. obc%segment(n)%on_pe) cycle
    if (obc%segment(n)%direction == obc_direction_n) then
      j = obc%segment(n)%HI%JsdB
      if ((j >= js-1) .and. (j <= je)) then
        iss = max(is,obc%segment(n)%HI%isd) ; ies = min(ie,obc%segment(n)%HI%ied)
        do i=iss,ies ; hatvtot(i) = h(i,j,1) ; enddo
        do k=2,nz ; do i=iss,ies
          hatvtot(i) = hatvtot(i) + h(i,j,k)
        enddo ; enddo
        do i=iss,ies
          ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect)
        enddo
        do k=1,nz ; do i=iss,ies
          cs%frhatv(i,j,k) = h(i,j,k) * ihatvtot(i)
        enddo ; enddo
      endif
    elseif (obc%segment(n)%direction == obc_direction_s) then
      j = obc%segment(n)%HI%JsdB
      if ((j >= js-1) .and. (j <= je)) then
        iss = max(is,obc%segment(n)%HI%isd) ; ies = min(ie,obc%segment(n)%HI%ied)
        do i=iss,ies ; hatvtot(i) = h(i,j+1,1) ; enddo
        do k=2,nz ; do i=iss,ies
          hatvtot(i) = hatvtot(i) + h(i,j+1,k)
        enddo ; enddo
        do i=iss,ies
          ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect)
        enddo
        do k=1,nz ; do i=iss,ies
          cs%frhatv(i,j,k) = h(i,j+1,k) * ihatvtot(i)
        enddo ; enddo
      endif
    elseif (obc%segment(n)%direction == obc_direction_e) then
      i = obc%segment(n)%HI%IsdB
      if ((i >= is-1) .and. (i <= ie)) then
        do j = max(js,obc%segment(n)%HI%jsd), min(je,obc%segment(n)%HI%jed)
          htot = h(i,j,1)
          do k=2,nz ; htot = htot + h(i,j,k) ; enddo
          ihtot = g%mask2dCu(i,j) / (htot + h_neglect)
          do k=1,nz ; cs%frhatu(i,j,k) = h(i,j,k) * ihtot ; enddo
        enddo
      endif
    elseif (obc%segment(n)%direction == obc_direction_w) then
      i = obc%segment(n)%HI%IsdB
      if ((i >= is-1) .and. (i <= ie)) then
        do j = max(js,obc%segment(n)%HI%jsd), min(je,obc%segment(n)%HI%jed)
          htot = h(i+1,j,1)
          do k=2,nz ; htot = htot + h(i+1,j,k) ; enddo
          ihtot = g%mask2dCu(i,j) / (htot + h_neglect)
          do k=1,nz ; cs%frhatu(i,j,k) = h(i+1,j,k) * ihtot ; enddo
        enddo
      endif
    else
      call mom_error(fatal, "btcalc encountered and OBC segment of indeterminate direction.")
    endif
  enddo ; endif

  if (cs%debug) then
    call uvchksum("btcalc frhat[uv]", cs%frhatu, cs%frhatv, g%HI, 0, .true., .true.)
    if (present(h_u) .and. present(h_v)) &
      call uvchksum("btcalc h_[uv]", h_u, h_v, g%HI, 0, .true., .true., scale=gv%H_to_m)
    call hchksum(h, "btcalc h",g%HI, haloshift=1, scale=gv%H_to_m)
  endif

end subroutine btcalc

!> The function find_uhbt determines the zonal transport for a given velocity.
function find_uhbt(u, BTC) result(uhbt)
  real, intent(in) :: u    !< The local zonal velocity, in m s-1
  type(local_bt_cont_u_type), intent(in) :: BTC
  real :: uhbt !< The result

  if (u == 0.0) then
    uhbt = 0.0
  elseif (u < btc%uBT_EE) then
    uhbt = (u - btc%uBT_EE) * btc%FA_u_EE + btc%uh_EE
  elseif (u < 0.0) then
    uhbt = u * (btc%FA_u_E0 + btc%uh_crvE * u**2)
  elseif (u <= btc%uBT_WW) then
    uhbt = u * (btc%FA_u_W0 + btc%uh_crvW * u**2)
  else ! (u > BTC%uBT_WW)
    uhbt = (u - btc%uBT_WW) * btc%FA_u_WW + btc%uh_WW
  endif
end function find_uhbt

!> This function inverts the transport function to determine the barotopic
!! velocity that is consistent with a given transport.
function uhbt_to_ubt(uhbt, BTC, guess) result(ubt)
  real, intent(in) :: uhbt                      !< The barotropic zonal transport that should be inverted for,
                                                !! in units of H m2 s-1.
  type(local_bt_cont_u_type), intent(in) :: BTC !< A structure containing various fields that allow the
                                                !! barotropic transports to be calculated consistently with the
                                                !! layers' continuity equations.
  real, optional, intent(in) :: guess           !< A guess at what ubt will be.  The result is not allowed
                                                !! to be dramatically larger than guess.
  real :: ubt                                   !< The result - The velocity that gives uhbt transport, in m s-1.

  ! Local variables
  real :: ubt_min, ubt_max, uhbt_err, derr_du
  real :: uherr_min, uherr_max
  real, parameter :: tol = 1.0e-10
  real :: dvel, vsr  ! Temporary variables used in the limiting the velocity.
  real, parameter :: vs1 = 1.25  ! Nondimensional parameters used in limiting
  real, parameter :: vs2 = 2.0   ! the velocity, starting at vs1, with the
                                 ! maximum increase of vs2, both nondim.
  integer :: itt, max_itt = 20

  ! Find the value of ubt that gives uhbt.
  if (uhbt == 0.0) then
    ubt = 0.0
  elseif (uhbt < btc%uh_EE) then
    ubt = btc%uBT_EE + (uhbt - btc%uh_EE) / btc%FA_u_EE
  elseif (uhbt < 0.0) then
    ! Iterate to convergence with Newton's method (when bounded) and the
    ! false position method otherwise.  ubt will be negative.
    ubt_min = btc%uBT_EE ; uherr_min = btc%uh_EE - uhbt
    ubt_max = 0.0 ; uherr_max = -uhbt
    ! Use a false-position method first guess.
    ubt = btc%uBT_EE * (uhbt / btc%uh_EE)
    do itt = 1, max_itt
      uhbt_err = ubt * (btc%FA_u_E0 + btc%uh_crvE * ubt**2) - uhbt

      if (abs(uhbt_err) < tol*abs(uhbt)) exit
      if (uhbt_err > 0.0) then ; ubt_max = ubt ; uherr_max = uhbt_err ; endif
      if (uhbt_err < 0.0) then ; ubt_min = ubt ; uherr_min = uhbt_err ; endif

      derr_du = btc%FA_u_E0 + 3.0 * btc%uh_crvE * ubt**2
      if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
          (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0)) then
        ! Use a false-position method guess.
        ubt = ubt_max + (ubt_min-ubt_max) * (uherr_max / (uherr_max-uherr_min))
      else ! Use Newton's method.
        ubt = ubt - uhbt_err / derr_du
        if (abs(uhbt_err) < (0.01*tol)*abs(ubt_min*derr_du)) exit
      endif
    enddo
  elseif (uhbt <= btc%uh_WW) then
    ! Iterate to convergence with Newton's method.  ubt will be positive.
    ubt_min = 0.0 ; uherr_min = -uhbt
    ubt_max = btc%uBT_WW ; uherr_max = btc%uh_WW - uhbt
    ! Use a false-position method first guess.
    ubt = btc%uBT_WW * (uhbt / btc%uh_WW)
    do itt = 1, max_itt
      uhbt_err = ubt * (btc%FA_u_W0 + btc%uh_crvW * ubt**2) - uhbt

      if (abs(uhbt_err) < tol*abs(uhbt)) exit
      if (uhbt_err > 0.0) then ; ubt_max = ubt ; uherr_max = uhbt_err ; endif
      if (uhbt_err < 0.0) then ; ubt_min = ubt ; uherr_min = uhbt_err ; endif

      derr_du = btc%FA_u_W0 + 3.0 * btc%uh_crvW * ubt**2
      if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
          (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0)) then
        ! Use a false-position method guess.
        ubt = ubt_min + (ubt_max-ubt_min) * (-uherr_min / (uherr_max-uherr_min))
      else ! Use Newton's method.
        ubt = ubt - uhbt_err / derr_du
        if (abs(uhbt_err) < (0.01*tol)*(ubt_max*derr_du)) exit
      endif
    enddo
  else ! (uhbt > BTC%uh_WW)
    ubt = btc%uBT_WW + (uhbt - btc%uh_WW) / btc%FA_u_WW
  endif

  if (present(guess)) then
    dvel = abs(ubt) - vs1*abs(guess)
    if (dvel > 0.0) then ! Limit the velocity
      if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) ) then
        vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1)))
      else  ! The exp be less than 4e-18 anyway in this case, so neglect it.
        vsr = vs2
      endif
      ubt = sign(vsr * guess, ubt)
    endif
  endif

end function uhbt_to_ubt

!> The function find_vhbt determines the meridional transport for a given velocity.
function find_vhbt(v, BTC) result(vhbt)
  real, intent(in) :: v    !< The local meridional velocity, in m s-1
  type(local_bt_cont_v_type), intent(in) :: BTC
  real :: vhbt !< The result

  if (v == 0.0) then
    vhbt = 0.0
  elseif (v < btc%vBT_NN) then
    vhbt = (v - btc%vBT_NN) * btc%FA_v_NN + btc%vh_NN
  elseif (v < 0.0) then
    vhbt = v * (btc%FA_v_N0 + btc%vh_crvN * v**2)
  elseif (v <= btc%vBT_SS) then
    vhbt = v * (btc%FA_v_S0 + btc%vh_crvS * v**2)
  else ! (v > BTC%vBT_SS)
    vhbt = (v - btc%vBT_SS) * btc%FA_v_SS + btc%vh_SS
  endif
end function find_vhbt

!> This function inverts the transport function to determine the barotopic
!! velocity that is consistent with a given transport.
function vhbt_to_vbt(vhbt, BTC, guess) result(vbt)
  real, intent(in) :: vhbt                      !< The barotropic meridional transport that should be
                                                !! inverted for, in units of H m2 s-1.
  type(local_bt_cont_v_type), intent(in) :: BTC !< A structure containing various fields that allow the
                                                !! barotropic transports to be calculated consistently
                                                !! with the layers' continuity equations.
  real, optional, intent(in) :: guess           !< A guess at what vbt will be. The result is not allowed
                                                !! to be dramatically larger than guess.
  real :: vbt !< The result - The velocity that gives vhbt transport, in m s-1.

  ! Local variables
  real :: vbt_min, vbt_max, vhbt_err, derr_dv
  real :: vherr_min, vherr_max
  real, parameter :: tol = 1.0e-10
  real :: dvel, vsr  ! Temporary variables used in the limiting the velocity.
  real, parameter :: vs1 = 1.25  ! Nondimensional parameters used in limiting
  real, parameter :: vs2 = 2.0   ! the velocity, starting at vs1, with the
                                 ! maximum increase of vs2, both nondim.
  integer :: itt, max_itt = 20

  ! Find the value of vbt that gives vhbt.
  if (vhbt == 0.0) then
    vbt = 0.0
  elseif (vhbt < btc%vh_NN) then
    vbt = btc%vBT_NN + (vhbt - btc%vh_NN) / btc%FA_v_NN
  elseif (vhbt < 0.0) then
    ! Iterate to convergence with Newton's method (when bounded) and the
    ! false position method otherwise.  vbt will be negative.
    vbt_min = btc%vBT_NN ; vherr_min = btc%vh_NN - vhbt
    vbt_max = 0.0 ; vherr_max = -vhbt
    ! Use a false-position method first guess.
    vbt = btc%vBT_NN * (vhbt / btc%vh_NN)
    do itt = 1, max_itt
      vhbt_err = vbt * (btc%FA_v_N0 + btc%vh_crvN * vbt**2) - vhbt

      if (abs(vhbt_err) < tol*abs(vhbt)) exit
      if (vhbt_err > 0.0) then ; vbt_max = vbt ; vherr_max = vhbt_err ; endif
      if (vhbt_err < 0.0) then ; vbt_min = vbt ; vherr_min = vhbt_err ; endif

      derr_dv = btc%FA_v_N0 + 3.0 * btc%vh_crvN * vbt**2
      if ((vhbt_err >= derr_dv*(vbt - vbt_min)) .or. &
          (-vhbt_err >= derr_dv*(vbt_max - vbt)) .or. (derr_dv <= 0.0)) then
        ! Use a false-position method guess.
        vbt = vbt_max + (vbt_min-vbt_max) * (vherr_max / (vherr_max-vherr_min))
      else ! Use Newton's method.
        vbt = vbt - vhbt_err / derr_dv
        if (abs(vhbt_err) < (0.01*tol)*abs(derr_dv*vbt_min)) exit
      endif
    enddo
  elseif (vhbt <= btc%vh_SS) then
    ! Iterate to convergence with Newton's method.  vbt will be positive.
    vbt_min = 0.0 ; vherr_min = -vhbt
    vbt_max = btc%vBT_SS ; vherr_max = btc%vh_SS - vhbt
    ! Use a false-position method first guess.
    vbt = btc%vBT_SS * (vhbt / btc%vh_SS)
    do itt = 1, max_itt
      vhbt_err = vbt * (btc%FA_v_S0 + btc%vh_crvS * vbt**2) - vhbt

      if (abs(vhbt_err) < tol*abs(vhbt)) exit
      if (vhbt_err > 0.0) then ; vbt_max = vbt ; vherr_max = vhbt_err ; endif
      if (vhbt_err < 0.0) then ; vbt_min = vbt ; vherr_min = vhbt_err ; endif

      derr_dv = btc%FA_v_S0 + 3.0 * btc%vh_crvS * vbt**2
      if ((vhbt_err >= derr_dv*(vbt - vbt_min)) .or. &
          (-vhbt_err >= derr_dv*(vbt_max - vbt)) .or. (derr_dv <= 0.0)) then
        ! Use a false-position method guess.
        vbt = vbt_min + (vbt_max-vbt_min) * (-vherr_min / (vherr_max-vherr_min))
      else ! Use Newton's method.
        vbt = vbt - vhbt_err / derr_dv
        if (abs(vhbt_err) < (0.01*tol)*(vbt_max*derr_dv)) exit
      endif
    enddo
  else ! (vhbt > BTC%vh_SS)
    vbt = btc%vBT_SS + (vhbt - btc%vh_SS) / btc%FA_v_SS
  endif

  if (present(guess)) then
    dvel = abs(vbt) - vs1*abs(guess)
    if (dvel > 0.0) then ! Limit the velocity
      if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) ) then
        vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1)))
      else  ! The exp be less than 4e-18 anyway in this case, so neglect it.
        vsr = vs2
      endif
      vbt = sign(guess * vsr, vbt)
    endif
  endif

end function vhbt_to_vbt

!> This subroutine sets up reordered versions of the BT_cont type in the
!! local_BT_cont types, which have wide halos properly filled in.
subroutine set_local_bt_cont_types(BT_cont, BTCL_u, BTCL_v, G, MS, BT_Domain, halo)
  type(bt_cont_type),                                    intent(inout) :: BT_cont    !< The BT_cont_type input to the
                                                                                     !! barotropic solver.
  type(memory_size_type),                                intent(in)    :: MS         !< A type that describes the
                                                                                     !! memory sizes of the argument
                                                                                     !! arrays.
  type(local_bt_cont_u_type), dimension(SZIBW_(MS),SZJW_(MS)), intent(out) :: BTCL_u !< A structure with the u
                                                                                     !! information from BT_cont.
  type(local_bt_cont_v_type), dimension(SZIW_(MS),SZJBW_(MS)), intent(out) :: BTCL_v !< A structure with the v
                                                                                     !! information from BT_cont.
  type(ocean_grid_type),                                 intent(in)    :: G          !< The ocean's grid structure.
  type(mom_domain_type),                                 intent(inout) :: BT_Domain  !< The domain to use for updating
                                                                                     !! the halos of wide arrays.
  integer,                                     optional, intent(in)    :: halo       !< The extra halo size to use here.

  ! Local variables
  real, dimension(SZIBW_(MS),SZJW_(MS)) :: &
    u_polarity, uBT_EE, uBT_WW, FA_u_EE, FA_u_E0, FA_u_W0, FA_u_WW
  real, dimension(SZIW_(MS),SZJBW_(MS)) :: &
    v_polarity, vBT_NN, vBT_SS, FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS
  real, parameter :: C1_3 = 1.0/3.0
  integer :: i, j, is, ie, js, je, hs
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  hs = 1 ; if (present(halo)) hs = max(halo,0)

  ! Copy the BT_cont arrays into symmetric, potentially wide haloed arrays.
!$OMP parallel default(none) shared(is,ie,js,je,hs,u_polarity,uBT_EE,uBT_WW,FA_u_EE, &
!$OMP                               FA_u_E0,FA_u_W0,FA_u_WW,v_polarity,vBT_NN,vBT_SS,&
!$OMP                               FA_v_NN,FA_v_N0,FA_v_S0,FA_v_SS,BT_cont )
!$OMP do
  do j=js-hs,je+hs ; do i=is-hs-1,ie+hs
    u_polarity(i,j) = 1.0
    ubt_ee(i,j) = 0.0 ; ubt_ww(i,j) = 0.0
    fa_u_ee(i,j) = 0.0 ; fa_u_e0(i,j) = 0.0 ; fa_u_w0(i,j) = 0.0 ; fa_u_ww(i,j) = 0.0
  enddo ; enddo
!$OMP do
  do j=js-hs-1,je+hs ; do i=is-hs,ie+hs
    v_polarity(i,j) = 1.0
    vbt_nn(i,j) = 0.0 ; vbt_ss(i,j) = 0.0
    fa_v_nn(i,j) = 0.0 ; fa_v_n0(i,j) = 0.0 ; fa_v_s0(i,j) = 0.0 ; fa_v_ss(i,j) = 0.0
  enddo ; enddo
!$OMP do
  do j=js,je; do i=is-1,ie
    ubt_ee(i,j) = bt_cont%uBT_EE(i,j) ; ubt_ww(i,j) = bt_cont%uBT_WW(i,j)
    fa_u_ee(i,j) = bt_cont%FA_u_EE(i,j) ; fa_u_e0(i,j) = bt_cont%FA_u_E0(i,j)
    fa_u_w0(i,j) = bt_cont%FA_u_W0(i,j) ; fa_u_ww(i,j) = bt_cont%FA_u_WW(i,j)
  enddo ; enddo
!$OMP do
  do j=js-1,je; do i=is,ie
    vbt_nn(i,j) = bt_cont%vBT_NN(i,j) ; vbt_ss(i,j) = bt_cont%vBT_SS(i,j)
    fa_v_nn(i,j) = bt_cont%FA_v_NN(i,j) ; fa_v_n0(i,j) = bt_cont%FA_v_N0(i,j)
    fa_v_s0(i,j) = bt_cont%FA_v_S0(i,j) ; fa_v_ss(i,j) = bt_cont%FA_v_SS(i,j)
  enddo ; enddo
!$OMP end parallel

  if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
  if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
!--- begin setup for group halo update
  call create_group_pass(bt_cont%pass_polarity_BT, u_polarity, v_polarity, bt_domain)
  call create_group_pass(bt_cont%pass_polarity_BT, ubt_ee, vbt_nn, bt_domain)
  call create_group_pass(bt_cont%pass_polarity_BT, ubt_ww, vbt_ss, bt_domain)

  call create_group_pass(bt_cont%pass_FA_uv, fa_u_ee, fa_v_nn, bt_domain, to_all+scalar_pair)
  call create_group_pass(bt_cont%pass_FA_uv, fa_u_e0, fa_v_n0, bt_domain, to_all+scalar_pair)
  call create_group_pass(bt_cont%pass_FA_uv, fa_u_w0, fa_v_s0, bt_domain, to_all+scalar_pair)
  call create_group_pass(bt_cont%pass_FA_uv, fa_u_ww, fa_v_ss, bt_domain, to_all+scalar_pair)
!--- end setup for group halo update
  ! Do halo updates on BT_cont.
  call do_group_pass(bt_cont%pass_polarity_BT, bt_domain)
  call do_group_pass(bt_cont%pass_FA_uv, bt_domain)
  if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
  if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)

!$OMP parallel default(none) shared(is,ie,js,je,hs,BTCL_u,FA_u_EE,FA_u_E0,FA_u_W0, &
!$OMP                               FA_u_WW,uBT_EE,uBT_WW,u_polarity,BTCL_v,       &
!$OMP                               FA_v_NN,FA_v_N0,FA_v_S0,FA_v_SS,vBT_NN,vBT_SS, &
!$OMP                               v_polarity )
!$OMP do
  do j=js-hs,je+hs ; do i=is-hs-1,ie+hs
    btcl_u(i,j)%FA_u_EE = fa_u_ee(i,j) ; btcl_u(i,j)%FA_u_E0 = fa_u_e0(i,j)
    btcl_u(i,j)%FA_u_W0 = fa_u_w0(i,j) ; btcl_u(i,j)%FA_u_WW = fa_u_ww(i,j)
    btcl_u(i,j)%uBT_EE = ubt_ee(i,j)   ; btcl_u(i,j)%uBT_WW = ubt_ww(i,j)
    ! Check for reversed polarity in the tripolar halo regions.
    if (u_polarity(i,j) < 0.0) then
      call swap(btcl_u(i,j)%FA_u_EE, btcl_u(i,j)%FA_u_WW)
      call swap(btcl_u(i,j)%FA_u_E0, btcl_u(i,j)%FA_u_W0)
      call swap(btcl_u(i,j)%uBT_EE,  btcl_u(i,j)%uBT_WW)
    endif

    btcl_u(i,j)%uh_EE = btcl_u(i,j)%uBT_EE * &
        (c1_3 * (2.0*btcl_u(i,j)%FA_u_E0 + btcl_u(i,j)%FA_u_EE))
    btcl_u(i,j)%uh_WW = btcl_u(i,j)%uBT_WW * &
        (c1_3 * (2.0*btcl_u(i,j)%FA_u_W0 + btcl_u(i,j)%FA_u_WW))

    btcl_u(i,j)%uh_crvE = 0.0 ; btcl_u(i,j)%uh_crvW = 0.0
    if (abs(btcl_u(i,j)%uBT_WW) > 0.0) btcl_u(i,j)%uh_crvW = &
      (c1_3 * (btcl_u(i,j)%FA_u_WW - btcl_u(i,j)%FA_u_W0)) / btcl_u(i,j)%uBT_WW**2
    if (abs(btcl_u(i,j)%uBT_EE) > 0.0) btcl_u(i,j)%uh_crvE = &
      (c1_3 * (btcl_u(i,j)%FA_u_EE - btcl_u(i,j)%FA_u_E0)) / btcl_u(i,j)%uBT_EE**2
  enddo ; enddo
!$OMP do
  do j=js-hs-1,je+hs ; do i=is-hs,ie+hs
    btcl_v(i,j)%FA_v_NN = fa_v_nn(i,j) ; btcl_v(i,j)%FA_v_N0 = fa_v_n0(i,j)
    btcl_v(i,j)%FA_v_S0 = fa_v_s0(i,j) ; btcl_v(i,j)%FA_v_SS = fa_v_ss(i,j)
    btcl_v(i,j)%vBT_NN = vbt_nn(i,j)   ; btcl_v(i,j)%vBT_SS = vbt_ss(i,j)
    ! Check for reversed polarity in the tripolar halo regions.
    if (v_polarity(i,j) < 0.0) then
      call swap(btcl_v(i,j)%FA_v_NN, btcl_v(i,j)%FA_v_SS)
      call swap(btcl_v(i,j)%FA_v_N0, btcl_v(i,j)%FA_v_S0)
      call swap(btcl_v(i,j)%vBT_NN,  btcl_v(i,j)%vBT_SS)
    endif

    btcl_v(i,j)%vh_NN = btcl_v(i,j)%vBT_NN * &
        (c1_3 * (2.0*btcl_v(i,j)%FA_v_N0 + btcl_v(i,j)%FA_v_NN))
    btcl_v(i,j)%vh_SS = btcl_v(i,j)%vBT_SS * &
        (c1_3 * (2.0*btcl_v(i,j)%FA_v_S0 + btcl_v(i,j)%FA_v_SS))

    btcl_v(i,j)%vh_crvN = 0.0 ; btcl_v(i,j)%vh_crvS = 0.0
    if (abs(btcl_v(i,j)%vBT_SS) > 0.0) btcl_v(i,j)%vh_crvS = &
      (c1_3 * (btcl_v(i,j)%FA_v_SS - btcl_v(i,j)%FA_v_S0)) / btcl_v(i,j)%vBT_SS**2
    if (abs(btcl_v(i,j)%vBT_NN) > 0.0) btcl_v(i,j)%vh_crvN = &
      (c1_3 * (btcl_v(i,j)%FA_v_NN - btcl_v(i,j)%FA_v_N0)) / btcl_v(i,j)%vBT_NN**2
  enddo ; enddo
!$OMP end parallel
end subroutine set_local_bt_cont_types


!> Adjust_local_BT_cont_types sets up reordered versions of the BT_cont type
!! in the local_BT_cont types, which have wide halos properly filled in.
subroutine adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, &
                                      G, MS, halo)
  type(memory_size_type), intent(in)  :: MS   !< A type that describes the memory sizes of the argument arrays.
  real, dimension(SZIBW_(MS),SZJW_(MS)), &
                          intent(in)  :: ubt  !< The linearization zonal barotropic velocity in m s-1.
  real, dimension(SZIBW_(MS),SZJW_(MS)), &
                          intent(in)  :: uhbt !< The linearization zonal barotropic transport in H m2 s-1.
  real, dimension(SZIW_(MS),SZJBW_(MS)), &
                          intent(in)  :: vbt  !< The linearization meridional barotropic velocity in m s-1.
  real, dimension(SZIW_(MS),SZJBW_(MS)), &
                          intent(in)  :: vhbt !< The linearization meridional barotropic transport in H m2 s-1.
  type(local_bt_cont_u_type), dimension(SZIBW_(MS),SZJW_(MS)), &
                          intent(out) :: BTCL_u !< A structure with the u information from BT_cont.
  type(local_bt_cont_v_type), dimension(SZIW_(MS),SZJBW_(MS)), &
                          intent(out) :: BTCL_v !< A structure with the v information from BT_cont.
  type(ocean_grid_type),  intent(in)  :: G    !< The ocean's grid structure.
  integer,      optional, intent(in)  :: halo !< The extra halo size to use here.

  ! Local variables
  real, dimension(SZIBW_(MS),SZJW_(MS)) :: &
    u_polarity, uBT_EE, uBT_WW, FA_u_EE, FA_u_E0, FA_u_W0, FA_u_WW
  real, dimension(SZIW_(MS),SZJBW_(MS)) :: &
    v_polarity, vBT_NN, vBT_SS, FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS
  real, parameter :: C1_3 = 1.0/3.0
  integer :: i, j, is, ie, js, je, hs
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  hs = 1 ; if (present(halo)) hs = max(halo,0)

  !$OMP parallel do default(shared)
  do j=js-hs,je+hs ; do i=is-hs-1,ie+hs
    if ((ubt(i,j) > btcl_u(i,j)%uBT_WW) .and. (uhbt(i,j) > btcl_u(i,j)%uh_WW)) then
      ! Expand the cubic fit to use this new point.  ubt is negative.
      btcl_u(i,j)%ubt_WW = ubt(i,j)
      if (3.0*uhbt(i,j) < 2.0*ubt(i,j) * btcl_u(i,j)%FA_u_W0) then
        ! No further bounding is needed.
        btcl_u(i,j)%uh_crvW = (uhbt(i,j) - ubt(i,j) * btcl_u(i,j)%FA_u_W0) / ubt(i,j)**3
      else ! This should not happen often!
        btcl_u(i,j)%FA_u_W0 = 1.5*uhbt(i,j) / ubt(i,j)
        btcl_u(i,j)%uh_crvW = -0.5*uhbt(i,j) / ubt(i,j)**3
      endif
      btcl_u(i,j)%uh_WW = uhbt(i,j)
      ! I don't know whether this is helpful.
!     BTCL_u(I,j)%FA_u_WW = min(BTCL_u(I,j)%FA_u_WW, uhbt(I,j) / ubt(I,j))
    elseif ((ubt(i,j) < btcl_u(i,j)%uBT_EE) .and. (uhbt(i,j) < btcl_u(i,j)%uh_EE)) then
      ! Expand the cubic fit to use this new point.  ubt is negative.
      btcl_u(i,j)%ubt_EE = ubt(i,j)
      if (3.0*uhbt(i,j) < 2.0*ubt(i,j) * btcl_u(i,j)%FA_u_E0) then
        ! No further bounding is needed.
        btcl_u(i,j)%uh_crvE = (uhbt(i,j) - ubt(i,j) * btcl_u(i,j)%FA_u_E0) / ubt(i,j)**3
      else ! This should not happen often!
        btcl_u(i,j)%FA_u_E0 = 1.5*uhbt(i,j) / ubt(i,j)
        btcl_u(i,j)%uh_crvE = -0.5*uhbt(i,j) / ubt(i,j)**3
      endif
      btcl_u(i,j)%uh_EE = uhbt(i,j)
      ! I don't know whether this is helpful.
!     BTCL_u(I,j)%FA_u_EE = min(BTCL_u(I,j)%FA_u_EE, uhbt(I,j) / ubt(I,j))
    endif
  enddo ; enddo
  !$OMP parallel do default(shared)
  do j=js-hs-1,je+hs ; do i=is-hs,ie+hs
    if ((vbt(i,j) > btcl_v(i,j)%vBT_SS) .and. (vhbt(i,j) > btcl_v(i,j)%vh_SS)) then
      ! Nxpand the cubic fit to use this new point.  vbt is negative.
      btcl_v(i,j)%vbt_SS = vbt(i,j)
      if (3.0*vhbt(i,j) < 2.0*vbt(i,j) * btcl_v(i,j)%FA_v_S0) then
        ! No fvrther bovnding is needed.
        btcl_v(i,j)%vh_crvS = (vhbt(i,j) - vbt(i,j) * btcl_v(i,j)%FA_v_S0) / vbt(i,j)**3
      else ! This shovld not happen often!
        btcl_v(i,j)%FA_v_S0 = 1.5*vhbt(i,j) / vbt(i,j)
        btcl_v(i,j)%vh_crvS = -0.5*vhbt(i,j) / vbt(i,j)**3
      endif
      btcl_v(i,j)%vh_SS = vhbt(i,j)
      ! I don't know whether this is helpful.
!     BTCL_v(i,J)%FA_v_SS = min(BTCL_v(i,J)%FA_v_SS, vhbt(i,J) / vbt(i,J))
    elseif ((vbt(i,j) < btcl_v(i,j)%vBT_NN) .and. (vhbt(i,j) < btcl_v(i,j)%vh_NN)) then
      ! Nxpand the cubic fit to use this new point.  vbt is negative.
      btcl_v(i,j)%vbt_NN = vbt(i,j)
      if (3.0*vhbt(i,j) < 2.0*vbt(i,j) * btcl_v(i,j)%FA_v_N0) then
        ! No fvrther bovnding is needed.
        btcl_v(i,j)%vh_crvN = (vhbt(i,j) - vbt(i,j) * btcl_v(i,j)%FA_v_N0) / vbt(i,j)**3
      else ! This shovld not happen often!
        btcl_v(i,j)%FA_v_N0 = 1.5*vhbt(i,j) / vbt(i,j)
        btcl_v(i,j)%vh_crvN = -0.5*vhbt(i,j) / vbt(i,j)**3
      endif
      btcl_v(i,j)%vh_NN = vhbt(i,j)
      ! I don't know whether this is helpful.
!     BTCL_v(i,J)%FA_v_NN = min(BTCL_v(i,J)%FA_v_NN, vhbt(i,J) / vbt(i,J))
    endif
  enddo ; enddo

end subroutine adjust_local_bt_cont_types

!> This subroutine uses the BTCL types to find typical or maximum face
!! areas, which can then be used for finding wave speeds, etc.
subroutine bt_cont_to_face_areas(BT_cont, Datu, Datv, G, MS, halo, maximize)
  type(bt_cont_type),                         intent(inout) :: BT_cont    !< The BT_cont_type input to the
                                                                          !! barotropic solver.
  type(memory_size_type),                     intent(in)    :: MS         !< A type that describes the memory
                                                                          !! sizes of the argument arrays.
  real, dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), intent(out)   :: Datu
  real, dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), intent(out)   :: Datv
  type(ocean_grid_type),                      intent(in)  :: G            !< The ocean's grid structure.
  integer,                          optional, intent(in)  :: halo         !< The extra halo size to use here.
  logical,                          optional, intent(in)  :: maximize

  ! Local variables
  logical :: find_max
  integer :: i, j, is, ie, js, je, hs
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  hs = 1 ; if (present(halo)) hs = max(halo,0)
  find_max = .false. ; if (present(maximize)) find_max = maximize

  if (find_max) then
    do j=js-hs,je+hs ; do i=is-1-hs,ie+hs
      datu(i,j) = max(bt_cont%FA_u_EE(i,j), bt_cont%FA_u_E0(i,j), &
                      bt_cont%FA_u_W0(i,j), bt_cont%FA_u_WW(i,j))
    enddo ; enddo
    do j=js-1-hs,je+hs ; do i=is-hs,ie+hs
      datv(i,j) = max(bt_cont%FA_v_NN(i,j), bt_cont%FA_v_N0(i,j), &
                      bt_cont%FA_v_S0(i,j), bt_cont%FA_v_SS(i,j))
    enddo ; enddo
  else
    do j=js-hs,je+hs ; do i=is-1-hs,ie+hs
      datu(i,j) = 0.5 * (bt_cont%FA_u_E0(i,j) + bt_cont%FA_u_W0(i,j))
    enddo ; enddo
    do j=js-1-hs,je+hs ; do i=is-hs,ie+hs
      datv(i,j) = 0.5 * (bt_cont%FA_v_N0(i,j) + bt_cont%FA_v_S0(i,j))
    enddo ; enddo
  endif

end subroutine bt_cont_to_face_areas

subroutine swap(a,b)
  real, intent(inout) :: a, b
  real :: tmp
  tmp = a ; a = b ; b = tmp
end subroutine swap

!> This subroutine determines the open face areas of cells for calculating
!! the barotropic transport.
subroutine find_face_areas(Datu, Datv, G, GV, CS, MS, eta, halo, add_max)
  type(memory_size_type),                   intent(in) :: MS
!  (in)      MS - A type that describes the memory sizes of the argument arrays.
  real, dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), intent(out)   :: Datu !< The open zonal face area,
                                                                            !! in H m (m2 or kg m-1).
  real, dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), intent(out)   :: Datv !< The open meridional face area,
                                                                            !! in H m (m2 or kg m-1).
  type(ocean_grid_type),                    intent(in) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),                  intent(in) :: GV   !< The ocean's vertical grid structure.
  type(barotropic_cs),                      pointer    :: CS   !< The control structure returned by a previous
                                                                            !! call to barotropic_init.
  real, dimension(MS%isdw:MS%iedw,MS%jsdw:MS%jedw), optional, intent(in) :: eta !< The barotropic free surface
                                                                            !! height anomaly or column mass
                                                                            !! anomaly, in H (m or kg m-2).
  integer,                        optional, intent(in) :: halo              !< The halo size to use, default = 1.
  real,                           optional, intent(in) :: add_max           !< A value to add to the maximum
                                                                            !! depth (used to overestimate the
                                                                            !! external wave speed) in m.


  ! Local variables
  real :: H1, H2      ! Temporary total thicknesses, in m or kg m-2.
  integer :: i, j, is, ie, js, je, hs
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  hs = 1 ; if (present(halo)) hs = max(halo,0)

!$OMP parallel default(none) shared(is,ie,js,je,hs,eta,GV,CS,Datu,Datv,add_max) &
!$OMP                       private(H1,H2)
  if (present(eta)) then
    ! The use of harmonic mean thicknesses ensure positive definiteness.
    if (gv%Boussinesq) then
!$OMP do
      do j=js-hs,je+hs ; do i=is-1-hs,ie+hs
        h1 = cs%bathyT(i,j)*gv%m_to_H + eta(i,j) ; h2 = cs%bathyT(i+1,j)*gv%m_to_H + eta(i+1,j)
        datu(i,j) = 0.0 ; if ((h1 > 0.0) .and. (h2 > 0.0)) &
        datu(i,j) = cs%dy_Cu(i,j) * (2.0 * h1 * h2) / (h1 + h2)
!       Datu(I,j) = CS%dy_Cu(I,j) * 0.5 * (H1 + H2)
      enddo ; enddo
!$OMP do
      do j=js-1-hs,je+hs ; do i=is-hs,ie+hs
        h1 = cs%bathyT(i,j)*gv%m_to_H + eta(i,j) ; h2 = cs%bathyT(i,j+1)*gv%m_to_H + eta(i,j+1)
        datv(i,j) = 0.0 ; if ((h1 > 0.0) .and. (h2 > 0.0)) &
        datv(i,j) = cs%dx_Cv(i,j) * (2.0 * h1 * h2) / (h1 + h2)
!       Datv(i,J) = CS%dy_v(i,J) * 0.5 * (H1 + H2)
      enddo ; enddo
    else
!$OMP do
      do j=js-hs,je+hs ; do i=is-1-hs,ie+hs
        datu(i,j) = 0.0 ; if ((eta(i,j) > 0.0) .and. (eta(i+1,j) > 0.0)) &
        datu(i,j) = cs%dy_Cu(i,j) * (2.0 * eta(i,j) * eta(i+1,j)) / &
                                  (eta(i,j) + eta(i+1,j))
        ! Datu(I,j) = CS%dy_Cu(I,j) * 0.5 * (eta(i,j) + eta(i+1,j))
      enddo ; enddo
!$OMP do
      do j=js-1-hs,je+hs ; do i=is-hs,ie+hs
        datv(i,j) = 0.0 ; if ((eta(i,j) > 0.0) .and. (eta(i,j+1) > 0.0)) &
        datv(i,j) = cs%dx_Cv(i,j) * (2.0 * eta(i,j) * eta(i,j+1)) / &
                                  (eta(i,j) + eta(i,j+1))
        ! Datv(i,J) = CS%dy_v(i,J) * 0.5 * (eta(i,j) + eta(i,j+1))
      enddo ; enddo
    endif
  elseif (present(add_max)) then
!$OMP do
    do j=js-hs,je+hs ; do i=is-1-hs,ie+hs
      datu(i,j) = cs%dy_Cu(i,j) * gv%m_to_H * &
                  (max(cs%bathyT(i+1,j), cs%bathyT(i,j)) + add_max)
    enddo ; enddo
!$OMP do
    do j=js-1-hs,je+hs ; do i=is-hs,ie+hs
      datv(i,j) = cs%dx_Cv(i,j) * gv%m_to_H * &
                  (max(cs%bathyT(i,j+1), cs%bathyT(i,j)) + add_max)
    enddo ; enddo
  else
!$OMP do
    do j=js-hs,je+hs ; do i=is-1-hs,ie+hs
      datu(i, j) = 0.0
      !Would be "if (G%mask2dCu(I,j)>0.) &" is G was valid on BT domain
      if (cs%bathyT(i+1,j)+cs%bathyT(i,j)>0.) &
        datu(i,j) = 2.0*cs%dy_Cu(i,j) * gv%m_to_H * &
                  (cs%bathyT(i+1,j) * cs%bathyT(i,j)) / &
                  (cs%bathyT(i+1,j) + cs%bathyT(i,j))
    enddo ; enddo
!$OMP do
    do j=js-1-hs,je+hs ; do i=is-hs,ie+hs
      datv(i, j) = 0.0
      !Would be "if (G%mask2dCv(i,J)>0.) &" is G was valid on BT domain
      if (cs%bathyT(i,j+1)+cs%bathyT(i,j)>0.) &
        datv(i,j) = 2.0*cs%dx_Cv(i,j) * gv%m_to_H * &
                  (cs%bathyT(i,j+1) * cs%bathyT(i,j)) / &
                  (cs%bathyT(i,j+1) + cs%bathyT(i,j))
    enddo ; enddo
  endif
!$OMP end parallel

end subroutine find_face_areas

!> bt_mass_source determines the appropriately limited mass source for
!! the barotropic solver, along with a corrective fictitious mass source that
!! will drive the barotropic estimate of the free surface height toward the
!! baroclinic estimate.
subroutine bt_mass_source(h, eta, fluxes, set_cor, dt_therm, dt_since_therm, &
                          G, GV, CS)
  type(ocean_grid_type),              intent(in) :: G        !< The ocean's grid structure.
  type(verticalgrid_type),            intent(in) :: GV       !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h  !< Layer thicknesses, in H (usually m or kg m-2).
  real, dimension(SZI_(G),SZJ_(G)),   intent(in) :: eta      !< The free surface height that is to be corrected, in m.
  type(forcing),                      intent(in) :: fluxes   !< A structure containing pointers to any possible
                                                             !! forcing fields.  Unused fields have NULL ptrs.
  logical,                            intent(in) :: set_cor  !< A flag to indicate whether to set the corrective
                                                             !! fluxes (and update the slowly varying part of eta_cor)
                                                             !! (.true.) or whether to incrementally update the
                                                             !! corrective fluxes.
  real,                               intent(in) :: dt_therm !< The thermodynamic time step, in s.
  real,                               intent(in) :: dt_since_therm !< The elapsed time since mass forcing was
                                                             !! applied, s.
  type(barotropic_cs),                pointer    :: CS       !< The control structure returned by a previous call
                                                             !! to barotropic_init.

  ! Local variables
  real :: h_tot(szi_(g))      ! The sum of the layer thicknesses, in H.
  real :: eta_h(szi_(g))      ! The free surface height determined from
                              ! the sum of the layer thicknesses, in H.
  real :: d_eta               ! The difference between estimates of the total
                              ! thicknesses, in H.
  real :: limit_dt            ! The fractional mass-source limit divided by the
                              ! thermodynamic time step, in s-1.
  integer :: is, ie, js, je, nz, i, j, k
  real, parameter :: frac_cor = 0.25
  real, parameter :: slow_rate = 0.125

  if (.not.associated(cs)) call mom_error(fatal, "bt_mass_source: "// &
        "Module MOM_barotropic must be initialized before it is used.")
  if (.not.cs%split) return

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

!$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,h,set_cor,CS,dt_therm,  &
!$OMP                                  fluxes,eta,dt_since_therm)               &
!$OMP                          private(eta_h,h_tot,limit_dt,d_eta)
  do j=js,je
    do i=is,ie ; h_tot(i) = h(i,j,1) ; enddo
    if (gv%Boussinesq) then
      do i=is,ie ; eta_h(i) = h(i,j,1) - g%bathyT(i,j)*gv%m_to_H ; enddo
    else
      do i=is,ie ; eta_h(i) = h(i,j,1) ; enddo
    endif
    do k=2,nz ; do i=is,ie
      eta_h(i) = eta_h(i) + h(i,j,k)
      h_tot(i) = h_tot(i) + h(i,j,k)
    enddo ; enddo

    if (set_cor) then
      do i=is,ie ; cs%eta_source(i,j) = 0.0 ; enddo
      if (cs%eta_source_limit > 0.0) then
        limit_dt = cs%eta_source_limit/dt_therm
        if (associated(fluxes%lprec)) then ; do i=is,ie
          cs%eta_source(i,j) = cs%eta_source(i,j) + fluxes%lprec(i,j)
        enddo ; endif
        if (associated(fluxes%fprec)) then ; do i=is,ie
          cs%eta_source(i,j) = cs%eta_source(i,j) + fluxes%fprec(i,j)
        enddo ; endif
        if (associated(fluxes%vprec)) then ; do i=is,ie
          cs%eta_source(i,j) = cs%eta_source(i,j) + fluxes%vprec(i,j)
        enddo ; endif
        if (associated(fluxes%lrunoff)) then ; do i=is,ie
          cs%eta_source(i,j) = cs%eta_source(i,j) + fluxes%lrunoff(i,j)
        enddo ; endif
        if (associated(fluxes%frunoff)) then ; do i=is,ie
          cs%eta_source(i,j) = cs%eta_source(i,j) + fluxes%frunoff(i,j)
        enddo ; endif
        if (associated(fluxes%evap)) then ; do i=is,ie
          cs%eta_source(i,j) = cs%eta_source(i,j) + fluxes%evap(i,j)
        enddo ; endif
        do i=is,ie
          cs%eta_source(i,j) = cs%eta_source(i,j)*gv%kg_m2_to_H
          if (abs(cs%eta_source(i,j)) > limit_dt * h_tot(i)) then
            cs%eta_source(i,j) = sign(limit_dt * h_tot(i), cs%eta_source(i,j))
          endif
        enddo
      endif
    endif

    if (set_cor) then
      do i=is,ie
        d_eta = eta_h(i) - (eta(i,j) - dt_since_therm*cs%eta_source(i,j))
        cs%eta_cor(i,j) = d_eta
      enddo
    else
      do i=is,ie
        d_eta = eta_h(i) - (eta(i,j) - dt_since_therm*cs%eta_source(i,j))
        cs%eta_cor(i,j) = cs%eta_cor(i,j) + d_eta
      enddo
    endif
  enddo

end subroutine bt_mass_source

!> barotropic_init initializes a number of time-invariant fields used in the
!! barotropic calculation and initializes any barotropic fields that have not
!! already been initialized.
subroutine barotropic_init(u, v, h, eta, Time, G, GV, param_file, diag, CS, &
                           restart_CS, BT_cont, tides_CSp)
  type(ocean_grid_type),            intent(inout) :: G          !< The ocean's grid structure.
  type(verticalgrid_type),          intent(in)    :: GV         !< The ocean's vertical grid structure.
  real, intent(in), dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u    !< The zonal velocity, in m s-1.
  real, intent(in), dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v    !< The meridional velocity, in m s-1.
  real, intent(in), dimension(SZI_(G),SZJ_(G),SZK_(G))  :: h    !< Layer thicknesses, in H (usually m or kg m-2).
  real, intent(in), dimension(SZI_(G),SZJ_(G))    :: eta        !< Free surface height or column mass anomaly, in
                                                                !! m or kg m-2.
  type(time_type), target,          intent(in)    :: Time       !< The current model time.
  type(param_file_type),            intent(in)    :: param_file !< A structure to parse for run-time parameters.
  type(diag_ctrl), target,          intent(inout) :: diag       !< A structure that is used to regulate diagnostic
                                                                !! output.
  type(barotropic_cs),              pointer       :: CS         !< A pointer to the control structure for this module
                                                                !! that is set in register_barotropic_restarts.
  type(mom_restart_cs),             pointer       :: restart_CS !< A pointer to the restart control structure.
  type(bt_cont_type),     optional, pointer       :: BT_cont    !< A structure with elements that describe the
                                                                !! effective open face areas as a function of
                                                                !! barotropic flow.
  type(tidal_forcing_cs), optional, pointer       :: tides_CSp  !< A pointer to the control structure of the tide
                                                                !! module.

! This include declares and sets the variable "version".
#include "version_variable.h"
  ! Local variables
  character(len=40)  :: mdl = "MOM_barotropic"  ! This module's name.
  real :: Datu(szibs_(g),szj_(g)), Datv(szi_(g),szjbs_(g))
  real :: gtot_estimate ! Summing GV%g_prime gives an upper-bound estimate for pbce.
  real :: SSH_extra     ! An estimate of how much higher SSH might get, for use
                        ! in calculating the safe external wave speed.
  real :: dtbt_input
  type(memory_size_type) :: MS
  type(group_pass_type) :: pass_static_data, pass_q_D_Cor
  type(group_pass_type) :: pass_bt_hbt_btav, pass_a_polarity
  logical :: apply_bt_drag, use_BT_cont_type
  character(len=48) :: thickness_units, flux_units
  character*(40) :: hvel_str
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  integer :: isdw, iedw, jsdw, jedw
  integer :: i, j, k
  integer :: wd_halos(2), bt_halo_sz
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  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
  ms%isdw = g%isd ; ms%iedw = g%ied ; ms%jsdw = g%jsd ; ms%jedw = g%jed

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

  cs%diag => diag ; cs%Time => time
  if (present(tides_csp)) then
    if (associated(tides_csp)) cs%tides_CSp => tides_csp
  endif

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "SPLIT", cs%split, &
                 "Use the split time stepping if true.", default=.true.)
  if (.not.cs%split) return

  call get_param(param_file, mdl, "BOUND_BT_CORRECTION", cs%bound_BT_corr, &
                 "If true, the corrective pseudo mass-fluxes into the \n"//&
                 "barotropic solver are limited to values that require \n"//&
                 "less than maxCFL_BT_cont to be accommodated.",default=.false.)
  call get_param(param_file, mdl, "BT_CONT_CORR_BOUNDS", cs%BT_cont_bounds, &
                 "If true, and BOUND_BT_CORRECTION is true, use the \n"//&
                 "BT_cont_type variables to set limits determined by \n"//&
                 "MAXCFL_BT_CONT on the CFL number of the velocites \n"//&
                 "that are likely to be driven by the corrective mass fluxes.", &
                 default=.true.) !, do_not_log=.not.CS%bound_BT_corr)
  call get_param(param_file, mdl, "ADJUST_BT_CONT", cs%adjust_BT_cont, &
                 "If true, adjust the curve fit to the BT_cont type \n"//&
                 "that is used by the barotropic solver to match the \n"//&
                 "transport about which the flow is being linearized.", default=.false.)
  call get_param(param_file, mdl, "GRADUAL_BT_ICS", cs%gradual_BT_ICs, &
                 "If true, adjust the initial conditions for the \n"//&
                 "barotropic solver to the values from the layered \n"//&
                 "solution over a whole timestep instead of instantly. \n"//&
                 "This is a decent approximation to the inclusion of \n"//&
                 "sum(u dh_dt) while also correcting for truncation errors.", &
                 default=.false.)
  call get_param(param_file, mdl, "BT_USE_VISC_REM_U_UH0", cs%visc_rem_u_uh0, &
                 "If true, use the viscous remnants when estimating the \n"//&
                 "barotropic velocities that were used to calculate uh0 \n"//&
                 "and vh0.  False is probably the better choice.", default=.false.)
  call get_param(param_file, mdl, "BT_USE_WIDE_HALOS", cs%use_wide_halos, &
                 "If true, use wide halos and march in during the \n"//&
                 "barotropic time stepping for efficiency.", default=.true., &
                 layoutparam=.true.)
  call get_param(param_file, mdl, "BTHALO", bt_halo_sz, &
                 "The minimum halo size for the barotropic solver.", default=0, &
                 layoutparam=.true.)
#ifdef STATIC_MEMORY_
  if ((bt_halo_sz > 0) .and. (bt_halo_sz /= bthalo_)) call mom_error(fatal, &
      "barotropic_init: Run-time values of BTHALO must agree with the \n"//&
      "macro BTHALO_ with STATIC_MEMORY_.")
  wd_halos(1) = whaloi_+nihalo_ ; wd_halos(2) = whaloj_+njhalo_
#else
  wd_halos(1) = bt_halo_sz; wd_halos(2) =  bt_halo_sz
#endif
  call log_param(param_file, mdl, "!BT x-halo", wd_halos(1), &
                 "The barotropic x-halo size that is actually used.", &
                 layoutparam=.true.)
  call log_param(param_file, mdl, "!BT y-halo", wd_halos(2), &
                 "The barotropic y-halo size that is actually used.", &
                 layoutparam=.true.)

  call get_param(param_file, mdl, "USE_BT_CONT_TYPE", use_bt_cont_type, &
               "If true, use a structure with elements that describe \n"//&
               "effective face areas from the summed continuity solver \n"//&
               "as a function the barotropic flow in coupling between \n"//&
               "the barotropic and baroclinic flow.  This is only used \n"//&
               "if SPLIT is true. \n", default=.true.)
  call get_param(param_file, mdl, "NONLINEAR_BT_CONTINUITY", &
                                cs%Nonlinear_continuity, &
                 "If true, use nonlinear transports in the barotropic \n"//&
                 "continuity equation.  This does not apply if \n"//&
                 "USE_BT_CONT_TYPE is true.", default=.false.)
  cs%Nonlin_cont_update_period = 1
  if (cs%Nonlinear_continuity) &
    call get_param(param_file, mdl, "NONLIN_BT_CONT_UPDATE_PERIOD", &
                                  cs%Nonlin_cont_update_period, &
                 "If NONLINEAR_BT_CONTINUITY is true, this is the number \n"//&
                 "of barotropic time steps between updates to the face \n"//&
                 "areas, or 0 to update only before the barotropic stepping.",&
                 units="nondim", default=1)
  call get_param(param_file, mdl, "BT_MASS_SOURCE_LIMIT", cs%eta_source_limit, &
                 "The fraction of the initial depth of the ocean that can \n"//&
                 "be added to or removed from the bartropic solution \n"//&
                 "within a thermodynamic time step.  By default this is 0 \n"//&
                 "for no correction.", units="nondim", default=0.0)
  call get_param(param_file, mdl, "BT_PROJECT_VELOCITY", cs%BT_project_velocity,&
                 "If true, step the barotropic velocity first and project \n"//&
                 "out the velocity tendancy by 1+BEBT when calculating the \n"//&
                 "transport.  The default (false) is to use a predictor \n"//&
                 "continuity step to find the pressure field, and then \n"//&
                 "to do a corrector continuity step using a weighted \n"//&
                 "average of the old and new velocities, with weights \n"//&
                 "of (1-BEBT) and BEBT.", default=.false.)

  call get_param(param_file, mdl, "DYNAMIC_SURFACE_PRESSURE", cs%dynamic_psurf, &
                 "If true, add a dynamic pressure due to a viscous ice \n"//&
                 "shelf, for instance.", default=.false.)
  if (cs%dynamic_psurf) then
    call get_param(param_file, mdl, "ICE_LENGTH_DYN_PSURF", cs%ice_strength_length, &
                 "The length scale at which the Rayleigh damping rate due \n"//&
                 "to the ice strength should be the same as if a Laplacian \n"//&
                 "were applied, if DYNAMIC_SURFACE_PRESSURE is true.", &
                 units="m", default=1.0e4)
    call get_param(param_file, mdl, "DEPTH_MIN_DYN_PSURF", cs%Dmin_dyn_psurf, &
                  "The minimum depth to use in limiting the size of the \n"//&
                  "dynamic surface pressure for stability, if \n"//&
                  "DYNAMIC_SURFACE_PRESSURE is true..", units="m", &
                  default=1.0e-6)
    call get_param(param_file, mdl, "CONST_DYN_PSURF", cs%const_dyn_psurf, &
                 "The constant that scales the dynamic surface pressure, \n"//&
                 "if DYNAMIC_SURFACE_PRESSURE is true.  Stable values \n"//&
                 "are < ~1.0.", units="nondim", default=0.9)
  endif

  call get_param(param_file, mdl, "TIDES", cs%tides, &
                 "If true, apply tidal momentum forcing.", default=.false.)
  call get_param(param_file, mdl, "SADOURNY", cs%Sadourny, &
                 "If true, the Coriolis terms are discretized with the \n"//&
                 "Sadourny (1975) energy conserving scheme, otherwise \n"//&
                 "the Arakawa & Hsu scheme is used.  If the internal \n"//&
                 "deformation radius is not resolved, the Sadourny scheme \n"//&
                 "should probably be used.", default=.true.)

  call get_param(param_file, mdl, "BT_THICK_SCHEME", hvel_str, &
                 "A string describing the scheme that is used to set the \n"//&
                 "open face areas used for barotropic transport and the \n"//&
                 "relative weights of the accelerations. Valid values are:\n"//&
                 "\t ARITHMETIC - arithmetic mean layer thicknesses \n"//&
                 "\t HARMONIC - harmonic mean layer thicknesses \n"//&
                 "\t HYBRID (the default) - use arithmetic means for \n"//&
                 "\t    layers above the shallowest bottom, the harmonic \n"//&
                 "\t    mean for layers below, and a weighted average for \n"//&
                 "\t    layers that straddle that depth \n"//&
                 "\t FROM_BT_CONT - use the average thicknesses kept \n"//&
                 "\t    in the h_u and h_v fields of the BT_cont_type", &
                 default=bt_cont_string)
  select case (hvel_str)
    case (hybrid_string) ; cs%hvel_scheme = hybrid
    case (harmonic_string) ; cs%hvel_scheme = harmonic
    case (arithmetic_string) ; cs%hvel_scheme = arithmetic
    case (bt_cont_string) ; cs%hvel_scheme = from_bt_cont
    case default
      call mom_mesg('barotropic_init: BT_THICK_SCHEME ="'//trim(hvel_str)//'"', 0)
      call mom_error(fatal, "barotropic_init: Unrecognized setting "// &
            "#define BT_THICK_SCHEME "//trim(hvel_str)//" found in input file.")
  end select
  if ((cs%hvel_scheme == from_bt_cont) .and. .not.use_bt_cont_type) &
    call mom_error(fatal, "barotropic_init: BT_THICK_SCHEME FROM_BT_CONT "//&
                           "can only be used if USE_BT_CONT_TYPE is defined.")

  call get_param(param_file, mdl, "BT_STRONG_DRAG", cs%strong_drag, &
                 "If true, use a stronger estimate of the retarding \n"//&
                 "effects of strong bottom drag, by making it implicit \n"//&
                 "with the barotropic time-step instead of implicit with \n"//&
                 "the baroclinic time-step and dividing by the number of \n"//&
                 "barotropic steps.", default=.false.)

  call get_param(param_file, mdl, "CLIP_BT_VELOCITY", cs%clip_velocity, &
                 "If true, limit any velocity components that exceed \n"//&
                 "CFL_TRUNCATE.  This should only be used as a desperate \n"//&
                 "debugging measure.", default=.false.)
  call get_param(param_file, mdl, "CFL_TRUNCATE", cs%CFL_trunc, &
                 "The value of the CFL number that will cause velocity \n"//&
                 "components to be truncated; instability can occur past 0.5.", &
                 units="nondim", default=0.5, do_not_log=.not.cs%clip_velocity)
  call get_param(param_file, mdl, "MAXVEL", cs%maxvel, &
                 "The maximum velocity allowed before the velocity \n"//&
                 "components are truncated.", units="m s-1", default=3.0e8, &
                 do_not_log=.not.cs%clip_velocity)
  call get_param(param_file, mdl, "MAXCFL_BT_CONT", cs%maxCFL_BT_cont, &
                 "The maximum permitted CFL number associated with the \n"//&
                 "barotropic accelerations from the summed velocities \n"//&
                 "times the time-derivatives of thicknesses.", units="nondim", &
                 default=0.25)

  call get_param(param_file, mdl, "DT_BT_FILTER", cs%dt_bt_filter, &
                 "A time-scale over which the barotropic mode solutions \n"//&
                 "are filtered, in seconds if positive, or as a fraction \n"//&
                 "of DT if negative. When used this can never be taken to \n"//&
                 "be longer than 2*dt.  Set this to 0 to apply no filtering.", &
                 units="sec or nondim", default=-0.25)
  call get_param(param_file, mdl, "G_BT_EXTRA", cs%G_extra, &
                 "A nondimensional factor by which gtot is enhanced.", &
                 units="nondim", default=0.0)
  call get_param(param_file, mdl, "SSH_EXTRA", ssh_extra, &
                 "An estimate of how much higher SSH might get, for use \n"//&
                 "in calculating the safe external wave speed. The \n"//&
                 "default is the minimum of 10 m or 5% of MAXIMUM_DEPTH.", &
                 units="m", default=min(10.0,0.05*g%max_depth))

  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_BT", cs%debug_bt, &
                 "If true, write out verbose debugging data within the \n"//&
                 "barotropic time-stepping loop. The data volume can be \n"//&
                 "quite large if this is true.", default=cs%debug)

  cs%linearized_BT_PV = .true.
  call get_param(param_file, mdl, "BEBT", cs%bebt, &
                 "BEBT determines whether the barotropic time stepping \n"//&
                 "uses the forward-backward time-stepping scheme or a \n"//&
                 "backward Euler scheme. BEBT is valid in the range from \n"//&
                 "0 (for a forward-backward treatment of nonrotating \n"//&
                 "gravity waves) to 1 (for a backward Euler treatment). \n"//&
                 "In practice, BEBT must be greater than about 0.05.", &
                 units="nondim", default=0.1)
  call get_param(param_file, mdl, "DTBT", cs%dtbt, &
                 "The barotropic time step, in s. DTBT is only used with \n"//&
                 "the split explicit time stepping. To set the time step \n"//&
                 "automatically based the maximum stable value use 0, or \n"//&
                 "a negative value gives the fraction of the stable value. \n"//&
                 "Setting DTBT to 0 is the same as setting it to -0.98. \n"//&
                 "The value of DTBT that will actually be used is an \n"//&
                 "integer fraction of DT, rounding down.", units="s or nondim",&
                 default = -0.98)

  ! Initialize a version of the MOM domain that is specific to the barotropic solver.
  call clone_mom_domain(g%Domain, cs%BT_Domain, min_halo=wd_halos, symmetric=.true.)
#ifdef STATIC_MEMORY_
  if (wd_halos(1) /= whaloi_+nihalo_) call mom_error(fatal, "barotropic_init: "//&
          "Barotropic x-halo sizes are incorrectly resized with STATIC_MEMORY_.")
  if (wd_halos(2) /= whaloj_+njhalo_) call mom_error(fatal, "barotropic_init: "//&
          "Barotropic y-halo sizes are incorrectly resized with STATIC_MEMORY_.")
#else
  if (bt_halo_sz > 0) then
    if (wd_halos(1) > bt_halo_sz) &
      call mom_mesg("barotropic_init: barotropic x-halo size increased.", 3)
    if (wd_halos(2) > bt_halo_sz) &
      call mom_mesg("barotropic_init: barotropic y-halo size increased.", 3)
  endif
#endif

  cs%isdw = g%isc-wd_halos(1) ; cs%iedw = g%iec+wd_halos(1)
  cs%jsdw = g%jsc-wd_halos(2) ; cs%jedw = g%jec+wd_halos(2)
  isdw = cs%isdw ; iedw = cs%iedw ; jsdw = cs%jsdw ; jedw = cs%jedw

  alloc_(cs%frhatu(isdb:iedb,jsd:jed,nz)) ; alloc_(cs%frhatv(isd:ied,jsdb:jedb,nz))
  alloc_(cs%eta_source(isd:ied,jsd:jed)) ; alloc_(cs%eta_cor(isd:ied,jsd:jed))
  if (cs%bound_BT_corr) then
    alloc_(cs%eta_cor_bound(isd:ied,jsd:jed)) ; cs%eta_cor_bound(:,:) = 0.0
  endif
  alloc_(cs%IDatu(isdb:iedb,jsd:jed)) ; alloc_(cs%IDatv(isd:ied,jsdb:jedb))

  alloc_(cs%ua_polarity(isdw:iedw,jsdw:jedw))
  alloc_(cs%va_polarity(isdw:iedw,jsdw:jedw))

  cs%frhatu(:,:,:) = 0.0 ; cs%frhatv(:,:,:) = 0.0
  cs%eta_source(:,:) = 0.0 ; cs%eta_cor(:,:) = 0.0
  cs%IDatu(:,:) = 0.0 ; cs%IDatv(:,:) = 0.0

  cs%ua_polarity(:,:) = 1.0 ; cs%va_polarity(:,:) = 1.0
  call create_group_pass(pass_a_polarity, cs%ua_polarity, cs%va_polarity, cs%BT_domain, to_all, agrid)
  call do_group_pass(pass_a_polarity, cs%BT_domain)

  if (use_bt_cont_type) &
    call alloc_bt_cont_type(bt_cont, g, (cs%hvel_scheme == from_bt_cont))

  if (cs%debug) then ! Make a local copy of loop ranges for chksum calls
    allocate(cs%debug_BT_HI)
    cs%debug_BT_HI%isc=g%isc
    cs%debug_BT_HI%iec=g%iec
    cs%debug_BT_HI%jsc=g%jsc
    cs%debug_BT_HI%jec=g%jec
    cs%debug_BT_HI%IscB=g%isc-1
    cs%debug_BT_HI%IecB=g%iec
    cs%debug_BT_HI%JscB=g%jsc-1
    cs%debug_BT_HI%JecB=g%jec
    cs%debug_BT_HI%isd=cs%isdw
    cs%debug_BT_HI%ied=cs%iedw
    cs%debug_BT_HI%jsd=cs%jsdw
    cs%debug_BT_HI%jed=cs%jedw
    cs%debug_BT_HI%IsdB=cs%isdw-1
    cs%debug_BT_HI%IedB=cs%iedw
    cs%debug_BT_HI%JsdB=cs%jsdw-1
    cs%debug_BT_HI%JedB=cs%jedw
  endif

  ! IareaT, IdxCu, and IdyCv need to be allocated with wide halos.
  alloc_(cs%IareaT(cs%isdw:cs%iedw,cs%jsdw:cs%jedw)) ; cs%IareaT(:,:) = 0.0
  alloc_(cs%bathyT(cs%isdw:cs%iedw,cs%jsdw:cs%jedw)) ; cs%bathyT(:,:) = gv%Angstrom_z !### Change to 0.0?
  alloc_(cs%IdxCu(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw)) ; cs%IdxCu(:,:) = 0.0
  alloc_(cs%IdyCv(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw)) ; cs%IdyCv(:,:) = 0.0
  alloc_(cs%dy_Cu(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw)) ; cs%dy_Cu(:,:) = 0.0
  alloc_(cs%dx_Cv(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw)) ; cs%dx_Cv(:,:) = 0.0
  do j=g%jsd,g%jed ; do i=g%isd,g%ied
    cs%IareaT(i,j) = g%IareaT(i,j)
    cs%bathyT(i,j) = g%bathyT(i,j)
  enddo ; enddo
  ! Note: G%IdxCu & G%IdyCv may be smaller than CS%IdxCu & CS%IdyCv, even without
  !   wide halos.
  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
    cs%IdxCu(i,j) = g%IdxCu(i,j) ; cs%dy_Cu(i,j) = g%dy_Cu(i,j)
  enddo ; enddo
  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
    cs%IdyCv(i,j) = g%IdyCv(i,j) ; cs%dx_Cv(i,j) = g%dx_Cv(i,j)
  enddo ; enddo
  call create_group_pass(pass_static_data, cs%IareaT, cs%BT_domain, to_all)
  call create_group_pass(pass_static_data, cs%bathyT, cs%BT_domain, to_all)
  call create_group_pass(pass_static_data, cs%IdxCu, cs%IdyCv, cs%BT_domain, &
                         to_all+scalar_pair)
  call create_group_pass(pass_static_data, cs%dy_Cu, cs%dx_Cv, cs%BT_domain, &
                         to_all+scalar_pair)
  call do_group_pass(pass_static_data, cs%BT_domain)

  if (cs%linearized_BT_PV) then
    alloc_(cs%q_D(cs%isdw-1:cs%iedw,cs%jsdw-1:cs%jedw))
    alloc_(cs%D_u_Cor(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw))
    alloc_(cs%D_v_Cor(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw))
    cs%q_D(:,:) = 0.0 ; cs%D_u_Cor(:,:) = 0.0 ; cs%D_v_Cor(:,:) = 0.0
    do j=js,je ; do i=is-1,ie
      cs%D_u_Cor(i,j) = 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      cs%D_v_Cor(i,j) = 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j))
    enddo ; enddo
    do j=js-1,je ; do i=is-1,ie
      if (g%mask2dT(i,j)+g%mask2dT(i,j+1)+g%mask2dT(i+1,j)+g%mask2dT(i+1,j+1)>0.) then
        cs%q_D(i,j) = 0.25 * g%CoriolisBu(i,j) * &
           ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
           ((g%areaT(i,j) * g%bathyT(i,j) + g%areaT(i+1,j+1) * g%bathyT(i+1,j+1)) + &
            (g%areaT(i+1,j) * g%bathyT(i+1,j) + g%areaT(i,j+1) * g%bathyT(i,j+1)))
      else ! All four h points are masked out so q_D(I,J) will is meaningless
        cs%q_D(i,j) = 0.
      endif
    enddo ; enddo
    ! With very wide halos, q and D need to be calculated on the available data
    ! domain and then updated onto the full computational domain.
    call create_group_pass(pass_q_d_cor, cs%q_D, cs%BT_Domain, to_all, position=corner)
    call create_group_pass(pass_q_d_cor, cs%D_u_Cor, cs%D_v_Cor, cs%BT_Domain, &
                           to_all+scalar_pair)
    call do_group_pass(pass_q_d_cor, cs%BT_Domain)
  endif

  ! Estimate the maximum stable barotropic time step.
  dtbt_input = cs%dtbt
  cs%dtbt_fraction = 0.98 ; if (cs%dtbt < 0.0) cs%dtbt_fraction = -cs%dtbt
  gtot_estimate = 0.0
  do k=1,g%ke ; gtot_estimate = gtot_estimate + gv%g_prime(k) ; enddo
  call set_dtbt(g, gv, cs, gtot_est = gtot_estimate, ssh_add = ssh_extra)
  if (dtbt_input > 0.0) cs%dtbt = dtbt_input

  call log_param(param_file, mdl, "DTBT as used", cs%dtbt)
  call log_param(param_file, mdl, "estimated maximum DTBT", cs%dtbt_max)

  ! ubtav, vbtav, ubt_IC, vbt_IC, uhbt_IC, and vhbt_IC are allocated and
  ! initialized in register_barotropic_restarts.

  if (gv%Boussinesq) then
    thickness_units = "meter" ; flux_units = "meter3 second-1"
  else
    thickness_units = "kilogram meter-2" ; flux_units = "kilogram second-1"
  endif

  cs%id_PFu_bt = register_diag_field('ocean_model', 'PFuBT', diag%axesCu1, time, &
      'Zonal Anomalous Barotropic Pressure Force Force Acceleration', 'meter second-2')
  cs%id_PFv_bt = register_diag_field('ocean_model', 'PFvBT', diag%axesCv1, time, &
      'Meridional Anomalous Barotropic Pressure Force Acceleration', 'meter second-2')
  cs%id_Coru_bt = register_diag_field('ocean_model', 'CoruBT', diag%axesCu1, time, &
      'Zonal Barotropic Coriolis Acceleration', 'meter second-2')
  cs%id_Corv_bt = register_diag_field('ocean_model', 'CorvBT', diag%axesCv1, time, &
      'Meridional Barotropic Coriolis Acceleration', 'meter second-2')
  cs%id_uaccel = register_diag_field('ocean_model', 'u_accel_bt', diag%axesCu1, time, &
      'Barotropic zonal acceleration', 'meter second-2')
  cs%id_vaccel = register_diag_field('ocean_model', 'v_accel_bt', diag%axesCv1, time, &
      'Barotropic meridional acceleration', 'meter second-2')
  cs%id_ubtforce = register_diag_field('ocean_model', 'ubtforce', diag%axesCu1, time, &
      'Barotropic zonal acceleration from baroclinic terms', 'meter second-2')
  cs%id_vbtforce = register_diag_field('ocean_model', 'vbtforce', diag%axesCv1, time, &
      'Barotropic meridional acceleration from baroclinic terms', 'meter second-2')

  cs%id_eta_bt = register_diag_field('ocean_model', 'eta_bt', diag%axesT1, time, &
      'Barotropic end SSH', thickness_units)
  cs%id_ubt = register_diag_field('ocean_model', 'ubt', diag%axesCu1, time, &
      'Barotropic end zonal velocity', 'meter second-1')
  cs%id_vbt = register_diag_field('ocean_model', 'vbt', diag%axesCv1, time, &
      'Barotropic end meridional velocity', 'meter second-1')
  cs%id_eta_st = register_diag_field('ocean_model', 'eta_st', diag%axesT1, time, &
      'Barotropic start SSH', thickness_units)
  cs%id_ubt_st = register_diag_field('ocean_model', 'ubt_st', diag%axesCu1, time, &
      'Barotropic start zonal velocity', 'meter second-1')
  cs%id_vbt_st = register_diag_field('ocean_model', 'vbt_st', diag%axesCv1, time, &
      'Barotropic start meridional velocity', 'meter second-1')
  cs%id_ubtav = register_diag_field('ocean_model', 'ubtav', diag%axesCu1, time, &
      'Barotropic time-average zonal velocity', 'meter second-1')
  cs%id_vbtav = register_diag_field('ocean_model', 'vbtav', diag%axesCv1, time, &
      'Barotropic time-average meridional velocity', 'meter second-1')
  cs%id_eta_cor = register_diag_field('ocean_model', 'eta_cor', diag%axesT1, time, &
      'Corrective mass flux', 'meter second-1')
  cs%id_visc_rem_u = register_diag_field('ocean_model', 'visc_rem_u', diag%axesCuL, time, &
      'Viscous remnant at u', 'Nondim')
  cs%id_visc_rem_v = register_diag_field('ocean_model', 'visc_rem_v', diag%axesCvL, time, &
      'Viscous remnant at v', 'Nondim')
  cs%id_gtotn = register_diag_field('ocean_model', 'gtot_n', diag%axesT1, time, &
      'gtot to North', 'm s-2')
  cs%id_gtots = register_diag_field('ocean_model', 'gtot_s', diag%axesT1, time, &
      'gtot to South', 'm s-2')
  cs%id_gtote = register_diag_field('ocean_model', 'gtot_e', diag%axesT1, time, &
      'gtot to East', 'm s-2')
  cs%id_gtotw = register_diag_field('ocean_model', 'gtot_w', diag%axesT1, time, &
      'gtot to West', 'm s-2')
  cs%id_eta_hifreq = register_diag_field('ocean_model', 'eta_hifreq', diag%axesT1, time, &
      'High Frequency Barotropic SSH', thickness_units)
  cs%id_ubt_hifreq = register_diag_field('ocean_model', 'ubt_hifreq', diag%axesCu1, time, &
      'High Frequency Barotropic zonal velocity', 'meter second-1')
  cs%id_vbt_hifreq = register_diag_field('ocean_model', 'vbt_hifreq', diag%axesCv1, time, &
      'High Frequency Barotropic meridional velocity', 'meter second-1')
  cs%id_eta_pred_hifreq = register_diag_field('ocean_model', 'eta_pred_hifreq', diag%axesT1, time, &
      'High Frequency Predictor Barotropic SSH', thickness_units)
  cs%id_uhbt_hifreq = register_diag_field('ocean_model', 'uhbt_hifreq', diag%axesCu1, time, &
      'High Frequency Barotropic zonal transport', 'meter3 second-1')
  cs%id_vhbt_hifreq = register_diag_field('ocean_model', 'vhbt_hifreq', diag%axesCv1, time, &
      'High Frequency Barotropic meridional transport', 'meter3 second-1')
  cs%id_frhatu = register_diag_field('ocean_model', 'frhatu', diag%axesCuL, time, &
      'Fractional thickness of layers in u-columns', 'Nondim')
  cs%id_frhatv = register_diag_field('ocean_model', 'frhatv', diag%axesCvL, time, &
      'Fractional thickness of layers in v-columns', 'Nondim')
  cs%id_frhatu1 = register_diag_field('ocean_model', 'frhatu1', diag%axesCuL, time, &
      'Predictor Fractional thickness of layers in u-columns', 'Nondim')
  cs%id_frhatv1 = register_diag_field('ocean_model', 'frhatv1', diag%axesCvL, time, &
      'Predictor Fractional thickness of layers in v-columns', 'Nondim')
  cs%id_uhbt = register_diag_field('ocean_model', 'uhbt', diag%axesCu1, time, &
      'Barotropic zonal transport averaged over a baroclinic step', 'meter3 second-1')
  cs%id_vhbt = register_diag_field('ocean_model', 'vhbt', diag%axesCv1, time, &
      'Barotropic meridional transport averaged over a baroclinic step', 'meter3 second-1')

  if (use_bt_cont_type) then
    cs%id_BTC_FA_u_EE = register_diag_field('ocean_model', 'BTC_FA_u_EE', diag%axesCu1, time, &
        'BTCont type far east face area', 'meter2')
    cs%id_BTC_FA_u_E0 = register_diag_field('ocean_model', 'BTC_FA_u_E0', diag%axesCu1, time, &
        'BTCont type near east face area', 'meter2')
    cs%id_BTC_FA_u_WW = register_diag_field('ocean_model', 'BTC_FA_u_WW', diag%axesCu1, time, &
        'BTCont type far west face area', 'meter2')
    cs%id_BTC_FA_u_W0 = register_diag_field('ocean_model', 'BTC_FA_u_W0', diag%axesCu1, time, &
        'BTCont type near west face area', 'meter2')
    cs%id_BTC_ubt_EE = register_diag_field('ocean_model', 'BTC_ubt_EE', diag%axesCu1, time, &
        'BTCont type far east velocity', 'meter second-1')
    cs%id_BTC_ubt_WW = register_diag_field('ocean_model', 'BTC_ubt_WW', diag%axesCu1, time, &
        'BTCont type far west velocity', 'meter second-1')
    cs%id_BTC_FA_v_NN = register_diag_field('ocean_model', 'BTC_FA_v_NN', diag%axesCv1, time, &
        'BTCont type far north face area', 'meter2')
    cs%id_BTC_FA_v_N0 = register_diag_field('ocean_model', 'BTC_FA_v_N0', diag%axesCv1, time, &
        'BTCont type near north face area', 'meter2')
    cs%id_BTC_FA_v_SS = register_diag_field('ocean_model', 'BTC_FA_v_SS', diag%axesCv1, time, &
        'BTCont type far south face area', 'meter2')
    cs%id_BTC_FA_v_S0 = register_diag_field('ocean_model', 'BTC_FA_v_S0', diag%axesCv1, time, &
        'BTCont type near south face area', 'meter2')
    cs%id_BTC_vbt_NN = register_diag_field('ocean_model', 'BTC_vbt_NN', diag%axesCv1, time, &
        'BTCont type far north velocity', 'meter second-1')
    cs%id_BTC_vbt_SS = register_diag_field('ocean_model', 'BTC_vbt_SS', diag%axesCv1, time, &
        'BTCont type far south velocity', 'meter second-1')
  endif
  cs%id_uhbt0 = register_diag_field('ocean_model', 'uhbt0', diag%axesCu1, time, &
      'Barotropic zonal transport difference', 'meter3 second-1')
  cs%id_vhbt0 = register_diag_field('ocean_model', 'vhbt0', diag%axesCv1, time, &
      'Barotropic meridional transport difference', 'meter3 second-1')

  if (cs%id_frhatu1 > 0) call safe_alloc_ptr(cs%frhatu1, isdb,iedb,jsd,jed,nz)
  if (cs%id_frhatv1 > 0) call safe_alloc_ptr(cs%frhatv1, isd,ied,jsdb,jedb,nz)

  if (.NOT.query_initialized(cs%ubtav,"ubtav",restart_cs) .or. &
      .NOT.query_initialized(cs%vbtav,"vbtav",restart_cs)) then
    call btcalc(h, g, gv, cs, may_use_default=.true.)
    cs%ubtav(:,:) = 0.0 ; cs%vbtav(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      cs%ubtav(i,j) = cs%ubtav(i,j) + cs%frhatu(i,j,k) * u(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      cs%vbtav(i,j) = cs%vbtav(i,j) + cs%frhatv(i,j,k) * v(i,j,k)
    enddo ; enddo ; enddo
  endif

  if (.NOT.query_initialized(cs%ubt_IC,"ubt_IC",restart_cs) .or. &
      .NOT.query_initialized(cs%vbt_IC,"vbt_IC",restart_cs)) then
    do j=js,je ; do i=is-1,ie ; cs%ubt_IC(i,j) = cs%ubtav(i,j) ; enddo ; enddo
    do j=js-1,je ; do i=is,ie ; cs%vbt_IC(i,j) = cs%vbtav(i,j) ; enddo ; enddo
  endif

!   Calculate other constants which are used for btstep.

  ! The following is only valid with the Boussinesq approximation.
! if (GV%Boussinesq) then
    do j=js,je ; do i=is-1,ie
      if (g%mask2dCu(i,j)>0.) then
        cs%IDatu(i,j) = g%mask2dCu(i,j) * 2.0 / (g%bathyT(i+1,j) + g%bathyT(i,j))
      else ! Both neighboring H points are masked out so IDatu(I,j) is meaningless
        cs%IDatu(i,j) = 0.
      endif
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      if (g%mask2dCv(i,j)>0.) then
        cs%IDatv(i,j) = g%mask2dCv(i,j) * 2.0 / (g%bathyT(i,j+1) + g%bathyT(i,j))
      else ! Both neighboring H points are masked out so IDatu(I,j) is meaningless
        cs%IDatv(i,j) = 0.
      endif
    enddo ; enddo
! else
!   do j=js,je ; do I=is-1,ie
!     CS%IDatu(I,j) = G%mask2dCu(I,j) * 2.0 / (GV%Rho0*(G%bathyT(i+1,j) + G%bathyT(i,j)))
!   enddo ; enddo
!   do J=js-1,je ; do i=is,ie
!     CS%IDatv(i,J) = G%mask2dCv(i,J) * 2.0 / (GV%Rho0*(G%bathyT(i,j+1) + G%bathyT(i,j)))
!   enddo ; enddo
! endif

  call find_face_areas(datu, datv, g, gv, cs, ms, halo=1)
  if (cs%bound_BT_corr) then
    ! ### Consider replacing maxvel with G%dxT(i,j) * (CS%maxCFL_BT_cont*Idt)
    ! ### and G%dyT(i,j) * (CS%maxCFL_BT_cont*Idt)
    do j=js,je ; do i=is,ie
      cs%eta_cor_bound(i,j) = gv%m_to_H * g%IareaT(i,j) * 0.1 * cs%maxvel * &
         ((datu(i-1,j) + datu(i,j)) + (datv(i,j) + datv(i,j-1)))
    enddo ; enddo
  endif

  if (.NOT.query_initialized(cs%uhbt_IC,"uhbt_IC",restart_cs) .or. &
      .NOT.query_initialized(cs%vhbt_IC,"vhbt_IC",restart_cs)) then
    do j=js,je ; do i=is-1,ie ; cs%uhbt_IC(i,j) = cs%ubtav(i,j) * datu(i,j) ; enddo ; enddo
    do j=js-1,je ; do i=is,ie ; cs%vhbt_IC(i,j) = cs%vbtav(i,j) * datv(i,j) ; enddo ; enddo
  endif

  call create_group_pass(pass_bt_hbt_btav, cs%ubt_IC, cs%vbt_IC, g%Domain)
  call create_group_pass(pass_bt_hbt_btav, cs%uhbt_IC, cs%vhbt_IC, g%Domain)
  call create_group_pass(pass_bt_hbt_btav, cs%ubtav, cs%vbtav, g%Domain)
  call do_group_pass(pass_bt_hbt_btav, g%Domain)

!  id_clock_pass = cpu_clock_id('(Ocean BT halo updates)', grain=CLOCK_ROUTINE)
  id_clock_calc_pre  = cpu_clock_id('(Ocean BT pre-calcs only)', grain=clock_routine)
  id_clock_pass_pre = cpu_clock_id('(Ocean BT pre-step halo updates)', grain=clock_routine)
  id_clock_calc = cpu_clock_id('(Ocean BT stepping calcs only)', grain=clock_routine)
  id_clock_pass_step = cpu_clock_id('(Ocean BT stepping halo updates)', grain=clock_routine)
  id_clock_calc_post = cpu_clock_id('(Ocean BT post-calcs only)', grain=clock_routine)
  id_clock_pass_post = cpu_clock_id('(Ocean BT post-step halo updates)', grain=clock_routine)
  if (dtbt_input <= 0.0) &
    id_clock_sync = cpu_clock_id('(Ocean BT global synch)', grain=clock_routine)

end subroutine barotropic_init

!> Clean up the barotropic control structure.
subroutine barotropic_end(CS)
  type(barotropic_cs), pointer :: CS  !< Control structure to clear out.
  dealloc_(cs%frhatu)   ; dealloc_(cs%frhatv)
  dealloc_(cs%IDatu)    ; dealloc_(cs%IDatv)
  dealloc_(cs%ubtav)    ; dealloc_(cs%vbtav)
  dealloc_(cs%eta_cor)  ; dealloc_(cs%eta_source)
  dealloc_(cs%ua_polarity) ; dealloc_(cs%va_polarity)
  if (cs%bound_BT_corr) then
    dealloc_(cs%eta_cor_bound)
  endif

  call destroy_bt_obc(cs%BT_OBC)

  deallocate(cs)
end subroutine barotropic_end

!> This subroutine is used to register any fields from MOM_barotropic.F90
!! that should be written to or read from the restart file.
subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS)
  type(hor_index_type),    intent(in) :: HI         !< A horizontal index type structure.
  type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters.
  type(barotropic_cs),     pointer    :: CS         !< A pointer that is set to point to the control
                                                    !! structure for this module.
  type(verticalgrid_type), intent(in) :: GV         !< The ocean's vertical grid structure.
  type(mom_restart_cs),    pointer    :: restart_CS !< A pointer to the restart control structure.

  ! Local variables
  type(vardesc) :: vd(3)
  real :: slow_rate
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed
  isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB

  if (associated(cs)) then
    call mom_error(warning, "register_barotropic_restarts called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)

  alloc_(cs%ubtav(isdb:iedb,jsd:jed))      ; cs%ubtav(:,:) = 0.0
  alloc_(cs%vbtav(isd:ied,jsdb:jedb))      ; cs%vbtav(:,:) = 0.0
  alloc_(cs%ubt_IC(isdb:iedb,jsd:jed))     ; cs%ubt_IC(:,:) = 0.0
  alloc_(cs%vbt_IC(isd:ied,jsdb:jedb))     ; cs%vbt_IC(:,:) = 0.0
  alloc_(cs%uhbt_IC(isdb:iedb,jsd:jed))    ; cs%uhbt_IC(:,:) = 0.0
  alloc_(cs%vhbt_IC(isd:ied,jsdb:jedb))    ; cs%vhbt_IC(:,:) = 0.0

  vd(2) = var_desc("ubtav","meter second-1","Time mean barotropic zonal velocity", &
                hor_grid='u', z_grid='1')
  vd(3) = var_desc("vbtav","meter second-1","Time mean barotropic meridional velocity",&
                hor_grid='v', z_grid='1')
  call register_restart_field(cs%ubtav, vd(2), .false., restart_cs)
  call register_restart_field(cs%vbtav, vd(3), .false., restart_cs)

  vd(2) = var_desc("ubt_IC", "meter second-1", &
              longname="Next initial condition for the barotropic zonal velocity", &
              hor_grid='u', z_grid='1')
  vd(3) = var_desc("vbt_IC", "meter second-1", &
              longname="Next initial condition for the barotropic meridional velocity",&
              hor_grid='v', z_grid='1')
  call register_restart_field(cs%ubt_IC, vd(2), .false., restart_cs)
  call register_restart_field(cs%vbt_IC, vd(3), .false., restart_cs)

  if (gv%Boussinesq) then
   vd(2) = var_desc("uhbt_IC", "meter3 second-1", &
                longname="Next initial condition for the barotropic zonal transport", &
                hor_grid='u', z_grid='1')
    vd(3) = var_desc("vhbt_IC", "meter3 second-1", &
                longname="Next initial condition for the barotropic meridional transport",&
                hor_grid='v', z_grid='1')
  else
    vd(2) = var_desc("uhbt_IC", "kg second-1", &
                longname="Next initial condition for the barotropic zonal transport", &
                hor_grid='u', z_grid='1')
    vd(3) = var_desc("vhbt_IC", "kg second-1", &
                longname="Next initial condition for the barotropic meridional transport",&
                hor_grid='v', z_grid='1')
  endif
  call register_restart_field(cs%uhbt_IC, vd(2), .false., restart_cs)
  call register_restart_field(cs%vhbt_IC, vd(3), .false., restart_cs)

end subroutine register_barotropic_restarts

end module mom_barotropic