MOM_legacy_barotropic.F90

Go to the documentation of this file.

module mom_legacy_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 : pass_var, pass_vector, min_across_pes, clone_mom_domain
use mom_domains, only : pass_var_start, pass_var_complete
use mom_domains, only : pass_vector_start, pass_vector_complete
use mom_domains, only : to_all, scalar_pair, agrid, corner, mom_domain_type
use mom_domains, only : pass_var_start, pass_var_complete
use mom_domains, only : pass_vector_start, pass_vector_complete
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, obc_flather
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 legacy_btcalc, legacy_bt_mass_source, legacy_btstep, legacy_barotropic_init
public legacy_barotropic_end, register_legacy_barotropic_restarts, legacy_set_dtbt

type, public :: legacy_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_) :: &
    datu_res, &     ! A nondimensional factor by which the zonal face areas
                    ! are to be rescaled to account for the effective face
                    ! areas of the layers, if rescale_D_bt is true.
                    ! Datu_res is set in btcalc.
    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_) :: &
    datv_res, &     ! A nondimensional factor by which the meridional face
                    ! areas are to be rescaled to account for the effective
                    ! face areas of the layers, if rescale_D_bt is true.
                    ! Datv_res is set in btcalc.
    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, frhatv1 ! Predictor values.

  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 0.1*MAXVEL.
  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.
  real    :: drag_amp        ! A nondimensional value (presumably between 0 and
                             ! 1) scaling the magnitude of the bottom drag
                             ! applied in the barotropic solver, 1 by default.
  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 :: rescale_d_bt    ! If true, the open face areas are rescaled by a
                             ! function of the ratio of the summed harmonic
                             ! mean thicknesses to the harmonic mean of the
                             ! summed thicknesses.
  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
                             ! exceed maxvel.  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    :: 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.
  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.

  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_datu_res = -1, id_datv_res = -1
  integer :: id_uhbt = -1, id_frhatu = -1, id_vhbt = -1, id_frhatv = -1
  integer :: id_frhatu1 = -1, id_frhatv1 = -1
end type legacy_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


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.
  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
end type bt_obc_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
logical :: apply_u_obcs, apply_v_obcs

! 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

subroutine legacy_btstep(use_fluxes, 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, uhbt_out, vhbt_out, 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.
  logical,                    intent(in)    :: use_fluxes !< A logical indicating whether velocities
                        !! (false) or fluxes (true) are used to initialize the barotropic velocities.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                              intent(in)    :: U_in       !< The initial (3-D) zonal velocity or
                     !! volume or mass fluxes,depending on flux_form, in m s-1 or m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                              intent(in)    :: V_in       !< The initial (3-D) meridional velocity
                     !! or volume/mass fluxes, depending on flux_form, in m s-1 or m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJ_(G)),          &
                              intent(in)    :: eta_in     !< The initial barotropic free surface
                                        !! height anomaly or column mass anomaly, in m or kg m-2.
  real,                                      &
                              intent(in)    :: dt         !< The time increment over which to
                                                          !! integrate, in s.
  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.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                              intent(in)    :: U_Cor      !< The (3-D) zonal- and meridional-
                        !! velocities or volume or mass fluxes used to calculate the Coriolis
                        !! terms in bc_accel_u and!! bc_accel_v, in m s-1 or m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                              intent(in)    :: V_Cor      !< The (3-D) zonal- and meridional-
                        !! velocities or volume or mass fluxes used to calculate the Coriolis
                        !! terms in bc_accel_u and bc_accel_v, in m s-1 or m3 s-1 or kg s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                              intent(out)   :: accel_layer_u !< The accelerations 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 accelerations of each layer
                                                 !! due to the barotropic calculation, in m s-2.
  real, dimension(SZI_(G),SZJ_(G)),          &
                              intent(inout) :: 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(legacy_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      :: 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 (_u) and meridional (_v) directions.
        !! Nondimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                  intent(in), optional      :: visc_rem_v !< 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 (_u) and meridional (_v) directions.
        !! Nondimensional between 0 (at the bottom) and 1 (far above the bottom).
  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.
  real, dimension(SZIB_(G),SZJ_(G)),         &
                 intent(out), optional      :: uhbt_out   !< 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), optional      :: vhbt_out   !< The barotropic meridional volume or
                      !! mass fluxes averaged through the barotropic steps, in m3 s-1 or kg s-1.
  type(ocean_obc_type),                      &
                 pointer,     optional      :: OBC        !< An open boundary condition type, which
                                   !! contains the values associated with open boundary conditions.
  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

! Arguments: use_fluxes - A logical indicating whether velocities (false) or
!                         fluxes (true) are used to initialize the barotropic
!                         velocities.
!  (in)      U_in - The initial (3-D) zonal velocity or volume or mass fluxes,
!                   depending on flux_form, in m s-1 or m3 s-1 or kg s-1.
!  (in)      V_in - The initial (3-D) meridional velocity or volume/mass fluxes,
!                   depending on flux_form, in m s-1 or m3 s-1 or kg s-1.
!  (in)      eta_in - The initial barotropic free surface height anomaly or
!                     column mass anomaly, in m or kg m-2.
!  (in)      dt - The time increment to integrate over.
!  (in)      bc_accel_u - The zonal baroclinic accelerations, in m s-2.
!  (in)      bc_accel_v - The meridional baroclinic accelerations, in m s-2.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      pbce - The baroclinic pressure anomaly in each layer
!                   due to free surface height anomalies, in m2 H-1 s-2.
!  (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 halos.
!  (in)      U_Cor - The (3-D) zonal- and meridional- velocities or volume or
!  (in)      V_Cor   mass fluxes used to calculate the Coriolis terms in
!                    bc_accel_u and bc_accel_v, in m s-1 or m3 s-1 or kg s-1.
!  (out)     accel_layer_u - The accelerations of each layer due to the
!  (out)     accel_layer_v - barotropic calculation, in m s-2.
!  (out)     eta_out - The final barotropic free surface height anomaly or
!                      column mass anomaly, in m or kg m-2.
!  (out)     uhbtav - the barotropic zonal volume or mass fluxes averaged
!                     through the barotropic steps, in m3 s-1 or kg s-1.
!  (out)     vhbtav - the barotropic meridional volume or mass fluxes averaged
!                     through the barotropic steps, in m3 s-1 or kg s-1.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 barotropic_init.
!  (in,opt)  visc_rem_u - Both the fraction of the momentum originally in a
!  (in,opt)  visc_rem_v - 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 (_u) and
!                         meridional (_v) directions.  Nondimensional between
!                         0 (at the bottom) and 1 (far above the bottom).
!  (out,opt) etaav - The free surface height or column mass averaged over the
!                    barotropic integration, in m or kg m-2.
!  (out,opt) uhbt_out - The barotropic zonal volume or mass fluxes at the end
!                     of the barotropic steps, in m3 s-1 or kg s-1.
!  (out,opt) vhbt_out - The barotropic meridional volume or mass fluxes at the
!                     end of the barotropic steps, in m3 s-1 or kg s-1.
!  (in,opt)  OBC - An open boundary condition type, which contains the values
!                  associated with open boundary conditions.
!  (in,opt)  BT_cont - A structure with elements that describe the effective
!                      open face areas as a function of barotropic flow.
!  (in,opt)  eta_PF_start - The eta field consistent with the pressure gradient
!                      at the start of the barotropic stepping, in m or kg m-2.
!  (in,opt)  taux_bot - The zonal bottom frictional stress from ocean to the
!                       seafloor, in Pa.
!  (in,opt)  tauy_bot - The meridional bottom frictional stress from ocean to
!                       the seafloor, in Pa.

!    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.
!    Depending on the value of use_fluxes, the initial velocities are determined
!  from input velocites or volume (Boussinesq) or mass (non-Boussinesq) fluxes.

  real :: ubt_Cor(szib_(g),szj_(g)) ! The barotropic velocities that had been
  real :: vbt_Cor(szi_(g),szjb_(g)) ! use 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, if it is present. ND.
    tmp_u         ! A temporary array at u points.
  real, dimension(SZI_(G),SZJB_(G)) :: &
    av_rem_v, &   ! The weighted average of visc_rem_u, if it is present. ND.
    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_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_ref_u, &  ! The zonal barotropic Coriolis acceleration due
                  ! to the reference velocities, 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_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_ref_v, &  ! The meridional barotropic Coriolis acceleration due
                  ! to the reference velocities, 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 m or kg m-2.
    eta_pred      ! A predictor value of eta, in 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 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 :: 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 :: vel_trans   ! The combination of the previous and current velocity
                      ! that does the mass transport, 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 :: Cor, gradP  ! The Coriolis and pressure gradient accelerations, m s-1.
  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_visc_rem, use_BT_cont
  logical :: 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, 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.
  integer :: nfilter

  logical :: apply_OBCs, apply_OBC_flather
  type(bt_obc_type) :: BT_OBC  ! A structure with all of this module's fields
                               ! for applying open boundary conditions.
  type(memory_size_type) :: MS
  character(len=200) :: mesg
  integer :: pid_ubt, pid_eta, pid_e_anom, pid_etaav, pid_uhbtav, pid_ubtav
  integer :: pid_q, pid_eta_PF, pid_dyn_coef_eta, pid_eta_src
  integer :: pid_DCor_u, pid_Datu_res, pid_tmp_u, pid_gtot_E, pid_gtot_W
  integer :: pid_bt_rem_u, pid_Datu, pid_BT_force_u, pid_Cor_ref
  integer :: pid_eta_PF_1, pid_d_eta_PF, pid_uhbt0
  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

  if (.not.associated(cs)) call mom_error(fatal, &
      "legacy_btstep: Module MOM_legacy_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)
  use_visc_rem = present(visc_rem_u)
  if ((use_visc_rem) .neqv. present(visc_rem_v)) call mom_error(fatal, &
      "btstep: Either both visc_rem_u and visc_rem_v or neither"// &
       " one must be present in call to btstep.")
  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, &
      "legacy_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, &
      "legacy_btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.")
  if (add_uh0 .and. use_fluxes) call mom_error(warning, &
      "legacy_btstep: with use_fluxes, add_uh0 does nothing!")
  if (use_fluxes) add_uh0 = .false.

  ! 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. ; apply_u_obcs = .false. ; apply_v_obcs = .false.
  apply_obc_flather = .false.
  if (present(obc)) then ; if (associated(obc)) then
    apply_u_obcs = obc%Flather_u_BCs_exist_globally .or. obc%specified_u_BCs_exist_globally
    apply_v_obcs = obc%Flather_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_obcs = obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. apply_obc_flather

    if (apply_obc_flather .and. .not.gv%Boussinesq) call mom_error(fatal, &
      "legacy_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,'("legacy_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
  do_ave = query_averaging_enabled(cs%diag)

  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

!   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(none) shared(jsvf,jevf,isvf,ievf,q,CS)
    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(none) shared(jsvf,jevf,isvf,ievf,DCor_u,CS)
    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(none) shared(jsvf,jevf,isvf,ievf,DCor_v,CS)
    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.
    !  D here should be replaced with D+eta(Bous) or eta(non-Bous).
!$OMP parallel do default(none) shared(js,je,is,ie,DCor_u,G)
    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(none) shared(js,je,is,ie,DCor_v,G)
    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(none) shared(js,je,is,ie,q,G)
    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
      pid_q = pass_var_start(q, cs%BT_Domain, to_all, position=corner)
      pid_dcor_u = pass_vector_start(dcor_u, dcor_v, cs%BT_Domain, to_all+scalar_pair)
    else
      call pass_var(q, cs%BT_Domain, to_all, position=corner)
      call pass_vector(dcor_u, dcor_v, cs%BT_Domain, to_all+scalar_pair)
    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
    ! Start all halo updates that are ready to go.
  !###   if (use_BT_cont) call start_set_local_BT_cont_types( ... )
    if ((.not.use_bt_cont) .and. cs%rescale_D_bt .and. (ievf>ie)) 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)
      pid_datu_res = pass_vector_start(cs%Datu_res, cs%Datv_res, cs%BT_Domain, &
                                       to_all+scalar_pair)
      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
  endif

  ! Zero out various wide-halo arrays.
  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.
  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
  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.  eta_in and eta_PF_in
  ! have the correct values in their halo regions.
  if (interp_eta_pf) then
    do j=g%jsd,g%jed ; do i=g%isd,g%ied
      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
    do j=g%jsd,g%jed ; do i=g%isd,g%ied
      eta(i,j) = eta_in(i,j)
      eta_pf(i,j) = eta_pf_in(i,j)
    enddo ; enddo
  endif

  if (use_visc_rem) then
!$OMP parallel do default(none) shared(Isq,Ieq,js,je,nz,visc_rem_u,Instep,wt_u,CS) private(visc_rem)
    do k=1,nz ; do j=js-1,je+1 ; do i=isq-1,ieq+1
      ! 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(none) shared(is,ie,Jsq,Jeq,nz,visc_rem_v,Instep,wt_v,CS) private(visc_rem)
    do k=1,nz ; do j=jsq-1,jeq+1 ; do i=is-1,ie+1
      ! 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
  else
    do k=1,nz ; do j=js-1,je+1 ; do i=isq-1,ieq+1
      wt_u(i,j,k) = cs%frhatu(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=jsq-1,jeq+1 ; do i=is-1,ie+1
      wt_v(i,j,k) = cs%frhatv(i,j,k)
    enddo ; enddo ; enddo
  endif

  ! 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.
  do k=1,nz
    do j=js,je ; 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
    do j=js-1,je ; 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 pass_var_complete(pid_q, q, cs%BT_Domain, to_all, position=corner)
      call pass_vector_complete(pid_dcor_u, dcor_u, dcor_v, cs%BT_Domain, to_all+scalar_pair)
    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%rescale_D_bt .and. (ievf>ie)) then
      ! Datu_res was previously calculated in btcalc, and will be used in find_face_areas.
      ! This halo update is needed for wider halos than 1.  The complete goes here.
      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 pass_vector_complete(pid_datu_res, cs%Datu_res, cs%Datv_res, &
                                  cs%BT_Domain, to_all+scalar_pair)
      else
        call pass_vector(cs%Datu_res, cs%Datv_res, cs%BT_Domain, to_all+scalar_pair)
      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 (cs%Nonlinear_continuity) then
      call find_face_areas(datu, datv, g, gv, cs, ms, cs%rescale_D_bt, eta, 1)
    else
      call find_face_areas(datu, datv, g, gv, cs, ms, cs%rescale_D_bt, halo=1)
    endif
  endif

  ! Set up fields related to the open boundary conditions.
  if (apply_obcs) then
    call set_up_bt_obc(obc, eta, 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.
  !  ### Should IDatu here be replaced with 1/D+eta(Bous) or 1/eta(non-Bous)?
  if (use_visc_rem) then
    do j=js,je ; do i=is-1,ie
      bt_force_u(i,j) = fluxes%taux(i,j) * i_rho0*cs%IDatu(i,j)*visc_rem_u(i,j,1)
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      bt_force_v(i,j) = fluxes%tauy(i,j) * i_rho0*cs%IDatv(i,j)*visc_rem_v(i,j,1)
    enddo ; enddo
  else
    do j=js,je ; do i=is-1,ie
      bt_force_u(i,j) = fluxes%taux(i,j) * i_rho0 * cs%IDatu(i,j)
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      bt_force_v(i,j) = fluxes%tauy(i,j) * i_rho0 * cs%IDatv(i,j)
    enddo ; enddo
  endif
  if (present(taux_bot) .and. present(tauy_bot)) then
    if (associated(taux_bot) .and. associated(tauy_bot)) then
      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
      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.
  do k=1,nz ; do j=js,je ; 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
  do k=1,nz ; do j=jsq,jeq ; 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
    do j=js,je ; do i=is-1,ie ; uhbt(i,j) = 0.0 ; ubt(i,j) = 0.0 ; enddo ; enddo
    do j=js-1,je ; do i=is,ie ; vhbt(i,j) = 0.0 ; vbt(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js,je ; 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
    do k=1,nz ; do j=js-1,je ; 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
    if (use_bt_cont) then
      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
      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
      do j=js,je ; do i=is-1,ie
        uhbt0(i,j) = uhbt(i,j) - datu(i,j)*ubt(i,j)
      enddo ; enddo
      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.
  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
  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
  if (use_fluxes) then
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      uhbt(i,j) = uhbt(i,j) + u_in(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      vhbt(i,j) = vhbt(i,j) + v_in(i,j,k)
    enddo ; enddo ; enddo
    if (use_bt_cont) then
      do j=js,je ; do i=is-1,ie
        ubt(i,j) = uhbt_to_ubt(uhbt(i,j),btcl_u(i,j), guess=cs%ubt_IC(i,j))
      enddo ; enddo
      do j=js-1,je ; do i=is,ie
        vbt(i,j) = vhbt_to_vbt(vhbt(i,j),btcl_v(i,j), guess=cs%vbt_IC(i,j))
      enddo ; enddo
    else
      do j=js,je ; do i=is-1,ie
        ubt(i,j) = 0.0 ; if (datu(i,j) > 0.0) ubt(i,j) = uhbt(i,j) / datu(i,j)
      enddo ; enddo
      do j=js-1,je ; do i=is,ie
        vbt(i,j) = 0.0 ; if (datv(i,j) > 0.0) vbt(i,j) = vhbt(i,j) / datv(i,j)
      enddo ; enddo
    endif
  else
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_in(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_in(i,j,k)
    enddo ; enddo ; enddo
  endif

  if (cs%gradual_BT_ICs) then
    if (use_fluxes) then
      if (use_bt_cont) then
        do j=js,je ; do i=is-1,ie
          vel_tmp = uhbt_to_ubt(cs%uhbt_IC(i,j),btcl_u(i,j), guess=cs%ubt_IC(i,j))
          bt_force_u(i,j) = bt_force_u(i,j) + (ubt(i,j) - vel_tmp) * idt
          ubt(i,j) = vel_tmp
        enddo ; enddo
        do j=js-1,je ; do i=is,ie
          vel_tmp = vhbt_to_vbt(cs%vhbt_IC(i,j),btcl_v(i,j), guess=cs%vbt_IC(i,j))
          bt_force_v(i,j) = bt_force_v(i,j) + (vbt(i,j) - vel_tmp) * idt
          vbt(i,j) = vel_tmp
        enddo ; enddo
      else
        do j=js,je ; do i=is-1,ie
          vel_tmp = 0.0 ; if (datu(i,j) > 0.0) vel_tmp = cs%uhbt_IC(i,j) / datu(i,j)
          bt_force_u(i,j) = bt_force_u(i,j) + (ubt(i,j) - vel_tmp) * idt
          ubt(i,j) = vel_tmp
        enddo ; enddo
        do j=js-1,je ; do i=is,ie
          vel_tmp = 0.0 ; if (datv(i,j) > 0.0) vel_tmp = cs%vhbt_IC(i,j) / datv(i,j)
          bt_force_v(i,j) = bt_force_v(i,j) + (vbt(i,j) - vel_tmp) * idt
          vbt(i,j) = vel_tmp
        enddo ; enddo
      endif
    else
      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
      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
  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
      pid_tmp_u = pass_vector_start(tmp_u, tmp_v, g%Domain)
    else
      call pass_vector(tmp_u, tmp_v, g%Domain, complete=.true.)
      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)
    pid_gtot_e = pass_vector_start(gtot_e, gtot_n, cs%BT_Domain, &
                     to_all+scalar_pair, agrid, complete=.false.)
    pid_gtot_w = pass_vector_start(gtot_w, gtot_s, cs%BT_Domain, &
                     to_all+scalar_pair, agrid)
    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(none) shared(isvf,ievf,jsvf,jevf,amer,bmer,cmer,dmer,DCor_u,q,CS)
  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(none) shared(isvf,ievf,jsvf,jevf,azon,bzon,czon,dzon,DCor_v,q,CS)
  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

  !   If they are present, use u_Cor and v_Cor as the reference values for the
  ! Coriolis terms, including the viscous remnant if it is present.
  if (use_fluxes) then
    do j=js-1,je+1 ; do i=is-1,ie ; uhbt(i,j) = 0.0 ; enddo ; enddo
    do j=js-1,je ; do i=is-1,ie+1 ; vhbt(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie
      uhbt(i,j) = uhbt(i,j) + u_cor(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=js-1,je ; do i=is-1,ie+1
      vhbt(i,j) = vhbt(i,j) + v_cor(i,j,k)
    enddo ; enddo ; enddo
    if (use_bt_cont) then
      do j=js-1,je+1 ; do i=is-1,ie
        ubt_cor(i,j) = uhbt_to_ubt(uhbt(i,j), btcl_u(i,j), guess=cs%ubtav(i,j))
      enddo ; enddo
      do j=js-1,je ; do i=is-1,ie+1
        vbt_cor(i,j) = vhbt_to_vbt(vhbt(i,j), btcl_v(i,j), guess=cs%vbtav(i,j))
      enddo ; enddo
    else
      do j=js-1,je+1 ; do i=is-1,ie
        ubt_cor(i,j) = 0.0 ; if (datu(i,j)>0.0) ubt_cor(i,j) = uhbt(i,j)/datu(i,j)
      enddo ; enddo
      do j=js-1,je ; do i=is-1,ie+1
        vbt_cor(i,j) = 0.0 ; if (datv(i,j)>0.0) vbt_cor(i,j) = vhbt(i,j)/datv(i,j)
      enddo ; enddo
    endif
  else
    do j=js-1,je+1 ; do i=is-1,ie ; ubt_cor(i,j) = 0.0 ; enddo ; enddo
    do j=js-1,je ; do i=is-1,ie+1 ; vbt_cor(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js-1,je+1 ; 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
    do k=1,nz ; do j=js-1,je ; do i=is-1,ie+1
      vbt_cor(i,j) = vbt_cor(i,j) + wt_v(i,j,k) * v_cor(i,j,k)
    enddo ; enddo ; enddo
  endif

!$OMP parallel do default(none) shared(is,ie,js,je,Cor_ref_u,azon,bzon,czon,dzon,vbt_Cor)
  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(none) shared(is,ie,js,je,Cor_ref_v,amer,bmer,cmer,dmer,ubt_Cor)
  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

! 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 pass_vector_complete(pid_tmp_u, tmp_u, tmp_v, 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 pass_vector_complete(pid_gtot_e, gtot_e, gtot_n, cs%BT_Domain, to_all+scalar_pair, agrid)
    call pass_vector_complete(pid_gtot_w, gtot_w, gtot_s, cs%BT_Domain, to_all+scalar_pair, agrid)
  else
    call pass_vector(gtot_e, gtot_n, cs%BT_Domain, to_all+scalar_pair, agrid, complete=.false.)
    call pass_vector(gtot_w, gtot_s, cs%BT_Domain, to_all+scalar_pair, agrid, complete=.true.)
  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

  ! Now start new halo updates.
  if (nonblock_setup) then
    if (.not.use_bt_cont) &
      pid_datu = pass_vector_start(datu, datv, 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)) &
      pid_bt_force_u = pass_vector_start(bt_force_u, bt_force_v, cs%BT_Domain, &
                                         complete=.false.)
    if (add_uh0) pid_uhbt0 = pass_vector_start(uhbt0, vhbt0, cs%BT_Domain)
    pid_cor_ref = pass_vector_start(cor_ref_u, cor_ref_v, 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)

  do j=js-1,je+1 ; do i=is-1,ie ; bt_rem_u(i,j) = g%mask2dCu(i,j) ; enddo ; enddo
  do j=js-1,je ; do i=is-1,ie+1 ; bt_rem_v(i,j) = g%mask2dCv(i,j) ; enddo ; enddo
  if ((use_visc_rem) .and. (cs%drag_amp > 0.0)) then
    do j=js-1,je+1 ; do i=is-1,ie ; av_rem_u(i,j) = 0.0 ; enddo ; enddo
    do j=js-1,je ; do i=is-1,ie+1 ; av_rem_v(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js-1,je+1 ; 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
    do k=1,nz ; do j=js-1,je ; do i=is-1,ie+1
      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
      do j=js-1,je+1 ; do i=is-1,ie
        bt_rem_u(i,j) = g%mask2dCu(i,j) * cs%drag_amp * &
           ((nstep * av_rem_u(i,j)) / (1.0 + (nstep-1)*av_rem_u(i,j)))
      enddo ; enddo
      do j=js-1,je ; do i=is-1,ie+1
        bt_rem_v(i,j) = g%mask2dCv(i,j) * cs%drag_amp * &
           ((nstep * av_rem_v(i,j)) / (1.0 + (nstep-1)*av_rem_v(i,j)))
      enddo ; enddo
    else
      do j=js-1,je+1 ; 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) * cs%drag_amp * (av_rem_u(i,j)**instep)
      enddo ; enddo
      do j=js-1,je ; do i=is-1,ie+1
        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) * cs%drag_amp * (av_rem_v(i,j)**instep)
      enddo ; enddo
    endif
  endif

  ! Zero out the arrays for various time-averaged quantities.
  if (find_etaav) then ; 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 ; do j=jsvf-1,jevf+1 ; do i=isvf-1,ievf+1
    eta_wtd(i,j) = 0.0
  enddo ; enddo ; endif
  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
  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.
  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 ; 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
  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

  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 ; 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
    if (cs%dynamic_psurf) pid_dyn_coef_eta = &
      pass_var_start(dyn_coef_eta, cs%BT_Domain, complete=.false.)
    if (interp_eta_pf) then
      pid_eta_pf_1 = pass_var_start(eta_pf_1, cs%BT_Domain, complete=.false.)
      pid_d_eta_pf = pass_var_start(d_eta_pf, cs%BT_Domain, complete=.false.)
    else
      if ((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) &
        pid_eta_pf = pass_var_start(eta_pf, cs%BT_Domain, complete=.false.)
    endif
    pid_eta_src = pass_var_start(eta_src, cs%BT_Domain)

    ! The following halo update is not needed without wide halos.  RWH
    if (ievf > ie) &
      pid_bt_rem_u = pass_vector_start(bt_rem_u, bt_rem_v, cs%BT_Domain, &
                                       to_all+scalar_pair)
  else
    if (interp_eta_pf) then
      call pass_var(eta_pf_1, cs%BT_Domain, complete=.false.)
      call pass_var(d_eta_pf, cs%BT_Domain, complete=.false.)
    else
      !   eta_PF_in had correct values in its halos, so only update eta_PF with
      ! extra-wide halo arrays.  This could have started almost immediately.
      if ((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) &
        call pass_var(eta_pf, cs%BT_Domain, complete=.false.)
    endif
    if (cs%dynamic_psurf) call pass_var(dyn_coef_eta, cs%BT_Domain, complete=.false.)
    ! The following halo update is not needed without wide halos.  RWH
    if (ievf > ie) &
      call pass_vector(bt_rem_u, bt_rem_v, cs%BT_Domain, to_all+scalar_pair)
    if (.not.use_bt_cont) &
      call pass_vector(datu, datv, cs%BT_Domain, to_all+scalar_pair)
    if (((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) .or. (isq <= is-1) .or. (jsq <= js-1)) &
      call pass_vector(bt_force_u, bt_force_v, cs%BT_Domain, complete=.false.)
    if (g%nonblocking_updates) then ! Passing needs to be completed now.
      call pass_var(eta_src, cs%BT_Domain, complete=.true.)
      if (add_uh0) call pass_vector(uhbt0, vhbt0, cs%BT_Domain, complete=.false.)
      call pass_vector(cor_ref_u, cor_ref_v, cs%BT_Domain, complete=.true.)
    else
      call pass_var(eta_src, cs%BT_Domain, complete=.false.)
      if (add_uh0) call pass_vector(uhbt0, vhbt0, cs%BT_Domain, complete=.false.)
      call pass_vector(cor_ref_u, cor_ref_v, cs%BT_Domain, complete=.false.)
    endif
  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 pass_vector_complete(pid_datu, datu, datv, 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 pass_vector_complete(pid_bt_force_u, bt_force_u, bt_force_v, cs%BT_Domain)
    if (add_uh0) call pass_vector_complete(pid_uhbt0, uhbt0, vhbt0, cs%BT_Domain)
    call pass_vector_complete(pid_cor_ref, cor_ref_u, cor_ref_v, cs%BT_Domain)

    if (cs%dynamic_psurf) &
      call pass_var_complete(pid_dyn_coef_eta, dyn_coef_eta, cs%BT_Domain)
    if (interp_eta_pf) then
      call pass_var_complete(pid_eta_pf_1, eta_pf_1, cs%BT_Domain)
      call pass_var_complete(pid_d_eta_pf, d_eta_pf, cs%BT_Domain)
    else
      if ((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) &
        call pass_var_complete(pid_eta_pf, eta_pf, cs%BT_Domain)
    endif
    call pass_var_complete(pid_eta_src, eta_src, cs%BT_Domain)

    if (ievf > ie) &
      call pass_vector_complete(pid_bt_rem_u, bt_rem_u, bt_rem_v, cs%BT_Domain,&
      to_all+scalar_pair)

    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)
    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,haloshift=1)
    call uvchksum("BT frhat[uv]", cs%frhatu, cs%frhatv, g%HI, haloshift=1)
    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)
    if (use_visc_rem) then
      call uvchksum("BT visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI, haloshift=1)
    endif
  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)

  ! 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 (abs(ubt(i,j)) > cs%maxvel) then
          ! Add some error handling later.
          ubt(i,j) = sign(cs%maxvel,ubt(i,j))
        endif
      enddo ; enddo
      do j=jsv-1,jev ; do i=isv,iev
        if (abs(vbt(i,j)) > cs%maxvel) then
          ! Add some error handling later.
          vbt(i,j) = sign(cs%maxvel,vbt(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)
      if (g%nonblocking_updates) then
        pid_ubt = pass_vector_start(ubt, vbt, cs%BT_Domain)
        pid_eta = pass_var_start(eta, cs%BT_Domain)
        call pass_vector_complete(pid_ubt, ubt, vbt, cs%BT_Domain)
        call pass_var_complete(pid_eta, eta, cs%BT_Domain)
      else
        call pass_var(eta, cs%BT_Domain)
        call pass_vector(ubt, vbt, cs%BT_Domain)
      endif
      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, cs%rescale_D_bt, eta, 1+iev-ie)
    endif

    if (cs%dynamic_psurf .or. .not.project_velocity) then
      if (use_bt_cont) then
!$OMP parallel do default(none) shared(isv,iev,jsv,jev,uhbt,ubt,BTCL_u,uhbt0)
        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
        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
        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
!$OMP parallel do default(none) shared(isv,iev,jsv,jev,eta_pred,eta,eta_src,dtbt,&
!$OMP                                  CS,Datu,ubt,uhbt0,Datv,vhbt0,vbt)
        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 ; 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 ; 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.
      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) then
!$OMP parallel do default(none) shared(isv,iev,jsv,jev,ubt_old,ubt)
      do j=jsv,jev ; do i=isv-2,iev+1
        ubt_old(i,j) = ubt(i,j)
      enddo; enddo
!$OMP parallel do default(none) shared(isv,iev,jsv,jev,vbt_old,vbt)
      do j=jsv-2,jev+1 ; do i=isv,iev
        vbt_old(i,j) = vbt(i,j)
      enddo ; enddo
    endif

!$OMP parallel default(none) shared(isv,iev,jsv,jev,G,amer,ubt,cmer,bmer,dmer,eta_PF_BT, &
!$OMP                               eta_PF,gtot_N,gtot_S,dgeo_de,CS,p_surf_dyn,n,        &
!$OMP                               v_accel_bt,wt_accel,PFv_bt_sum,wt_accel2,            &
!$OMP                               Corv_bt_sum,BT_OBC,vbt,bt_rem_v,BT_force_v,vhbt,     &
!$OMP                               Cor_ref_v,find_PF,find_Cor,apply_v_OBCs,dtbt,        &
!$OMP                               project_velocity,be_proj,bebt,use_BT_cont,BTCL_v,    &
!$OMP                               vhbt0,Datv,vbt_sum,wt_trans,vhbt_sum,vbt_wtd,wt_vel, &
!$OMP                               azon,bzon,czon,dzon,Cor_ref_u,gtot_E,gtot_W,OBC,     &
!$OMP                               u_accel_bt,PFu_bt_sum,Coru_bt_sum,apply_u_OBCs,      &
!$OMP                               bt_rem_u,BT_force_u,uhbt,BTCL_u,uhbt0,Datu,ubt_sum,  &
!$OMP                               uhbt_sum,ubt_wtd)                                    &
!$OMP                       private(Cor, gradP, vel_prev, vel_trans )
    if (mod(n+g%first_direction,2)==1) then
      ! On odd-steps, update v first.
!$OMP do
      do j=jsv-1,jev ; do i=isv-1,iev+1
        cor = -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)
        gradp = ((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)
        if (cs%dynamic_psurf) &
          gradp = gradp + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
        v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor + gradp)
        if (find_pf)  pfv_bt_sum(i,j)  = pfv_bt_sum(i,j) + wt_accel2(n) * gradp
        if (find_cor) corv_bt_sum(i,j) = corv_bt_sum(i,j) + wt_accel2(n) * cor

        if (apply_v_obcs) then ; if (obc%segnum_v(i,j) /= obc_none) cycle ; endif

        vel_prev = vbt(i,j)
        vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
                        dtbt * ((bt_force_v(i,j) + cor) + gradp))
        if (project_velocity) then
          vel_trans = (1.0 + be_proj)*vbt(i,j) - be_proj*vel_prev
        else
          vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
        endif
        if (use_bt_cont) then
          vhbt(i,j) = find_vhbt(vel_trans,btcl_v(i,j)) + vhbt0(i,j)
        else
          vhbt(i,j) = datv(i,j)*vel_trans + vhbt0(i,j)
        endif

        vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vel_trans
        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

!$OMP do
      do j=jsv,jev ; do i=isv-1,iev
        cor = ((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)
        gradp = ((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)
        if (cs%dynamic_psurf) &
          gradp = gradp + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
        u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor + gradp)
        if (find_pf)  pfu_bt_sum(i,j)  = pfu_bt_sum(i,j) + wt_accel2(n) * gradp
        if (find_cor) coru_bt_sum(i,j) = coru_bt_sum(i,j) + wt_accel2(n) * cor

        if (apply_u_obcs) then ; if (obc%segnum_u(i,j) /= obc_none) cycle ; endif

        vel_prev = ubt(i,j)
        ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
                        dtbt * ((bt_force_u(i,j) + cor) + gradp))
        if (project_velocity) then
          vel_trans = (1.0 + be_proj)*ubt(i,j) - be_proj*vel_prev
        else
          vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
        endif
        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

        ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * vel_trans
        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

    else
      ! On even steps, update u first.
!$OMP do
      do j=jsv-1,jev+1 ; do i=isv-1,iev
        cor = ((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)
        gradp = ((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)
        if (cs%dynamic_psurf) &
          gradp = gradp + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
        u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor + gradp)
        if (find_pf)  pfu_bt_sum(i,j)  = pfu_bt_sum(i,j) + wt_accel2(n) * gradp
        if (find_cor) coru_bt_sum(i,j) = coru_bt_sum(i,j) + wt_accel2(n) * cor

        if (apply_u_obcs) then ; if (obc%segnum_u(i,j) /= obc_none) cycle ; endif

        vel_prev = ubt(i,j)
        ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
                        dtbt * ((bt_force_u(i,j) + cor) + gradp))
        if (project_velocity) then
          vel_trans = (1.0 + be_proj)*ubt(i,j) - be_proj*vel_prev
        else
          vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
        endif
        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

        ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * vel_trans
        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

!$OMP do
      do j=jsv-1,jev ; do i=isv,iev
        cor = -1.0*((amer(i-1,j) * ubt(i-1,j) + bmer(i,j) * ubt(i,j)) + &
                (cmer(i,j+1) * ubt(i,j+1) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
        gradp = ((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)
        if (cs%dynamic_psurf) &
          gradp = gradp + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
        v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor + gradp)
        if (find_pf)  pfv_bt_sum(i,j)  = pfv_bt_sum(i,j) + wt_accel2(n) * gradp
        if (find_cor) corv_bt_sum(i,j) = corv_bt_sum(i,j) + wt_accel2(n) * cor

        if (apply_v_obcs) then ; if (obc%segnum_v(i,j) /= obc_none) cycle ; endif

        vel_prev = vbt(i,j)
        vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
                        dtbt * ((bt_force_v(i,j) + cor) + gradp))
        if (project_velocity) then
          vel_trans = (1.0 + be_proj)*vbt(i,j) - be_proj*vel_prev
        else
          vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
        endif
         if (use_bt_cont) then
          vhbt(i,j) = find_vhbt(vel_trans,btcl_v(i,j)) + vhbt0(i,j)
        else
          vhbt(i,j) = datv(i,j)*vel_trans + vhbt0(i,j)
        endif

        vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vel_trans
        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
    endif
!$OMP end parallel

    if (apply_obcs) then
      call apply_velocity_obcs(obc, ubt, vbt, uhbt, vhbt, &
           ubt_trans, vbt_trans, eta, ubt_old, vbt_old, bt_obc, &
           g, ms, iev-ie, dtbt, bebt, use_bt_cont, datu, datv, btcl_u, btcl_v, &
           uhbt0, vhbt0)
      if (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 (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)
    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)
    enddo ; enddo
    if (apply_obcs) call apply_eta_obcs(obc, eta, ubt_old, vbt_old, 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

  ! 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
    pid_e_anom = pass_var_start(e_anom, g%Domain)
  else
    if (find_etaav) call pass_var(etaav, g%Domain, complete=.false.)
    call pass_var(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 (present(uhbt_out)) then
    uhbt_out(:,:) = 0.0
    if (use_bt_cont) then ; do j=js,je ; do i=is-1,ie
      uhbt_out(i,j) = find_uhbt(ubt_wtd(i,j), btcl_u(i,j)) + uhbt0(i,j)
    enddo ; enddo ; else ; do j=js,je ; do i=is-1,ie
      uhbt_out(i,j) = ubt_wtd(i,j) * datu(i,j) + uhbt0(i,j)
    enddo ; enddo ; endif

    vhbt_out(:,:) = 0.0
    if (use_bt_cont) then ; do j=js-1,je ; do i=is,ie
      vhbt_out(i,j) = find_vhbt(vbt_wtd(i,j), btcl_v(i,j)) + vhbt0(i,j)
    enddo ; enddo ; else ; do j=js-1,je ; do i=is,ie
      vhbt_out(i,j) = vbt_wtd(i,j) * datv(i,j) + vhbt0(i,j)
    enddo ; enddo ; endif
  endif

  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 pass_var_complete(pid_e_anom, e_anom, g%Domain)

    if (find_etaav) pid_etaav = pass_var_start(etaav, g%Domain)
    pid_uhbtav = pass_vector_start(uhbtav, vhbtav, g%Domain, complete=.false.)
    pid_ubtav = pass_vector_start(cs%ubtav, cs%vbtav, g%Domain)
  else
    call pass_vector(cs%ubtav, cs%vbtav, g%Domain, complete=.false.)
    call pass_vector(uhbtav, vhbtav, 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.
!$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 (present(uhbt_out)) then
      do j=js,je ; do i=is-1,ie ; cs%uhbt_IC(i,j) = uhbt_out(i,j) ; enddo ; enddo
      do j=js-1,je ; do i=is,ie ; cs%vhbt_IC(i,j) = vhbt_out(i,j) ; enddo ; enddo
    elseif (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_Datu_res > 0) call post_data(cs%id_Datu_res, cs%Datu_res(isdb:iedb,jsd:jed), cs%diag)
    if (cs%id_Datv_res > 0) call post_data(cs%id_Datv_res, cs%Datv_res(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 (use_visc_rem) then
      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)
    endif

    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_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)
  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 pass_var_complete(pid_etaav, etaav, g%Domain)
    call pass_vector_complete(pid_uhbtav, uhbtav, vhbtav, g%Domain)
    call pass_vector_complete(pid_ubtav, cs%ubtav, cs%vbtav, g%Domain)
  endif

  if (apply_obcs) call destroy_bt_obc(bt_obc)

end subroutine legacy_btstep

subroutine legacy_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(legacy_barotropic_cs),           pointer       :: CS       !< The control structure returned
                                                          !! by a previous call to barotropic_init.
  real, dimension(SZI_(G),SZJ_(G)),         &
                            intent(in), optional      :: eta      !< The barotropic free surface
                                        !! height anomaly or column mass anomaly, in m or kg m-2.
  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

! Arguments: G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 barotropic_init.
!  (in,opt)  eta - The barotropic free surface height anomaly or
!                  column mass anomaly, in m or kg m-2.
!  (in,opt)  pbce - The baroclinic pressure anomaly in each layer
!                   due to free surface height anomalies, in m2 H-1 s-2.
!  (in,opt)  BT_cont - A structure with elements that describe the effective
!                      open face areas as a function of barotropic flow.
!  (in,opt)  gtot_est - An estimate of the total gravitational acceleration,
!                       in m s-2.
!  (in,opt)  SSH_add - An additional contribution to SSH to provide a margin of
!                      error when calculating the external wave speed, in m.
  ! This subroutine automatically determines an optimal value for dtbt based
  ! on some state of the ocean.
  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, &
      "legacy_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, &
      "legacy_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 legacy_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 ubt 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 u points.
  type(local_bt_cont_u_type), dimension(SZIBW_(MS),SZJW_(MS)), &
                                         intent(in)    :: BTCL_u  !< Structures 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  !< Structures 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
!   This subroutine applies the open boundary conditions on barotropic
! velocities and mass transports, as developed by Mehmet Ilicak.

! Arguments: OBC - an associated pointer to an OBC type.
!  (inout)   ubt - the zonal barotropic velocity, in m s-1.
!  (inout)   vbt - the meridional barotropic velocity, in m s-1.
!  (inout)   uhbt - the zonal barotropic transport, in H m2 s-1.
!  (inout)   vhbt - the meridional barotropic transport, in H m2 s-1.
!  (inout)   ubt_trans - the zonal barotropic velocity used in transport, m s-1.
!  (inout)   vbt_trans - the meridional BT velocity used in transports, m s-1.
!  (in)      eta - The barotropic free surface height anomaly or
!                  column mass anomaly, in m or kg m-2.
!  (in)      ubt_old - The starting value of ubt in a barotropic step, m s-1.
!  (in)      vbt_old - The starting value of ubt in a barotropic step, m s-1.
!  (in)      BT_OBC - A structure with the private barotropic arrays related
!                     to the open boundary conditions, set by set_up_BT_OBC.
!  (in)      G - The ocean's grid structure.
!  (in)      MS - A type that describes the memory sizes of the argument arrays.
!  (in)      halo - The extra halo size to use here.
!  (in)      dtbt - The time step, in s.
!  (in)      bebt - The fractional weighting of the future velocity in
!                   determining the transport.
!  (in)      use_BT_cont - If true, use the BT_cont_types to calculate transports.
!  (in)      Datu - A fixed estimate of the face areas at u points.
!  (in)      Datv - A fixed estimate of the face areas at u points.
!  (in)      BCTL_u - Structures of information used for a dynamic estimate
!  (in)      BCTL_v - of the face areas at u- and v- points.
  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
  integer :: i, j, is, ie, js, je
  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo

  if (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))%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_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))%direction == obc_direction_w) 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))  ! 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) * (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))%direction == obc_direction_n) then
          if ((vbt(i,j-1)+vbt(i+1,j-1)) > 0.0) then
            ubt(i,j) = 2.0*ubt(i,j-1)-ubt(i,j-2)
          else
            ubt(i,j) = bt_obc%ubt_outer(i,j)
          endif
          vel_trans = ubt(i,j)
        elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_s) then
          if ((vbt(i,j)+vbt(i+1,j)) < 0.0) then
            ubt(i,j) = 2.0*ubt(i,j+1)-ubt(i,j+2)
          else
            ubt(i,j) = bt_obc%ubt_outer(i,j)
          endif
          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 (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))%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_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))%direction == obc_direction_s) 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))%direction == obc_direction_e) then
          if ((ubt(i-1,j)+ubt(i-1,j+1)) > 0.0) then
            vbt(i,j) = 2.0*vbt(i-1,j)-vbt(i-2,j)
          else
            vbt(i,j) = bt_obc%vbt_outer(i,j)
          endif
          vel_trans = vbt(i,j)
!!!!!!!!!!!!!!!!!!! CLAMPED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!       cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J)           !
!       vbt(i,J) = (vbt(i-1,J) + CFL*vbt(i,J)) / (1.0 + CFL)  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_w) then
          if ((ubt(i,j)+ubt(i,j+1)) < 0.0) then
            vbt(i,j) = 2.0*vbt(i+1,j)-vbt(i+2,j)
          else
            vbt(i,j) = bt_obc%vbt_outer(i,j)
          endif
          vel_trans = vbt(i,j)
!!!!!!!!!!!!!!!!!! CLAMPED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!       cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J)           !
!       vbt(i,J) = (vbt(i-1,J) + CFL*vbt(i,J)) / (1.0 + CFL)  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        endif
      endif

      if (obc%segnum_v(i,j) /= obc_simple) 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.
!   This subroutine applies the open boundary conditions on the free surface
! height, as coded by Mehmet Ilicak.

! Arguments: OBC - an associated pointer to an OBC type.
!  (inout)   eta - The barotropic free surface height anomaly or
!                  column mass anomaly, in m or kg m-2.
!  (inout)   ubt - the zonal barotropic velocity, in m s-1.
!  (inout)   vbt - the meridional barotropic velocity, in m s-1.
!  (in)      G - The ocean's grid structure.
!  (in)      MS - A type that describes the memory sizes of the argument arrays.
!  (in)      halo - The extra halo size to use here.
!  (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%Flather_u_BCs_exist_globally) .and. 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
      endif
    endif ; enddo ; enddo
  endif

  if ((obc%Flather_v_BCs_exist_globally) .and. 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
      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 v points.
  type(local_bt_cont_u_type), dimension(SZIBW_(MS),SZJW_(MS)), &
                                         intent(in)    :: BTCL_u !< Structures 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 !< Structures of information used
                                         !! for a dynamic estimate of the face areas at v- points.
!   This subroutine sets up the private structure used to apply the open
! boundary conditions, as developed by Mehmet Ilicak.

! Arguments: OBC - an associated pointer to an OBC type.
!  (in)      eta - The barotropic free surface height anomaly or
!                  column mass anomaly, in m or kg m-2.
!  (in)      BT_OBC - A structure with the private barotropic arrays related
!                     to the open boundary conditions, set by set_up_BT_OBC.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      MS - A type that describes the memory sizes of the argument arrays.
!  (in)      halo - The extra halo size to use here.
!  (in)      dtbt - The time step, in s.
!  (in)      use_BT_cont - If true, use the BT_cont_types to calculate transports.
!  (in)      Datu - A fixed estimate of the face areas at u points.
!  (in)      Datv - A fixed estimate of the face areas at u points.
!  (in)      BCTL_u - Structures of information used for a dynamic estimate
!  (in)      BCTL_v - of the face areas at u- and v- points.

  integer :: i, j, k, is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq, n
  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

  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

  if (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 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
      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
        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
          bt_obc%H_u(i,j) = 0.5*((g%bathyT(i,j) + eta(i,j)) + &
                                 (g%bathyT(i+1,j) + eta(i+1,j)))
        else
          bt_obc%H_u(i,j) = 0.5*(eta(i,j) + eta(i+1,j))
        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 (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 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
      if (obc%segnum_v(i,j) == obc_simple) 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
        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
          bt_obc%H_v(i,j) = 0.5*((g%bathyT(i,j) + eta(i,j)) + &
                                 (g%bathyT(i,j+1) + eta(i,j+1)))
        else
          bt_obc%H_v(i,j) = 0.5*(eta(i,j) + eta(i,j+1))
        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

subroutine destroy_bt_obc(BT_OBC)
  type(bt_obc_type), intent(inout) :: BT_OBC

  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)
end subroutine destroy_bt_obc

!> This subroutine 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 legacy_btcalc(h, G, GV, CS, h_u, h_v, may_use_default)
  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(legacy_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-
                                                                 !! and v- points, in m or kg m-2.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                intent(in), optional      :: h_v !< The specified thicknesses at u-
                                                                 !! and v- points, in m or kg m-2.
  logical,                      intent(in), optional      :: may_use_default
!   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.

! Arguments: h - Layer thickness, in m or kg m-2 (H in later comments).
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 barotropic_init.
!  (in,opt)  h_u, h_v - The specified thicknesses at u- and v- points, in m or kg m-2.

! 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 :: htot(szi_(g),szj_(g)) !   The sum of the layer thicknesses, in H.
  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.

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

!    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, &
      "legacy_btcalc: Module MOM_legacy_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, &
        "legacy_btcalc: Inconsistent settings of optional arguments and hvel_scheme.")
  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-1,je+1
    if (present(h_u)) then
      do i=is-2,ie+1 ; hatutot(i) = h_u(i,j,1) ; enddo
      do k=2,nz ; do i=is-2,ie+1
        hatutot(i) = hatutot(i) + h_u(i,j,k)
      enddo ; enddo
      do i=is-2,ie+1 ; ihatutot(i) = 1.0 / (hatutot(i) + h_neglect) ; enddo
      do k=1,nz ; do i=is-2,ie+1
        cs%frhatu(i,j,k) = h_u(i,j,k) * ihatutot(i)
      enddo ; enddo
    else
      if (cs%hvel_scheme == arithmetic) then
        do i=is-2,ie+1
          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-2,ie+1
          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-2,ie+1
          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-2,ie+1
          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-2,ie+1
          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-2,ie+1
          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-2,ie+1 ; ihatutot(i) = 1.0 / (hatutot(i) + h_neglect) ; enddo
      do k=1,nz ; do i=is-2,ie+1
        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-2,je+1
    if (present(h_v)) then
      do i=is-1,ie+1 ; hatvtot(i) = h_v(i,j,1) ; enddo
      do k=2,nz ; do i=is-1,ie+1
        hatvtot(i) = hatvtot(i) + h_v(i,j,k)
      enddo ; enddo
      do i=is-1,ie+1 ; ihatvtot(i) = 1.0 / (hatvtot(i) + h_neglect) ; enddo
      do k=1,nz ; do i=is-1,ie+1
        cs%frhatv(i,j,k) = h_v(i,j,k) * ihatvtot(i)
      enddo ; enddo
    else
      if (cs%hvel_scheme == arithmetic) then
        do i=is-1,ie+1
          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-1,ie+1
          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-1,ie+1
          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-1,ie+1
          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-1,ie+1
          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-1,ie+1
          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-1,ie+1 ; ihatvtot(i) = 1.0 / (hatvtot(i) + h_neglect) ; enddo
      do k=1,nz ; do i=is-1,ie+1
        cs%frhatv(i,j,k) = cs%frhatv(i,j,k) * ihatvtot(i)
      enddo ; enddo
    endif
  enddo

  if (cs%rescale_D_bt) then
    do j=js-2,je+2 ; do i=is-2,ie+2 ; htot(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
      htot(i,j) = htot(i,j) + h(i,j,k)
    enddo ; enddo ; enddo

!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,htot,h_neglect,CS) &
!$OMP                          private(hatutot, h_harm, Rh)
    do j=js-1,je+1
      do i=is-2,ie+1 ; hatutot(i) = 0.0 ; enddo
      do k=1,nz ; do i=is-2,ie+1
        h_harm = 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) + h_harm
      enddo ; enddo
      do i=is-2,ie+1
        rh = hatutot(i) * (htot(i+1,j) + htot(i,j)) / &
             (2.0*(htot(i+1,j) * htot(i,j) + h_neglect**2))
        cs%Datu_res(i,j) = 1.0
        ! This curve satisfies F(1/2) = 1; F'(1/2) = 0; F(x) ~ x for x<<1.
        if (rh < 0.5) cs%Datu_res(i,j) = rh * ((4.0-7.0*rh) / (2.0-3.0*rh)**2)
      enddo
    enddo

!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,htot,h_neglect,CS) &
!$OMP                          private(hatvtot, h_harm, Rh)
    do j=js-2,je+1
      do i=is-1,ie+1 ; hatvtot(i) = 0.0 ; enddo
      do k=1,nz ; do i=is-1,ie+1
        h_harm = 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) + h_harm
      enddo ; enddo
      do i=is-1,ie+1
        rh = hatvtot(i) * (htot(i,j+1) + htot(i,j)) / &
             (2.0*(htot(i,j+1) * htot(i,j) + h_neglect**2))
        cs%Datv_res(i,j) = 1.0
        ! This curve satisfies F(1/2) = 1; F'(1/2) = 0; F(x) ~ x for x<<1.
        if (rh < 0.5) cs%Datv_res(i,j) = rh * ((4.0-7.0*rh) / (2.0-3.0*rh)**2)
      enddo
    enddo
  endif

  if (cs%debug) then
    call uvchksum("btcalc frhat[uv]", cs%frhatu, cs%frhatv, g%HI, haloshift=1)
    call hchksum(h, "btcalc h", g%HI, haloshift=1, scale=gv%H_to_m)
  endif

end subroutine legacy_btcalc

function find_uhbt(u, BTC) result(uhbt)
  real, intent(in) :: u    !< The zonal velocity, in m s-1
  type(local_bt_cont_u_type), intent(in) :: BTC
  real :: uhbt ! The result
  ! This function evaluates the zonal transport function.

  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

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
  !   This function inverts the transport function to determine the barotopic
  ! velocity that is consistent with a given transport.
  ! Arguments: uhbt - The barotropic zonal transport that should be inverted
  !                   for, in units of H m2 s-1.
  !  (in)      BTC - A structure containing various fields that allow the
  !                  barotropic transports to be calculated consistently with
  !                  the layers' continuity equations.
  !  (in,opt)  FA_rat_EE - The current value of the far-east face area divided
  !                        by its value when uhbt was originally calculated, ND.
  !  (in,opt)  FA_rat_WW - The current value of the far-west face area divided
  !                        by its value when uhbt was originally calculated, ND.
  !  (in, opt) guess - A guess at what ubt will be.  The result is not allowed
  !                    to be dramatically larger than guess.
  ! result: ubt - The velocity that gives uhbt transport, in m s-1.
  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

function find_vhbt(v, BTC) result(vhbt)
  real, intent(in) :: v    !< The meridional velocity, in m s-1
  type(local_bt_cont_v_type), intent(in) :: BTC
  real :: vhbt ! The result
  ! This function evaluates the meridional transport function.

  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

function vhbt_to_vbt(vhbt, BTC, guess) result(vbt)
  real,                       intent(in) :: vhbt
  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.
  !   This function inverts the transport function to determine the barotopic
  ! velocity that is consistent with a given transport.
  ! Arguments: vhbt_in - The barotropic meridional transport that should be
  !                      inverted for, in units of H m2 s-1.
  !  (in)      BTC - A structure containing various fields that allow the
  !                  barotropic transports to be calculated consistently with
  !                  the layers' continuity equations.
  !  (in,opt)  FA_rat_NN - The current value of the far-north face area divided
  !                        by its value when vhbt was originally calculated, ND.
  !  (in,opt)  FA_rat_SS - The current value of the far-south face area divided
  !                        by its value when vhbt was originally calculated, ND.
  !  (in, opt) guess - A guess at what vbt will be.  The result is not allowed
  !                    to be dramatically larger than guess.
  ! result: vbt - The velocity that gives vhbt transport, in m s-1.
  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


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(inout) :: 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.
!   This subroutine sets up reordered versions of the BT_cont type in the
! local_BT_cont types, which have wide halos properly filled in.
! Arguments: BT_cont - The BT_cont_type input to the barotropic solver.
!  (out)     BTCL_u - A structure with the u information from BT_cont.
!  (out)     BTCL_v - A structure with the v information from BT_cont.
!  (in)      G - The ocean's grid structure.
!  (in)      MS - A type that describes the memory sizes of the argument arrays.
!  (in)      BT_Domain - The domain to use for updating the halos of wide arrays.
!  (in)      halo - The extra halo size to use here.

  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.
  u_polarity(:,:) = 1.0
  ubt_ee(:,:) = 0.0 ; ubt_ww(:,:) = 0.0
  fa_u_ee(:,:) = 0.0 ; fa_u_e0(:,:) = 0.0 ; fa_u_w0(:,:) = 0.0 ; fa_u_ww(:,:) = 0.0
  do i=is-1,ie ; do j=js,je
    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

  v_polarity(:,:) = 1.0
  vbt_nn(:,:) = 0.0 ; vbt_ss(:,:) = 0.0
  fa_v_nn(:,:) = 0.0 ; fa_v_n0(:,:) = 0.0 ; fa_v_s0(:,:) = 0.0 ; fa_v_ss(:,:) = 0.0
  do i=is,ie ; do j=js-1,je
    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

  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)
  ! Do halo updates on BT_cont.
  call pass_vector(u_polarity, v_polarity, bt_domain, complete=.false.)
  call pass_vector(ubt_ee, vbt_nn, bt_domain, complete=.false.)
  call pass_vector(ubt_ww, vbt_ss, bt_domain, complete=.true.)

  call pass_vector(fa_u_ee, fa_v_nn, bt_domain, to_all+scalar_pair, complete=.false.)
  call pass_vector(fa_u_e0, fa_v_n0, bt_domain, to_all+scalar_pair, complete=.false.)
  call pass_vector(fa_u_w0, fa_v_s0, bt_domain, to_all+scalar_pair, complete=.false.)
  call pass_vector(fa_u_ww, fa_v_ss, bt_domain, to_all+scalar_pair, complete=.true.)
  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)

  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

  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

end subroutine set_local_bt_cont_types

subroutine bt_cont_to_face_areas(BT_cont, Datu, Datv, G, MS, halo, maximize)
  type(bt_cont_type),                         intent(inout) :: BT_cont
  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 !< 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.
  integer,                          optional, intent(in)    :: halo !< The halo size to use,
                                                                    !! default = 1.
  logical,                          optional, intent(in)    :: maximize
  !   This subroutine uses the BTCL types to find typical or maximum face
  ! areas, which can then be used for finding wave speeds, etc.
  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

subroutine find_face_areas(Datu, Datv, G, GV, CS, MS, rescale_faces, eta, halo, add_max)
  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 !< 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(legacy_barotropic_cs),        pointer     :: CS   !< The control structure returned by a
                                                         !! previous call to barotropic_init.
  logical,                 optional, intent(in)  :: rescale_faces !< If true, rescale the face areas
                                                                  !! by Datu_res, etc.
  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 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.
! Arguments: Datu - The open zonal face area, in H m (m2 or kg m-1).
!  (out)     Datv - The open meridional face area, in H m (m2 or kg m-1).
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 barotropic_init.
!  (in)      MS - A type that describes the memory sizes of the argument arrays.
!  (in, opt) rescale_faces - If true, rescale the face areas by Datu_res, etc.
!  (in, opt) eta - The barotropic free surface height anomaly or
!                  column mass anomaly, in m or kg m-2.
!  (in, opt) halo - The halo size to use, default = 1.
!  (in, opt) add_max - A value to add to the maximum depth (used to overestimate
!                      the external wave speed) in m.


!   This subroutine determines the open face areas of cells for calculating
! the barotropic transport.
  real :: H1, H2      ! Temporary total thicknesses, in m or kg m-2.
  logical :: rescale
  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)
  rescale = .false. ; if (present(rescale_faces)) rescale = rescale_faces

!$OMP parallel default(none) shared(is,ie,js,je,hs,eta,GV,CS,Datu,Datv,add_max,rescale) &
!$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) + eta(i,j) ; h2 = cs%bathyT(i+1,j) + 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) + eta(i,j) ; h2 = cs%bathyT(i,j+1) + 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) = 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) = 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

  if (rescale) then
!$OMP do
    do j=js-hs,je+hs ; do i=is-1-hs,ie+hs
      datu(i,j) = datu(i,j) * cs%Datu_res(i,j)
    enddo ; enddo
!$OMP do
    do j=js-1-hs,je+hs ; do i=is-hs,ie+hs
      datv(i,j) = datv(i,j) * cs%Datv_res(i,j)
    enddo ; enddo
  endif
!$OMP end parallel

end subroutine find_face_areas

subroutine legacy_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(legacy_barotropic_cs),  pointer    :: CS       !< The control structure returned by a
                                                      !! previous call to barotropic_init.

!   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.

! Arguments: h - Layer thickness, in m or kg m-2 (H).
!  (in)      eta - The free surface height that is to be corrected, in m.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (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.
!  (in)      dt_therm - The thermodynamic time step, in s.
!  (in)      dt_since_therm - The elapsed time since mass forcing was applied, s.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 barotropic_init.
  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) ; 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 legacy_bt_mass_source

subroutine legacy_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(legacy_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.
!   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.

! Arguments: u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m or kg m-2.
!  (in)      eta - Free surface height or column mass anomaly, in m or kg m-2.
!  (in)      Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (in/out)  CS - A pointer to the control structure for this module that is
!                 set in register_barotropic_restarts.
!  (in)      restart_CS - A pointer to the restart control structure.
!  (in,opt)  BT_cont - A structure with elements that describe the effective
!                      open face areas as a function of barotropic flow.
!  (in)      tides_CSp - a pointer to the control structure of the tide module.
! This include declares and sets the variable "version".
#include "version_variable.h"
  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
  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

  ! ### USE SOMETHING OTHER THAN MAXVEL FOR THIS...
  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 0.1*MAXVEL to be accommodated.",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_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, "RESCALE_BT_FACE_AREAS", cs%rescale_D_bt, &
                 "If true, the face areas used by the barotropic solver \n"//&
                 "are rescaled to approximately reflect the open face \n"//&
                 "areas of the interior layers.  This probably requires \n"//&
                 "FLUX_BT_COUPLING to work, and should not be used with \n"//&
                 "USE_BT_CONT_TYPE.", default=.false.)
  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, "APPLY_BT_DRAG", apply_bt_drag, &
                 "If defined, bottom drag is applied within the \n"//&
                 "barotropic solver.", default=.true.)
  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=.true.)

  call get_param(param_file, mdl, "CLIP_BT_VELOCITY", cs%clip_velocity, &
                 "If true, limit any velocity components that exceed \n"//&
                 "MAXVEL.  This should only be used as a desperate \n"//&
                 "debugging measure.", default=.false.)
  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.1)

  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)


  if (apply_bt_drag) then ; cs%drag_amp = 1.0 ; else ; cs%drag_amp = 0.0 ; endif

  ! 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%Datu_res(isdw-1:iedw,jsdw:jedw))
  alloc_(cs%Datv_res(isdw:iedw,jsdw-1:jedw))
  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%Datu_res(:,:) = 1.0 ; cs%Datv_res(:,:) = 1.0

  cs%ua_polarity(:,:) = 1.0 ; cs%va_polarity(:,:) = 1.0
  call pass_vector(cs%ua_polarity, cs%va_polarity, cs%BT_domain, to_all, agrid)

  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 !### Should this be 0 instead?
  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 pass_var(cs%IareaT, cs%BT_domain, to_all)
  call pass_var(cs%bathyT, cs%BT_domain, to_all)
  call pass_vector(cs%IdxCu, cs%IdyCv, cs%BT_domain, to_all+scalar_pair)
  call pass_vector(cs%dy_Cu, cs%dx_Cv, cs%BT_domain, to_all+scalar_pair)

  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
      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)))
    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 pass_var(cs%q_D, cs%BT_Domain, to_all, position=corner)
    call pass_vector(cs%D_u_Cor, cs%D_v_Cor, cs%BT_Domain, to_all+scalar_pair)
  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 legacy_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')
  if (cs%rescale_D_bt) then
    cs%id_Datu_res = register_diag_field('ocean_model', 'Datu_res', diag%axesCu1, time, &
      'Rescaling for zonal face area in barotropic continuity', 'Nondimensional')
    cs%id_Datv_res = register_diag_field('ocean_model', 'Datv_res', diag%axesCv1, time, &
      'Rescaling for meridional face area in barotropic continuity', 'Nondimensional')
  endif
  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 (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 legacy_btcalc(h, g, gv, cs, may_use_default=.true.)
    do j=js-1,je+1 ; do i=is-1,ie+1
      cs%ubtav(i,j) = 0.0 ; cs%vbtav(i,j) = 0.0
    enddo ; enddo
    do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
      cs%ubtav(i,j) = cs%ubtav(i,j) + cs%frhatu(i,j,k) * u(i,j,k)
      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
      cs%IDatu(i,j) = g%mask2dCu(i,j) * 2.0 / (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 / (g%bathyT(i,j+1) + g%bathyT(i,j))
    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
    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 pass_vector(cs%ubt_IC, cs%vbt_IC, g%Domain, complete=.false.)
  call pass_vector(cs%uhbt_IC, cs%vhbt_IC, g%Domain, complete=.false.)
  call pass_vector(cs%ubtav, cs%vbtav, 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 legacy_barotropic_init

subroutine legacy_barotropic_end(CS)
  type(legacy_barotropic_cs), pointer :: CS
  dealloc_(cs%frhatu)   ; dealloc_(cs%frhatv)
  dealloc_(cs%IDatu)    ; dealloc_(cs%IDatv)
  dealloc_(cs%Datu_res) ; dealloc_(cs%Datv_res)
  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

  deallocate(cs)
end subroutine legacy_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_legacy_barotropic_restarts(HI, GV, param_file, CS, restart_CS)
  type(hor_index_type),    intent(in) :: HI         !< A horizontal index type structure.
  type(verticalgrid_type), intent(in) :: GV         !< The ocean's vertical grid structure.
  type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters.
  type(legacy_barotropic_cs), pointer :: CS         !< A pointer that is set to point to the control
                                                    !! structure for this module.
  type(mom_restart_cs),    pointer    :: restart_CS !< A pointer to the restart control structure.
! This subroutine is used to register any fields from MOM_barotropic.F90
! that should be written to or read from the restart file.
! Arguments: HI - A horizontal index type structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
!  (in)      restart_CS - A pointer to the restart control structure.
  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_legacy_barotropic_restarts

end module mom_legacy_barotropic