MOM_continuity_PPM.F90

!> Solve the layer continuity equation using the PPM method for layer fluxes.
module mom_continuity_ppm

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

use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
use mom_diag_mediator, only : time_type, diag_ctrl
use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_grid, only : ocean_grid_type
use mom_open_boundary, only : ocean_obc_type, obc_segment_type, obc_none
use mom_open_boundary, only : obc_direction_e, obc_direction_w, obc_direction_n, obc_direction_s
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : bt_cont_type
use mom_verticalgrid, only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

public continuity_ppm, continuity_ppm_init, continuity_ppm_end, continuity_ppm_stencil

!>@{ CPU time clock IDs
integer :: id_clock_update, id_clock_correct
!>@}

!> Control structure for mom_continuity_ppm
type, public :: continuity_ppm_cs ; private
  type(diag_ctrl), pointer :: diag !< Diagnostics control structure.
  logical :: upwind_1st      !< If true, use a first-order upwind scheme.
  logical :: monotonic       !< If true, use the Colella & Woodward monotonic
                             !! limiter; otherwise use a simple positive
                             !! definite limiter.
  logical :: simple_2nd      !< If true, use a simple second order (arithmetic
                             !! mean) interpolation of the edge values instead
                             !! of the higher order interpolation.
  real :: tol_eta            !< The tolerance for free-surface height
                             !! discrepancies between the barotropic solution and
                             !! the sum of the layer thicknesses [H ~> m or kg m-2].
  real :: tol_vel            !< The tolerance for barotropic velocity
                             !! discrepancies between the barotropic solution and
                             !! the sum of the layer thicknesses [L T-1 ~> m s-1].
  real :: tol_eta_aux        !< The tolerance for free-surface height
                             !! discrepancies between the barotropic solution and
                             !! the sum of the layer thicknesses when calculating
                             !! the auxiliary corrected velocities [H ~> m or kg m-2].
  real :: cfl_limit_adjust   !< The maximum CFL of the adjusted velocities [nondim]
  logical :: aggress_adjust  !< If true, allow the adjusted velocities to have a
                             !! relative CFL change up to 0.5.  False by default.
  logical :: vol_cfl         !< If true, use the ratio of the open face lengths
                             !! to the tracer cell areas when estimating CFL
                             !! numbers.  Without aggress_adjust, the default is
                             !! false; it is always true with.
  logical :: better_iter     !< If true, stop corrective iterations using a
                             !! velocity-based criterion and only stop if the
                             !! iteration is better than all predecessors.
  logical :: use_visc_rem_max !< If true, use more appropriate limiting bounds
                             !! for corrections in strongly viscous columns.
  logical :: marginal_faces  !< If true, use the marginal face areas from the
                             !! continuity solver for use as the weights in the
                             !! barotropic solver.  Otherwise use the transport
                             !! averaged areas.
end type continuity_ppm_cs

!> A container for loop bounds
type :: loop_bounds_type ; private
  !>@{ Loop bounds
  integer :: ish, ieh, jsh, jeh
  !>@}
end type loop_bounds_type

contains

!> Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme,
!! based on Lin (1994).
subroutine continuity_ppm(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, &
                          visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont)
  type(ocean_grid_type),   intent(inout) :: g   !< The ocean's grid structure.
  type(continuity_ppm_cs), pointer       :: cs  !< Module's control structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u   !< Zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v   !< Meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in)    :: hin !< Initial layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(inout) :: h   !< Final layer thickness [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: uh  !< Zonal volume flux, u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(out)   :: vh  !< Meridional volume flux, v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].
  real,                    intent(in)    :: dt  !< Time increment [T ~> s].
  type(verticalgrid_type), intent(in)    :: gv  !< Vertical grid structure.
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G)), &
                 optional, intent(in)    :: uhbt !< The summed volume flux through zonal faces
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G)), &
                 optional, intent(in)    :: vhbt !< The summed volume flux through meridional faces
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  type(ocean_obc_type),  &
                 optional, pointer       :: obc !< Open boundaries control structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(in)    :: visc_rem_u
                             !< The fraction of zonal momentum originally
                             !! in a layer that remains after a time-step of viscosity, and the
                             !! fraction of a time-step's worth of a barotropic acceleration that
                             !! a layer experiences after viscosity is applied.
                             !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(in)    :: visc_rem_v
                             !< The fraction of meridional momentum originally
                             !! in a layer that remains after a time-step of viscosity, and the
                             !! fraction of a time-step's worth of a barotropic acceleration that
                             !! a layer experiences after viscosity is applied.
                             !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: u_cor
                             !< The zonal velocities that give uhbt as the depth-integrated transport [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(out)   :: v_cor
                             !< The meridional velocities that give vhbt as the depth-integrated
                             !! transport [L T-1 ~> m s-1].
  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.

  ! Local variables
  real :: h_min  ! The minimum layer thickness [H ~> m or kg m-2].  h_min could be 0.
  type(loop_bounds_type) :: lb
  integer :: is, ie, js, je, nz, stencil
  integer :: i, j, k

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

  h_min = gv%Angstrom_H

  if (.not.associated(cs)) call mom_error(fatal, &
         "MOM_continuity_PPM: Module must be initialized before it is used.")
  x_first = (mod(g%first_direction,2) == 0)

  if (present(visc_rem_u) .neqv. present(visc_rem_v)) call mom_error(fatal, &
      "MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither"// &
      " one must be present in call to continuity_PPM.")

  stencil = 3 ; if (cs%simple_2nd) stencil = 2 ; if (cs%upwind_1st) stencil = 1

  if (x_first) then
  !    First, advect zonally.
    lb%ish = g%isc ; lb%ieh = g%iec
    lb%jsh = g%jsc-stencil ; lb%jeh = g%jec+stencil
    call zonal_mass_flux(u, hin, uh, dt, g, gv, us, cs, lb, uhbt, obc, visc_rem_u, u_cor, bt_cont)

    call cpu_clock_begin(id_clock_update)
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
      h(i,j,k) = hin(i,j,k) - dt * g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k))
  !   Uncomment this line to prevent underflow.
  !   if (h(i,j,k) < h_min) h(i,j,k) = h_min
    enddo ; enddo ; enddo
    call cpu_clock_end(id_clock_update)

    lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec

    !    Now advect meridionally, using the updated thicknesses to determine
    !  the fluxes.
    call meridional_mass_flux(v, h, vh, dt, g, gv, us, cs, lb, vhbt, obc, visc_rem_v, v_cor, bt_cont)

    call cpu_clock_begin(id_clock_update)
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
      h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k))
  !   This line prevents underflow.
      if (h(i,j,k) < h_min) h(i,j,k) = h_min
    enddo ; enddo ; enddo
    call cpu_clock_end(id_clock_update)

  else  ! .not. x_first
  !    First, advect meridionally, so set the loop bounds accordingly.
    lb%ish = g%isc-stencil ; lb%ieh = g%iec+stencil
    lb%jsh = g%jsc ; lb%jeh = g%jec

    call meridional_mass_flux(v, hin, vh, dt, g, gv, us, cs, lb, vhbt, obc, visc_rem_v, v_cor, bt_cont)

    call cpu_clock_begin(id_clock_update)
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
      h(i,j,k) = hin(i,j,k) - dt * g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k))
    enddo ; enddo ; enddo
    call cpu_clock_end(id_clock_update)

  !    Now advect zonally, using the updated thicknesses to determine
  !  the fluxes.
    lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec
    call zonal_mass_flux(u, h, uh, dt, g, gv, us, cs, lb, uhbt, obc, visc_rem_u, u_cor, bt_cont)

    call cpu_clock_begin(id_clock_update)
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
      h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k))
      ! This line prevents underflow.
      if (h(i,j,k) < h_min) h(i,j,k) = h_min
    enddo ; enddo ; enddo
    call cpu_clock_end(id_clock_update)

  endif

end subroutine continuity_ppm

!> Calculates the mass or volume fluxes through the zonal faces, and other related quantities.
subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, &
                           visc_rem_u, u_cor, BT_cont)
  type(ocean_grid_type),   intent(inout) :: G    !< Ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: GV   !< Ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u    !< Zonal velocity [L T-1 ~> m s-1].
  real,  dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: uh   !< Volume flux through zonal faces = u*h*dy
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real,                    intent(in)    :: dt   !< Time increment [T ~> s].
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_ppm_cs), pointer       :: CS   !< This module's control structure.
  type(loop_bounds_type),  intent(in)    :: LB   !< Loop bounds structure.
  type(ocean_obc_type), &
                 optional, pointer       :: OBC  !< Open boundaries control structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(in)    :: visc_rem_u
                     !< The fraction of zonal momentum originally in a layer that remains after a
                     !! time-step of viscosity, and the fraction of a time-step's worth of a barotropic
                     !! acceleration that a layer experiences after viscosity is applied.
                     !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G),SZJ_(G)), &
                 optional, intent(in)    :: uhbt !< The summed volume flux through zonal faces
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: u_cor
                     !< The zonal velocitiess (u with a barotropic correction)
                     !! that give uhbt as the depth-integrated transport, m s-1.
  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.

  ! Local variables
  real, dimension(SZIB_(G),SZK_(G)) :: duhdu ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_L, h_R ! Left and right face thicknesses [H ~> m or kg m-2].
  real, dimension(SZIB_(G)) :: &
    du, &      ! Corrective barotropic change in the velocity [L T-1 ~> m s-1].
    du_min_CFL, & ! Min/max limits on du correction
    du_max_CFL, & ! to avoid CFL violations
    duhdu_tot_0, & ! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
    uh_tot_0, & ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
    visc_rem_max  ! The column maximum of visc_rem.
  logical, dimension(SZIB_(G)) :: do_I
  real, dimension(SZIB_(G),SZK_(G)) :: &
    visc_rem      ! A 2-D copy of visc_rem_u or an array of 1's.
  real, dimension(SZIB_(G)) :: FAuI  ! A list of sums of zonal face areas [H L ~> m2 or kg m-1].
  real :: FA_u    ! A sum of zonal face areas [H m ~> m2 or kg m-1].
  real :: I_vrm   ! 1.0 / visc_rem_max, nondim.
  real :: CFL_dt  ! The maximum CFL ratio of the adjusted velocities divided by
                  ! the time step [T-1 ~> s-1].
  real :: I_dt    ! 1.0 / dt [T-1 ~> s-1].
  real :: du_lim  ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1].
  real :: dx_E, dx_W ! Effective x-grid spacings to the east and west [L ~> m].
  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
  integer :: l_seg
  logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC
  logical :: local_Flather_OBC, local_open_BC, is_simple
  type(obc_segment_type), pointer :: segment => null()

  use_visc_rem = present(visc_rem_u)
  local_specified_bc = .false. ; set_bt_cont = .false. ; local_flather_obc = .false.
  local_open_bc = .false.
  if (present(bt_cont)) set_bt_cont = (associated(bt_cont))
  if (present(obc)) then ; if (associated(obc)) then
    local_specified_bc = obc%specified_u_BCs_exist_globally
    local_flather_obc = obc%Flather_u_BCs_exist_globally
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke

  cfl_dt = cs%CFL_limit_adjust / dt
  i_dt = 1.0 / dt
  if (cs%aggress_adjust) cfl_dt = i_dt

  call cpu_clock_begin(id_clock_update)
!$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,CS,h_L,h_in,h_R,G,GV,LB,visc_rem,OBC)
  do k=1,nz
    ! This sets h_L and h_R.
    if (cs%upwind_1st) then
      do j=jsh,jeh ; do i=ish-1,ieh+1
        h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
      enddo ; enddo
    else
      call ppm_reconstruction_x(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
                                2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
    endif
    do i=ish-1,ieh ; visc_rem(i,k) = 1.0 ; enddo
  enddo
  call cpu_clock_end(id_clock_update)

  call cpu_clock_begin(id_clock_correct)
!$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,u,h_in,h_L,h_R,use_visc_rem,visc_rem_u,  &
!$OMP                                  uh,dt,US,G,GV,CS,local_specified_BC,OBC,uhbt,set_BT_cont,    &
!$OMP                                  CFL_dt,I_dt,u_cor,BT_cont, local_Flather_OBC) &
!$OMP                          private(do_I,duhdu,du,du_max_CFL,du_min_CFL,uh_tot_0,duhdu_tot_0, &
!$OMP                                  is_simple,FAuI,visc_rem_max,I_vrm,du_lim,dx_E,dx_W, &
!$OMP                                  any_simple_OBC,l_seg) &
!$OMP      firstprivate(visc_rem)
  do j=jsh,jeh
    do i=ish-1,ieh ; do_i(i) = .true. ; visc_rem_max(i) = 0.0 ; enddo
    ! Set uh and duhdu.
    do k=1,nz
      if (use_visc_rem) then ; do i=ish-1,ieh
        visc_rem(i,k) = visc_rem_u(i,j,k)
        visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))
      enddo ; endif
      call zonal_flux_layer(u(:,j,k), h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), &
                            uh(:,j,k), duhdu(:,k), visc_rem(:,k), &
                            dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
      if (local_specified_bc) then
        do i=ish-1,ieh
          l_seg = obc%segnum_u(i,j)

          if (l_seg /= obc_none) then
            if (obc%segment(l_seg)%specified) &
              uh(i,j,k) = obc%segment(l_seg)%normal_trans(i,j,k)
          endif
        enddo
      endif
    enddo

    if ((.not.use_visc_rem).or.(.not.cs%use_visc_rem_max)) then ; do i=ish-1,ieh
      visc_rem_max(i) = 1.0
    enddo ; endif

    if (present(uhbt) .or. set_bt_cont) then
      !   Set limits on du that will keep the CFL number between -1 and 1.
      ! This should be adequate to keep the root bracketed in all cases.
      do i=ish-1,ieh
        i_vrm = 0.0
        if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)
        if (cs%vol_CFL) then
          dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
          dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
        else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif
        du_max_cfl(i) = 2.0* (cfl_dt * dx_w) * i_vrm
        du_min_cfl(i) = -2.0 * (cfl_dt * dx_e) * i_vrm
        uh_tot_0(i) = 0.0 ; duhdu_tot_0(i) = 0.0
      enddo
      do k=1,nz ; do i=ish-1,ieh
        duhdu_tot_0(i) = duhdu_tot_0(i) + duhdu(i,k)
        uh_tot_0(i) = uh_tot_0(i) + uh(i,j,k)
      enddo ; enddo
      if (use_visc_rem) then
        if (cs%aggress_adjust) then
          do k=1,nz ; do i=ish-1,ieh
            if (cs%vol_CFL) then
              dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
              dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
            else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

            du_lim = 0.499*((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k)))
            if (du_max_cfl(i) * visc_rem(i,k) > du_lim) &
              du_max_cfl(i) = du_lim / visc_rem(i,k)

            du_lim = 0.499*((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k)))
            if (du_min_cfl(i) * visc_rem(i,k) < du_lim) &
              du_min_cfl(i) = du_lim / visc_rem(i,k)
          enddo ; enddo
        else
          do k=1,nz ; do i=ish-1,ieh
            if (cs%vol_CFL) then
              dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
              dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
            else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

            if (du_max_cfl(i) * visc_rem(i,k) > dx_w*cfl_dt - u(i,j,k)) &
              du_max_cfl(i) = (dx_w*cfl_dt - u(i,j,k)) / visc_rem(i,k)
            if (du_min_cfl(i) * visc_rem(i,k) < -dx_e*cfl_dt - u(i,j,k)) &
              du_min_cfl(i) = -(dx_e*cfl_dt + u(i,j,k)) / visc_rem(i,k)
          enddo ; enddo
        endif
      else
        if (cs%aggress_adjust) then
          do k=1,nz ; do i=ish-1,ieh
            if (cs%vol_CFL) then
              dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
              dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
            else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

            du_max_cfl(i) = min(du_max_cfl(i), 0.499 * &
                        ((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k))) )
            du_min_cfl(i) = max(du_min_cfl(i), 0.499 * &
                        ((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k))) )
          enddo ; enddo
        else
          do k=1,nz ; do i=ish-1,ieh
            if (cs%vol_CFL) then
              dx_w = ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
              dx_e = ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
            else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ; endif

            du_max_cfl(i) = min(du_max_cfl(i), dx_w*cfl_dt - u(i,j,k))
            du_min_cfl(i) = max(du_min_cfl(i), -(dx_e*cfl_dt + u(i,j,k)))
          enddo ; enddo
        endif
      endif
      do i=ish-1,ieh
        du_max_cfl(i) = max(du_max_cfl(i),0.0)
        du_min_cfl(i) = min(du_min_cfl(i),0.0)
      enddo

      any_simple_obc = .false.
      if (present(uhbt) .or. set_bt_cont) then
        if (local_specified_bc .or. local_flather_obc) then ; do i=ish-1,ieh
          l_seg = obc%segnum_u(i,j)

          ! Avoid reconciling barotropic/baroclinic transports if transport is specified
          is_simple = .false.
          if (l_seg /= obc_none) &
            is_simple = obc%segment(l_seg)%specified
          do_i(i) = .not. (l_seg /= obc_none .and. is_simple)
          any_simple_obc = any_simple_obc .or. is_simple
        enddo ; else ; do i=ish-1,ieh
          do_i(i) = .true.
        enddo ; endif
      endif

      if (present(uhbt)) then
        call zonal_flux_adjust(u, h_in, h_l, h_r, uhbt(:,j), uh_tot_0, duhdu_tot_0, du, &
                               du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
                               j, ish, ieh, do_i, .true., uh, obc=obc)

        if (present(u_cor)) then ; do k=1,nz
          do i=ish-1,ieh ; u_cor(i,j,k) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
          if (local_specified_bc) then ; do i=ish-1,ieh
            l_seg = obc%segnum_u(i,j)

            if (l_seg /= obc_none) then
              if (obc%segment(l_seg)%specified) &
                u_cor(i,j,k) = obc%segment(l_seg)%normal_vel(i,j,k)
            endif
          enddo ; endif
        enddo ; endif ! u-corrected

      endif

      if (set_bt_cont) then
        call set_zonal_bt_cont(u, h_in, h_l, h_r, bt_cont, uh_tot_0, duhdu_tot_0,&
                               du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
                               visc_rem_max, j, ish, ieh, do_i)
        if (any_simple_obc) then
          do i=ish-1,ieh
            l_seg = obc%segnum_u(i,j)

            do_i(i) = .false.
            if (l_seg /= obc_none) &
              do_i(i) = obc%segment(l_seg)%specified

            if (do_i(i)) faui(i) = gv%H_subroundoff*g%dy_Cu(i,j)
          enddo
          ! NOTE: do_I(I) should prevent access to segment OBC_NONE
          do k=1,nz ; do i=ish-1,ieh ; if (do_i(i)) then
            if ((abs(obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)) > 0.0) .and. &
                (obc%segment(obc%segnum_u(i,j))%specified)) &
              faui(i) = faui(i) + obc%segment(obc%segnum_u(i,j))%normal_trans(i,j,k) / &
                                  obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
          endif ; enddo ; enddo
          do i=ish-1,ieh ; if (do_i(i)) then
            bt_cont%FA_u_W0(i,j) = faui(i) ; bt_cont%FA_u_E0(i,j) = faui(i)
            bt_cont%FA_u_WW(i,j) = faui(i) ; bt_cont%FA_u_EE(i,j) = faui(i)
            bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
          endif ; enddo
        endif
      endif ! set_BT_cont

    endif ! present(uhbt) or set_BT_cont

  enddo ! j-loop

  if (local_open_bc .and. set_bt_cont) then
    do n = 1, obc%number_of_segments
      if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then
        i = obc%segment(n)%HI%IsdB
        if (obc%segment(n)%direction == obc_direction_e) then
          do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
            fa_u = 0.0
            do k=1,nz ; fa_u = fa_u + h_in(i,j,k)*g%dy_Cu(i,j) ; enddo
            bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
            bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
            bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
          enddo
        else
          do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
            fa_u = 0.0
            do k=1,nz ; fa_u = fa_u + h_in(i+1,j,k)*g%dy_Cu(i,j) ; enddo
            bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
            bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
            bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
          enddo
        endif
      endif
    enddo
  endif
  call cpu_clock_end(id_clock_correct)

  if  (set_bt_cont) then ; if (allocated(bt_cont%h_u)) then
    if (present(u_cor)) then
      call zonal_face_thickness(u_cor, h_in, h_l, h_r, bt_cont%h_u, dt, g, us, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
    else
      call zonal_face_thickness(u, h_in, h_l, h_r, bt_cont%h_u, dt, g, us, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
    endif
  endif ; endif

end subroutine zonal_mass_flux

!> Evaluates the zonal mass or volume fluxes in a layer.
subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, &
                            ish, ieh, do_I, vol_CFL, OBC)
  type(ocean_grid_type),        intent(inout) :: G        !< Ocean's grid structure.
  real, dimension(SZIB_(G)),    intent(in)    :: u        !< Zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZIB_(G)),    intent(in)    :: visc_rem !< 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.
                        !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G)),     intent(in)    :: h        !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G)),     intent(in)    :: h_L      !< Left thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G)),     intent(in)    :: h_R      !< Right thickness [H ~> m or kg m-2].
  real, dimension(SZIB_(G)),    intent(inout) :: uh       !< Zonal mass or volume
                                                          !! transport [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G)),    intent(inout) :: duhdu    !< Partial derivative of uh
                                                          !! with u [H L ~> m2 or kg m-1].
  real,                         intent(in)    :: dt  !< Time increment [T ~> s].
  type(unit_scale_type),        intent(in)    :: US       !< A dimensional unit scaling type
  integer,                      intent(in)    :: j        !< Spatial index.
  integer,                      intent(in)    :: ish      !< Start of index range.
  integer,                      intent(in)    :: ieh      !< End of index range.
  logical, dimension(SZIB_(G)), intent(in)    :: do_I     !< Which i values to work on.
  logical,                      intent(in)    :: vol_CFL  !< If true, rescale the
          !! ratio of face areas to the cell areas when estimating the CFL number.
  type(ocean_obc_type), optional, pointer     :: OBC !< Open boundaries control structure.
  ! Local variables
  real :: CFL  ! The CFL number based on the local velocity and grid spacing [nondim]
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h_in.
  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
  integer :: i
  integer :: l_seg
  logical :: local_open_BC

  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif

  do i=ish-1,ieh ; if (do_i(i)) then
    ! Set new values of uh and duhdu.
    if (u(i) > 0.0) then
      if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
      else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif
      curv_3 = h_l(i) + h_r(i) - 2.0*h(i)
      uh(i) = g%dy_Cu(i,j) * u(i) * &
          (h_r(i) + cfl * (0.5*(h_l(i) - h_r(i)) + curv_3*(cfl - 1.5)))
      h_marg = h_r(i) + cfl * ((h_l(i) - h_r(i)) + 3.0*curv_3*(cfl - 1.0))
    elseif (u(i) < 0.0) then
      if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
      else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif
      curv_3 = h_l(i+1) + h_r(i+1) - 2.0*h(i+1)
      uh(i) = g%dy_Cu(i,j) * u(i) * &
          (h_l(i+1) + cfl * (0.5*(h_r(i+1)-h_l(i+1)) + curv_3*(cfl - 1.5)))
      h_marg = h_l(i+1) + cfl * ((h_r(i+1)-h_l(i+1)) + 3.0*curv_3*(cfl - 1.0))
    else
      uh(i) = 0.0
      h_marg = 0.5 * (h_l(i+1) + h_r(i))
    endif
    duhdu(i) = g%dy_Cu(i,j) * h_marg * visc_rem(i)
  endif ; enddo

  if (local_open_bc) then
    do i=ish-1,ieh ; if (do_i(i)) then
      l_seg = obc%segnum_u(i,j)

      if (l_seg /= obc_none) then
        if (obc%segment(l_seg)%open) then
          if (obc%segment(l_seg)%direction == obc_direction_e) then
            uh(i) = g%dy_Cu(i,j) * u(i) * h(i)
            duhdu(i) = g%dy_Cu(i,j) * h(i) * visc_rem(i)
          else
            uh(i) = g%dy_Cu(i,j) * u(i) * h(i+1)
            duhdu(i) = g%dy_Cu(i,j) * h(i+1) * visc_rem(i)
          endif
        endif
      endif
    endif ; enddo
  endif
end subroutine zonal_flux_layer

!> Sets the effective interface thickness at each zonal velocity point.
subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, US, LB, vol_CFL, &
                                marginal, visc_rem_u, OBC)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< Zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thickness used to
                                                                   !! calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_u  !< Thickness at zonal faces [H ~> m or kg m-2].
  real,                                      intent(in)    :: dt   !< Time increment [T ~> s].
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  type(loop_bounds_type),                    intent(in)    :: LB   !< Loop bounds structure.
  logical,                                   intent(in)    :: vol_CFL !< If true, rescale the ratio
                          !! of face areas to the cell areas when estimating the CFL number.
  logical,                                   intent(in)    :: marginal !< If true, report the
                          !! marginal face thicknesses; otherwise report transport-averaged thicknesses.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                   optional, intent(in)    :: visc_rem_u
                          !< Both the fraction of the momentum originally in a layer that remains after
                          !! a time-step of viscosity, and the fraction of a time-step's worth of a
                          !! barotropic acceleration that a layer experiences after viscosity is applied.
                          !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  type(ocean_obc_type),            optional, pointer       :: OBC !< Open boundaries control structure.

  ! Local variables
  real :: CFL  ! The CFL number based on the local velocity and grid spacing [nondim]
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h_in.
  real :: h_avg  ! The average thickness of a flux [H ~> m or kg m-2].
  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
  logical :: local_open_BC
  integer :: i, j, k, ish, ieh, jsh, jeh, nz, n
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke

  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
  do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
    if (u(i,j,k) > 0.0) then
      if (vol_cfl) then ; cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
      else ; cfl = u(i,j,k) * dt * g%IdxT(i,j) ; endif
      curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
      h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
      h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + 3.0*curv_3*(cfl - 1.0))
    elseif (u(i,j,k) < 0.0) then
      if (vol_cfl) then ; cfl = (-u(i,j,k)*dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
      else ; cfl = -u(i,j,k) * dt * g%IdxT(i+1,j) ; endif
      curv_3 = h_l(i+1,j,k) + h_r(i+1,j,k) - 2.0*h(i+1,j,k)
      h_avg = h_l(i+1,j,k) + cfl * (0.5*(h_r(i+1,j,k)-h_l(i+1,j,k)) + curv_3*(cfl - 1.5))
      h_marg = h_l(i+1,j,k) + cfl * ((h_r(i+1,j,k)-h_l(i+1,j,k)) + &
                                    3.0*curv_3*(cfl - 1.0))
    else
      h_avg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
      !   The choice to use the arithmetic mean here is somewhat arbitrariy, but
      ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
      h_marg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
 !    h_marg = (2.0 * h_L(i+1,j,k) * h_R(i,j,k)) / &
 !             (h_L(i+1,j,k) + h_R(i,j,k) + GV%H_subroundoff)
    endif

    if (marginal) then ; h_u(i,j,k) = h_marg
    else ; h_u(i,j,k) = h_avg ; endif
  enddo ; enddo ; enddo
  if (present(visc_rem_u)) then
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
      h_u(i,j,k) = h_u(i,j,k) * visc_rem_u(i,j,k)
    enddo ; enddo ; enddo
  endif

  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
  if (local_open_bc) then
    do n = 1, obc%number_of_segments
      if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then
        i = obc%segment(n)%HI%IsdB
        if (obc%segment(n)%direction == obc_direction_e) then
          if (present(visc_rem_u)) then ; do k=1,nz
            do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
              h_u(i,j,k) = h(i,j,k) * visc_rem_u(i,j,k)
            enddo
          enddo ; else ; do k=1,nz
            do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
              h_u(i,j,k) = h(i,j,k)
            enddo
          enddo ; endif
        else
          if (present(visc_rem_u)) then ; do k=1,nz
            do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
              h_u(i,j,k) = h(i+1,j,k) * visc_rem_u(i,j,k)
            enddo
          enddo ; else ; do k=1,nz
            do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
              h_u(i,j,k) = h(i+1,j,k)
            enddo
          enddo ; endif
        endif
      endif
    enddo
  endif

end subroutine zonal_face_thickness

!> Returns the barotropic velocity adjustment that gives the
!! desired barotropic (layer-summed) transport.
subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, &
                             du, du_max_CFL, du_min_CFL, dt, G, US, CS, visc_rem, &
                             j, ish, ieh, do_I_in, full_precision, uh_3d, OBC)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< Zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZK_(G)),         intent(in)    :: visc_rem !< 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.
                       !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G)),       optional, intent(in)    :: uhbt !< The summed volume flux
                       !! through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].

  real, dimension(SZIB_(G)),                 intent(in)    :: du_max_CFL  !< Maximum acceptable
                       !! value of du [L T-1 ~> m s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: du_min_CFL  !< Minimum acceptable
                       !! value of du [L T-1 ~> m s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: uh_tot_0    !< The summed transport
                       !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: duhdu_tot_0 !< The partial derivative
                       !! of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].
  real, dimension(SZIB_(G)),                 intent(out)   :: du !<
                       !! The barotropic velocity adjustment [L T-1 ~> m s-1].
  real,                                      intent(in)    :: dt   !< Time increment [T ~> s].
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_ppm_cs),                   pointer       :: CS   !< This module's control structure.
  integer,                                   intent(in)    :: j    !< Spatial index.
  integer,                                   intent(in)    :: ish  !< Start of index range.
  integer,                                   intent(in)    :: ieh  !< End of index range.
  logical, dimension(SZIB_(G)),              intent(in)    :: do_I_in     !<
                       !! A logical flag indicating which I values to work on.
  logical,                         optional, intent(in)    :: full_precision !<
                       !! A flag indicating how carefully to iterate.  The
                       !! default is .true. (more accurate).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), optional, intent(inout) :: uh_3d !<
                       !! Volume flux through zonal faces = u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].
  type(ocean_obc_type),            optional, pointer       :: OBC !< Open boundaries control structure.
  ! Local variables
  real, dimension(SZIB_(G),SZK_(G)) :: &
    uh_aux, &  ! An auxiliary zonal volume flux [H L2 s-1 ~> m3 s-1 or kg s-1].
    duhdu      ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
  real, dimension(SZIB_(G)) :: &
    uh_err, &  ! Difference between uhbt and the summed uh [H L2 T-1 ~> m3 s-1 or kg s-1].
    uh_err_best, & ! The smallest value of uh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
    u_new, &   ! The velocity with the correction added [L T-1 ~> m s-1].
    duhdu_tot,&! Summed partial derivative of uh with u [H L ~> m2 or kg m-1].
    du_min, &  ! Min/max limits on du correction based on CFL limits
    du_max     ! and previous iterations [L T-1 ~> m s-1].
  real :: du_prev ! The previous value of du [L T-1 ~> m s-1].
  real :: ddu    ! The change in du from the previous iteration [L T-1 ~> m s-1].
  real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
  real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].
  integer :: i, k, nz, itt, max_itts = 20
  logical :: full_prec, domore, do_I(SZIB_(G))

  nz = g%ke
  full_prec = .true. ; if (present(full_precision)) full_prec = full_precision

  uh_aux(:,:) = 0.0 ; duhdu(:,:) = 0.0

  if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh
    uh_aux(i,k) = uh_3d(i,j,k)
  enddo ; enddo ; endif

  do i=ish-1,ieh
    du(i) = 0.0 ; do_i(i) = do_i_in(i)
    du_max(i) = du_max_cfl(i) ; du_min(i) = du_min_cfl(i)
    uh_err(i) = uh_tot_0(i) - uhbt(i) ; duhdu_tot(i) = duhdu_tot_0(i)
    uh_err_best(i) = abs(uh_err(i))
  enddo

  do itt=1,max_itts
    if (full_prec) then
      select case (itt)
        case (:1) ; tol_eta = 1e-6 * cs%tol_eta
        case (2)  ; tol_eta = 1e-4 * cs%tol_eta
        case (3)  ; tol_eta = 1e-2 * cs%tol_eta
        case default ; tol_eta = cs%tol_eta
      end select
    else
      tol_eta = cs%tol_eta_aux ; if (itt<=1) tol_eta = 1e-6 * cs%tol_eta_aux
    endif
    tol_vel = cs%tol_vel

    do i=ish-1,ieh
      if (uh_err(i) > 0.0) then ; du_max(i) = du(i)
      elseif (uh_err(i) < 0.0) then ; du_min(i) = du(i)
      else ; do_i(i) = .false. ; endif
    enddo
    domore = .false.
    do i=ish-1,ieh ; if (do_i(i)) then
      if ((dt * min(g%IareaT(i,j),g%IareaT(i+1,j))*abs(uh_err(i)) > tol_eta) .or. &
          (cs%better_iter .and. ((abs(uh_err(i)) > tol_vel * duhdu_tot(i)) .or. &
                                 (abs(uh_err(i)) > uh_err_best(i))) )) then
        !   Use Newton's method, provided it stays bounded.  Otherwise bisect
        ! the value with the appropriate bound.
        if (full_prec) then
          ddu = -uh_err(i) / duhdu_tot(i)
          du_prev = du(i)
          du(i) = du(i) + ddu
          if (abs(ddu) < 1.0e-15*abs(du(i))) then
            do_i(i) = .false. ! ddu is small enough to quit.
          elseif (ddu > 0.0) then
            if (du(i) >= du_max(i)) then
              du(i) = 0.5*(du_prev + du_max(i))
              if (du_max(i) - du_prev < 1.0e-15*abs(du(i))) do_i(i) = .false.
            endif
          else ! ddu < 0.0
            if (du(i) <= du_min(i)) then
              du(i) = 0.5*(du_prev + du_min(i))
              if (du_prev - du_min(i) < 1.0e-15*abs(du(i))) do_i(i) = .false.
            endif
          endif
        else
          !   Use Newton's method, provided it stays bounded, just like above.
          du(i) = du(i) - uh_err(i) / duhdu_tot(i)
          if ((du(i) >= du_max(i)) .or. (du(i) <= du_min(i))) &
            du(i) = 0.5*(du_max(i) + du_min(i))
        endif
        if (do_i(i)) domore = .true.
      else
        do_i(i) = .false.
      endif
    endif ; enddo
    if (.not.domore) exit

    if ((itt < max_itts) .or. present(uh_3d)) then ; do k=1,nz
      do i=ish-1,ieh ; u_new(i) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
      call zonal_flux_layer(u_new, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), &
                            uh_aux(:,k), duhdu(:,k), visc_rem(:,k), &
                            dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
    enddo ; endif

    if (itt < max_itts) then
      do i=ish-1,ieh
        uh_err(i) = -uhbt(i) ; duhdu_tot(i) = 0.0
      enddo
      do k=1,nz ; do i=ish-1,ieh
        uh_err(i) = uh_err(i) + uh_aux(i,k)
        duhdu_tot(i) = duhdu_tot(i) + duhdu(i,k)
      enddo ; enddo
      do i=ish-1,ieh
        uh_err_best(i) = min(uh_err_best(i), abs(uh_err(i)))
      enddo
    endif
  enddo ! itt-loop
  ! If there are any faces which have not converged to within the tolerance,
  ! so-be-it, or else use a final upwind correction?
  ! This never seems to happen with 20 iterations as max_itt.

  if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh
    uh_3d(i,j,k) = uh_aux(i,k)
  enddo ; enddo ; endif

end subroutine zonal_flux_adjust

!> Sets a structure that describes the zonal barotropic volume or mass fluxes as a
!! function of barotropic flow to agree closely with the sum of the layer's transports.
subroutine set_zonal_bt_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, &
                             du_max_CFL, du_min_CFL, dt, G, US, CS, visc_rem, &
                             visc_rem_max, j, ish, ieh, do_I)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< Zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  type(bt_cont_type),                        intent(inout) :: BT_cont !< A structure with elements
                       !! that describe the effective open face areas as a function of barotropic flow.
  real, dimension(SZIB_(G)),                 intent(in)    :: uh_tot_0    !< The summed transport
                       !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: duhdu_tot_0 !< The partial derivative
                       !! of du_err with du at 0 adjustment [H L ~> m2 or kg m-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: du_max_CFL  !< Maximum acceptable
                       !! value of du [L T-1 ~> m s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: du_min_CFL  !< Minimum acceptable
                       !! value of du [L T-1 ~> m s-1].
  real,                                      intent(in)    :: dt   !< Time increment [T ~> s].
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_ppm_cs),                   pointer       :: CS   !< This module's control structure.
  real, dimension(SZIB_(G),SZK_(G)),         intent(in)    :: visc_rem !< 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.
                       !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G)),                 intent(in)    :: visc_rem_max !< Maximum allowable visc_rem.
  integer,                                   intent(in)    :: j        !< Spatial index.
  integer,                                   intent(in)    :: ish      !< Start of index range.
  integer,                                   intent(in)    :: ieh      !< End of index range.
  logical, dimension(SZIB_(G)),              intent(in)    :: do_I     !< A logical flag indicating
                       !! which I values to work on.
  ! Local variables
  real, dimension(SZIB_(G)) :: &
    du0, &        ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
    duL, duR, &   ! The barotropic velocity increments that give the westerly
                  ! (duL) and easterly (duR) test velocities [L T-1 ~> m s-1].
    zeros, &      ! An array of full of 0's.
    du_cfl, &     ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
    u_l, u_r, &   ! The westerly (u_L), easterly (u_R), and zero-barotropic
    u_0, &        ! transport (u_0) layer test velocities [L T-1 ~> m s-1].
    duhdu_l, &    ! The effective layer marginal face areas with the westerly
    duhdu_r, &    ! (_L), easterly (_R), and zero-barotropic (_0) test
    duhdu_0, &    ! velocities [H L ~> m2 or kg m-1].
    uh_l, uh_r, & ! The layer transports with the westerly (_L), easterly (_R),
    uh_0, &       ! and zero-barotropic (_0) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
    famt_l, famt_r, & ! The summed effective marginal face areas for the 3
    famt_0, &     ! test velocities [H L ~> m2 or kg m-1].
    uhtot_l, &    ! The summed transport with the westerly (uhtot_L) and
    uhtot_r       ! and easterly (uhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
  real :: FA_0    ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m].
  real :: FA_avg  ! The average effective face area [L H ~> m2 or kg m], nominally given by
                  ! the realized transport divided by the barotropic velocity.
  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]  This
                       ! limiting is necessary to keep the inverse of visc_rem
                       ! from leading to large CFL numbers.
  real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
                       ! in finding the barotropic velocity that changes the
                       ! flow direction.  This is necessary to keep the inverse
                       ! of visc_rem from leading to large CFL numbers.
  real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
                  ! flow is truly upwind [nondim]
  real :: Idt     ! The inverse of the time step [T-1 ~> s-1].
  logical :: domore
  integer :: i, k, nz

  nz = g%ke ; idt = 1.0 / dt
  min_visc_rem = 0.1 ; cfl_min = 1e-6

 ! Diagnose the zero-transport correction, du0.
  do i=ish-1,ieh ; zeros(i) = 0.0 ; enddo
  call zonal_flux_adjust(u, h_in, h_l, h_r, zeros, uh_tot_0, duhdu_tot_0, du0, &
                         du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
                         j, ish, ieh, do_i, .true.)

  ! Determine the westerly- and easterly- fluxes.  Choose a sufficiently
  ! negative velocity correction for the easterly-flux, and a sufficiently
  ! positive correction for the westerly-flux.
  domore = .false.
  do i=ish-1,ieh
    if (do_i(i)) domore = .true.
    du_cfl(i) = (cfl_min * idt) * g%dxCu(i,j)
    dur(i) = min(0.0,du0(i) - du_cfl(i))
    dul(i) = max(0.0,du0(i) + du_cfl(i))
    famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0
    uhtot_l(i) = 0.0 ; uhtot_r(i) = 0.0
  enddo

  if (.not.domore) then
    do k=1,nz ; do i=ish-1,ieh
      bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
      bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
      bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
    enddo ; enddo
    return
  endif

  do k=1,nz ; do i=ish-1,ieh ; if (do_i(i)) then
    visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
    if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
      if (u(i,j,k) + dur(i)*visc_rem_lim > -du_cfl(i)*visc_rem(i,k)) &
        dur(i) = -(u(i,j,k) + du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
      if (u(i,j,k) + dul(i)*visc_rem_lim < du_cfl(i)*visc_rem(i,k)) &
        dul(i) = -(u(i,j,k) - du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
    endif
  endif ; enddo ; enddo

  do k=1,nz
    do i=ish-1,ieh ; if (do_i(i)) then
      u_l(i) = u(i,j,k) + dul(i) * visc_rem(i,k)
      u_r(i) = u(i,j,k) + dur(i) * visc_rem(i,k)
      u_0(i) = u(i,j,k) + du0(i) * visc_rem(i,k)
    endif ; enddo
    call zonal_flux_layer(u_0, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_0, duhdu_0, &
                          visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
    call zonal_flux_layer(u_l, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_l, duhdu_l, &
                          visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
    call zonal_flux_layer(u_r, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_r, duhdu_r, &
                          visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
    do i=ish-1,ieh ; if (do_i(i)) then
      famt_0(i) = famt_0(i) + duhdu_0(i)
      famt_l(i) = famt_l(i) + duhdu_l(i)
      famt_r(i) = famt_r(i) + duhdu_r(i)
      uhtot_l(i) = uhtot_l(i) + uh_l(i)
      uhtot_r(i) = uhtot_r(i) + uh_r(i)
    endif ; enddo
  enddo
  do i=ish-1,ieh ; if (do_i(i)) then
    fa_0 = famt_0(i) ; fa_avg = famt_0(i)
    if ((dul(i) - du0(i)) /= 0.0) &
      fa_avg = uhtot_l(i) / (dul(i) - du0(i))
    if (fa_avg > max(fa_0, famt_l(i))) then ; fa_avg = max(fa_0, famt_l(i))
    elseif (fa_avg < min(fa_0, famt_l(i))) then ; fa_0 = fa_avg ; endif

    bt_cont%FA_u_W0(i,j) = fa_0 ; bt_cont%FA_u_WW(i,j) = famt_l(i)
    if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0) then ; bt_cont%uBT_WW(i,j) = 0.0 ; else
      bt_cont%uBT_WW(i,j) = (1.5 * (dul(i) - du0(i))) * &
                            ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))
    endif

    fa_0 = famt_0(i) ; fa_avg = famt_0(i)
    if ((dur(i) - du0(i)) /= 0.0) &
      fa_avg = uhtot_r(i) / (dur(i) - du0(i))
    if (fa_avg > max(fa_0, famt_r(i))) then ; fa_avg = max(fa_0, famt_r(i))
    elseif (fa_avg < min(fa_0, famt_r(i))) then ; fa_0 = fa_avg ; endif

    bt_cont%FA_u_E0(i,j) = fa_0 ; bt_cont%FA_u_EE(i,j) = famt_r(i)
    if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0) then ; bt_cont%uBT_EE(i,j) = 0.0 ; else
      bt_cont%uBT_EE(i,j) = (1.5 * (dur(i) - du0(i))) * &
                            ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
    endif
  else
    bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
    bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
    bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
  endif ; enddo

end subroutine set_zonal_bt_cont

!> Calculates the mass or volume fluxes through the meridional faces, and other related quantities.
subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, &
                                visc_rem_v, v_cor, BT_cont)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  type(verticalgrid_type),                   intent(in)    :: GV   !< Ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< Meridional velocity [L T-1 ~> m s-1].
  real,  dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: vh   !< Volume flux through meridional
                                                                   !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
  real,                                      intent(in)    :: dt   !< Time increment [T ~> s].
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_ppm_cs),                   pointer       :: CS   !< This module's control structure.G
  type(loop_bounds_type),                    intent(in)    :: LB   !< Loop bounds structure.
  type(ocean_obc_type),            optional, pointer       :: OBC  !< Open boundary condition type
                                   !! specifies whether, where, and what open boundary conditions are used.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                   optional, intent(in)    :: 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.  Nondimensional between
                                   !! 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJB_(G)), optional, intent(in)  :: vhbt  !< The summed volume flux through
                                   !< meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                     optional, intent(out) :: v_cor
                                   !< The meridional velocitiess (v with a barotropic correction)
                                   !! that give vhbt as the depth-integrated transport [L T-1 ~> m s-1].
  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.
  ! Local variables
  real, dimension(SZI_(G),SZK_(G)) :: &
    dvhdv      ! Partial derivative of vh with v [H L ~> m2 or kg m-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    h_L, h_R   ! Left and right face thicknesses [H ~> m or kg m-2].
  real, dimension(SZI_(G)) :: &
    dv, &      ! Corrective barotropic change in the velocity [L T-1 ~> m s-1].
    dv_min_CFL, & ! Min/max limits on dv correction
    dv_max_CFL, & ! to avoid CFL violations
    dvhdv_tot_0, & ! Summed partial derivative of vh with v [H L ~> m2 or kg m-1].
    vh_tot_0, &   ! Summed transport with no barotropic correction [H L2 T-1 ~> m3 s-1 or kg s-1].
    visc_rem_max  ! The column maximum of visc_rem.
  logical, dimension(SZI_(G)) :: do_I
  real, dimension(SZI_(G)) :: FAvi  ! A list of sums of meridional face areas [H L ~> m2 or kg m-1].
  real :: FA_v    ! A sum of meridional face areas [H m ~> m2 or kg m-1].
  real, dimension(SZI_(G),SZK_(G)) :: &
    visc_rem      ! A 2-D copy of visc_rem_v or an array of 1's.
  real :: I_vrm   ! 1.0 / visc_rem_max, nondim.
  real :: CFL_dt  ! The maximum CFL ratio of the adjusted velocities divided by
                  ! the time step [T-1 ~> s-1].
  real :: I_dt    ! 1.0 / dt [T-1 ~> s-1].
  real :: dv_lim  ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1].
  real :: dy_N, dy_S ! Effective y-grid spacings to the north and south [L ~> m].
  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
  integer :: l_seg
  logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC
  logical :: local_Flather_OBC, is_simple, local_open_BC
  type(obc_segment_type), pointer :: segment => null()

  use_visc_rem = present(visc_rem_v)
  local_specified_bc = .false. ; set_bt_cont = .false. ; local_flather_obc = .false.
  local_open_bc = .false.
  if (present(bt_cont)) set_bt_cont = (associated(bt_cont))
  if (present(obc)) then ; if (associated(obc)) then ; if (obc%OBC_pe) then
    local_specified_bc = obc%specified_v_BCs_exist_globally
    local_flather_obc = obc%Flather_v_BCs_exist_globally
    local_open_bc = obc%open_v_BCs_exist_globally
  endif ; endif ; endif
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke

  cfl_dt = cs%CFL_limit_adjust / dt
  i_dt = 1.0 / dt
  if (cs%aggress_adjust) cfl_dt = i_dt

  call cpu_clock_begin(id_clock_update)
!$OMP parallel do default(none) shared(nz,ish,ieh,jsh,jeh,h_in,h_L,h_R,G,GV,LB,CS,visc_rem,OBC)
  do k=1,nz
    ! This sets h_L and h_R.
    if (cs%upwind_1st) then
      do j=jsh-1,jeh+1 ; do i=ish,ieh
        h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
      enddo ; enddo
    else
      call ppm_reconstruction_y(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
                                2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
    endif
    do i=ish,ieh ; visc_rem(i,k) = 1.0 ; enddo
  enddo
  call cpu_clock_end(id_clock_update)

  call cpu_clock_begin(id_clock_correct)
!$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,v,h_in,h_L,h_R,vh,use_visc_rem, &
!$OMP                                  visc_rem_v,dt,US,G,GV,CS,local_specified_BC,OBC,vhbt, &
!$OMP                                  set_BT_cont,CFL_dt,I_dt,v_cor,BT_cont, local_Flather_OBC ) &
!$OMP                          private(do_I,dvhdv,dv,dv_max_CFL,dv_min_CFL,vh_tot_0,    &
!$OMP                                  dvhdv_tot_0,visc_rem_max,I_vrm,dv_lim,dy_N,      &
!$OMP                                  is_simple,FAvi,dy_S,any_simple_OBC,l_seg) &
!$OMP                     firstprivate(visc_rem)
  do j=jsh-1,jeh
    do i=ish,ieh ; do_i(i) = .true. ; visc_rem_max(i) = 0.0 ; enddo
    ! This sets vh and dvhdv.
    do k=1,nz
      if (use_visc_rem) then ; do i=ish,ieh
        visc_rem(i,k) = visc_rem_v(i,j,k)
        visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))
      enddo ; endif
      call merid_flux_layer(v(:,j,k), h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), &
                            vh(:,j,k), dvhdv(:,k), visc_rem(:,k), &
                            dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
      if (local_specified_bc) then
        do i=ish,ieh
          l_seg = obc%segnum_v(i,j)

          if (l_seg /= obc_none) then
            if (obc%segment(l_seg)%specified) &
              vh(i,j,k) = obc%segment(l_seg)%normal_trans(i,j,k)
          endif
        enddo
      endif
    enddo ! k-loop
    if ((.not.use_visc_rem) .or. (.not.cs%use_visc_rem_max)) then ; do i=ish,ieh
      visc_rem_max(i) = 1.0
    enddo ; endif

    if (present(vhbt) .or. set_bt_cont) then
      !   Set limits on dv that will keep the CFL number between -1 and 1.
      ! This should be adequate to keep the root bracketed in all cases.
      do i=ish,ieh
        i_vrm = 0.0
        if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)
        if (cs%vol_CFL) then
          dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
          dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
        else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
        dv_max_cfl(i) = 2.0 * (cfl_dt * dy_s) * i_vrm
        dv_min_cfl(i) = -2.0 * (cfl_dt * dy_n) * i_vrm
        vh_tot_0(i) = 0.0 ; dvhdv_tot_0(i) = 0.0
      enddo
      do k=1,nz ; do i=ish,ieh
        dvhdv_tot_0(i) = dvhdv_tot_0(i) + dvhdv(i,k)
        vh_tot_0(i) = vh_tot_0(i) + vh(i,j,k)
      enddo ; enddo

      if (use_visc_rem) then
        if (cs%aggress_adjust) then
          do k=1,nz ; do i=ish,ieh
            if (cs%vol_CFL) then
              dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
              dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
            else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
            dv_lim = 0.499*((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k)))
            if (dv_max_cfl(i) * visc_rem(i,k) > dv_lim) &
              dv_max_cfl(i) = dv_lim / visc_rem(i,k)

            dv_lim = 0.499*((-dy_n*cfl_dt - v(i,j,k)) + max(0.0,v(i,j+1,k)))
            if (dv_min_cfl(i) * visc_rem(i,k) < dv_lim) &
              dv_min_cfl(i) = dv_lim / visc_rem(i,k)
          enddo ; enddo
        else
          do k=1,nz ; do i=ish,ieh
            if (cs%vol_CFL) then
              dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
              dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
            else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
            if (dv_max_cfl(i) * visc_rem(i,k) > dy_s*cfl_dt - v(i,j,k)) &
              dv_max_cfl(i) = (dy_s*cfl_dt - v(i,j,k)) / visc_rem(i,k)
            if (dv_min_cfl(i) * visc_rem(i,k) < -dy_n*cfl_dt - v(i,j,k)) &
              dv_min_cfl(i) = -(dy_n*cfl_dt + v(i,j,k)) / visc_rem(i,k)
          enddo ; enddo
        endif
      else
        if (cs%aggress_adjust) then
          do k=1,nz ; do i=ish,ieh
            if (cs%vol_CFL) then
              dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
              dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
            else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
            dv_max_cfl(i) = min(dv_max_cfl(i), 0.499 * &
                        ((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k))) )
            dv_min_cfl(i) = max(dv_min_cfl(i), 0.499 * &
                        ((-dy_n*i_dt - v(i,j,k)) + max(0.0,v(i,j+1,k))) )
          enddo ; enddo
        else
          do k=1,nz ; do i=ish,ieh
            if (cs%vol_CFL) then
              dy_s = ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
              dy_n = ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
            else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ; endif
            dv_max_cfl(i) = min(dv_max_cfl(i), dy_s*cfl_dt - v(i,j,k))
            dv_min_cfl(i) = max(dv_min_cfl(i), -(dy_n*cfl_dt + v(i,j,k)))
          enddo ; enddo
        endif
      endif
      do i=ish,ieh
        dv_max_cfl(i) = max(dv_max_cfl(i),0.0)
        dv_min_cfl(i) = min(dv_min_cfl(i),0.0)
      enddo

      any_simple_obc = .false.
      if (present(vhbt) .or. set_bt_cont) then
        if (local_specified_bc .or. local_flather_obc) then ; do i=ish,ieh
          l_seg = obc%segnum_v(i,j)

          ! Avoid reconciling barotropic/baroclinic transports if transport is specified
          is_simple = .false.
          if (l_seg /= obc_none) &
            is_simple = obc%segment(l_seg)%specified
          do_i(i) = .not.(l_seg /= obc_none .and. is_simple)
          any_simple_obc = any_simple_obc .or. is_simple
        enddo ; else ; do i=ish,ieh
          do_i(i) = .true.
        enddo ; endif
      endif

      if (present(vhbt)) then
        call meridional_flux_adjust(v, h_in, h_l, h_r, vhbt(:,j), vh_tot_0, dvhdv_tot_0, dv, &
                               dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
                               j, ish, ieh, do_i, .true., vh, obc=obc)

        if (present(v_cor)) then ; do k=1,nz
          do i=ish,ieh ; v_cor(i,j,k) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
          if (local_specified_bc) then ; do i=ish,ieh
            l_seg = obc%segnum_v(i,j)

            if (l_seg /= obc_none) then
              if (obc%segment(obc%segnum_v(i,j))%specified) &
                v_cor(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
            endif
          enddo ; endif
        enddo ; endif ! v-corrected
      endif

      if (set_bt_cont) then
        call set_merid_bt_cont(v, h_in, h_l, h_r, bt_cont, vh_tot_0, dvhdv_tot_0,&
                               dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
                               visc_rem_max, j, ish, ieh, do_i)
        if (any_simple_obc) then
          do i=ish,ieh
            l_seg = obc%segnum_v(i,j)

            do_i(i) = .false.
            if(l_seg /= obc_none) &
              do_i(i) = (obc%segment(l_seg)%specified)

            if (do_i(i)) favi(i) = gv%H_subroundoff*g%dx_Cv(i,j)
          enddo
          ! NOTE: do_I(I) should prevent access to segment OBC_NONE
          do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then
            if ((abs(obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)) > 0.0) .and. &
                (obc%segment(obc%segnum_v(i,j))%specified)) &
              favi(i) = favi(i) + obc%segment(obc%segnum_v(i,j))%normal_trans(i,j,k) / &
                                  obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
          endif ; enddo ; enddo
          do i=ish,ieh ; if (do_i(i)) then
            bt_cont%FA_v_S0(i,j) = favi(i) ; bt_cont%FA_v_N0(i,j) = favi(i)
            bt_cont%FA_v_SS(i,j) = favi(i) ; bt_cont%FA_v_NN(i,j) = favi(i)
            bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
          endif ; enddo
        endif
      endif ! set_BT_cont

    endif ! present(vhbt) or set_BT_cont

  enddo ! j-loop

  if (local_open_bc .and. set_bt_cont) then
    do n = 1, obc%number_of_segments
      if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
        j = obc%segment(n)%HI%JsdB
        if (obc%segment(n)%direction == obc_direction_n) then
          do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
            fa_v = 0.0
            do k=1,nz ; fa_v = fa_v + h_in(i,j,k)*g%dx_Cv(i,j) ; enddo
            bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
            bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
            bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
          enddo
        else
          do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
            fa_v = 0.0
            do k=1,nz ; fa_v = fa_v + h_in(i,j+1,k)*g%dx_Cv(i,j) ; enddo
            bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
            bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
            bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
          enddo
        endif
      endif
    enddo
  endif
  call cpu_clock_end(id_clock_correct)

  if (set_bt_cont) then ; if (allocated(bt_cont%h_v)) then
    if (present(v_cor)) then
      call merid_face_thickness(v_cor, h_in, h_l, h_r, bt_cont%h_v, dt, g, us, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
    else
      call merid_face_thickness(v, h_in, h_l, h_r, bt_cont%h_v, dt, g, us, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
    endif
  endif ; endif

end subroutine meridional_mass_flux

!> Evaluates the meridional mass or volume fluxes in a layer.
subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, &
                            ish, ieh, do_I, vol_CFL, OBC)
  type(ocean_grid_type),        intent(inout) :: G        !< Ocean's grid structure.
  real, dimension(SZI_(G)),     intent(in)    :: v        !< Meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G)),     intent(in)    :: visc_rem !< 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.
         !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJ_(G)),  intent(in) :: h      !< Layer thickness used to calculate fluxes,
                                                          !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(in) :: h_L    !< Left thickness in the reconstruction
                                                          !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(in) :: h_R    !< Right thickness in the reconstruction
                                                          !! [H ~> m or kg m-2].
  real, dimension(SZI_(G)),     intent(inout) :: vh       !< Meridional mass or volume transport
                                                          !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G)),     intent(inout) :: dvhdv    !< Partial derivative of vh with v
                                                          !! [H L ~> m2 or kg m-1].
  real,                         intent(in)    :: dt       !< Time increment [T ~> s].
  type(unit_scale_type),        intent(in)    :: US       !< A dimensional unit scaling type
  integer,                      intent(in)    :: j        !< Spatial index.
  integer,                      intent(in)    :: ish      !< Start of index range.
  integer,                      intent(in)    :: ieh      !< End of index range.
  logical, dimension(SZI_(G)),  intent(in)    :: do_I     !< Which i values to work on.
  logical,                      intent(in)    :: vol_CFL  !< If true, rescale the
         !! ratio of face areas to the cell areas when estimating the CFL number.
  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
  ! Local variables
  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h, i.e. [H ~> m or kg m-2].
  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
  integer :: i
  integer :: l_seg
  logical :: local_open_BC

  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_v_BCs_exist_globally
  endif ; endif

  do i=ish,ieh ; if (do_i(i)) then
    if (v(i) > 0.0) then
      if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
      else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
      curv_3 = h_l(i,j) + h_r(i,j) - 2.0*h(i,j)
      vh(i) = g%dx_Cv(i,j) * v(i) * ( h_r(i,j) + cfl * &
          (0.5*(h_l(i,j) - h_r(i,j)) + curv_3*(cfl - 1.5)) )
      h_marg = h_r(i,j) + cfl * ((h_l(i,j) - h_r(i,j)) + &
                                  3.0*curv_3*(cfl - 1.0))
    elseif (v(i) < 0.0) then
      if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
      else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
      curv_3 = h_l(i,j+1) + h_r(i,j+1) - 2.0*h(i,j+1)
      vh(i) = g%dx_Cv(i,j) * v(i) * ( h_l(i,j+1) + cfl * &
          (0.5*(h_r(i,j+1)-h_l(i,j+1)) + curv_3*(cfl - 1.5)) )
      h_marg = h_l(i,j+1) + cfl * ((h_r(i,j+1)-h_l(i,j+1)) + &
                                    3.0*curv_3*(cfl - 1.0))
    else
      vh(i) = 0.0
      h_marg = 0.5 * (h_l(i,j+1) + h_r(i,j))
    endif
    dvhdv(i) = g%dx_Cv(i,j) * h_marg * visc_rem(i)
  endif ; enddo

  if (local_open_bc) then
    do i=ish,ieh ; if (do_i(i)) then
      l_seg = obc%segnum_v(i,j)

      if (l_seg /= obc_none) then
        if (obc%segment(l_seg)%open) then
          if (obc%segment(l_seg)%direction == obc_direction_n) then
            vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j)
            dvhdv(i) = g%dx_Cv(i,j) * h(i,j) * visc_rem(i)
          else
            vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j+1)
            dvhdv(i) = g%dx_Cv(i,j) * h(i,j+1) * visc_rem(i)
          endif
        endif
      endif
    endif ; enddo
  endif
end subroutine merid_flux_layer

!> Sets the effective interface thickness at each meridional velocity point.
subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, US, LB, vol_CFL, &
                                marginal, visc_rem_v, OBC)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< Meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thickness used to calculate fluxes,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the reconstruction,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the reconstruction,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: h_v  !< Thickness at meridional faces,
                                                                   !! [H ~> m or kg m-2].
  real,                                      intent(in)    :: dt   !< Time increment [T ~> s].
  type(loop_bounds_type),                    intent(in)    :: LB   !< Loop bounds structure.
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  logical,                                   intent(in)    :: vol_CFL !< If true, rescale the ratio
                          !! of face areas to the cell areas when estimating the CFL number.
  logical,                                   intent(in)    :: marginal !< If true, report the marginal
                          !! face thicknesses; otherwise report transport-averaged thicknesses.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), optional, intent(in) :: 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.
                          !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  type(ocean_obc_type),            optional, pointer :: OBC !< Open boundaries control structure.

  ! Local variables
  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h [H ~> m or kg m-2] .
  real :: h_avg  ! The average thickness of a flux [H ~> m or kg m-2].
  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
  logical :: local_open_BC
  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke

  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
  do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
    if (v(i,j,k) > 0.0) then
      if (vol_cfl) then ; cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
      else ; cfl = v(i,j,k) * dt * g%IdyT(i,j) ; endif
      curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
      h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
      h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + &
                                3.0*curv_3*(cfl - 1.0))
    elseif (v(i,j,k) < 0.0) then
      if (vol_cfl) then ; cfl = (-v(i,j,k)*dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
      else ; cfl = -v(i,j,k) * dt * g%IdyT(i,j+1) ; endif
      curv_3 = h_l(i,j+1,k) + h_r(i,j+1,k) - 2.0*h(i,j+1,k)
      h_avg = h_l(i,j+1,k) + cfl * (0.5*(h_r(i,j+1,k)-h_l(i,j+1,k)) + curv_3*(cfl - 1.5))
      h_marg = h_l(i,j+1,k) + cfl * ((h_r(i,j+1,k)-h_l(i,j+1,k)) + &
                                    3.0*curv_3*(cfl - 1.0))
    else
      h_avg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
      !   The choice to use the arithmetic mean here is somewhat arbitrariy, but
      ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
      h_marg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
 !    h_marg = (2.0 * h_L(i,j+1,k) * h_R(i,j,k)) / &
 !             (h_L(i,j+1,k) + h_R(i,j,k) + GV%H_subroundoff)
    endif

    if (marginal) then ; h_v(i,j,k) = h_marg
    else ; h_v(i,j,k) = h_avg ; endif
  enddo ; enddo ; enddo

  if (present(visc_rem_v)) then
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
      h_v(i,j,k) = h_v(i,j,k) * visc_rem_v(i,j,k)
    enddo ; enddo ; enddo
  endif

  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_v_BCs_exist_globally
  endif ; endif
  if (local_open_bc) then
    do n = 1, obc%number_of_segments
      if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
        j = obc%segment(n)%HI%JsdB
        if (obc%segment(n)%direction == obc_direction_n) then
          if (present(visc_rem_v)) then ; do k=1,nz
            do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
              h_v(i,j,k) = h(i,j,k) * visc_rem_v(i,j,k)
            enddo
          enddo ; else ; do k=1,nz
            do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
              h_v(i,j,k) = h(i,j,k)
            enddo
          enddo ; endif
        else
          if (present(visc_rem_v)) then ; do k=1,nz
            do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
              h_v(i,j,k) = h(i,j+1,k) * visc_rem_v(i,j,k)
            enddo
          enddo ; else ; do k=1,nz
            do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
              h_v(i,j,k) = h(i,j+1,k)
            enddo
          enddo ; endif
        endif
      endif
    enddo
  endif

end subroutine merid_face_thickness

!> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.
subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0, &
                             dv, dv_max_CFL, dv_min_CFL, dt, G, US, CS, visc_rem, &
                             j, ish, ieh, do_I_in, full_precision, vh_3d, OBC)
  type(ocean_grid_type),  intent(inout) :: G   !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                         intent(in)    :: v    !< Meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                         intent(in)    :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),&
                         intent(in)    :: h_L  !< Left thickness in the reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                         intent(in)    :: h_R  !< Right thickness in the reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)), intent(in) :: visc_rem
                             !< 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. Non-dimensional
                             !! between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G)), &
                  optional, intent(in)    :: vhbt !< The summed volume flux through meridional faces
                                                  !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G)), intent(in)    :: dv_max_CFL !< Maximum acceptable value of dv [L T-1 ~> m s-1].
  real, dimension(SZI_(G)), intent(in)    :: dv_min_CFL !< Minimum acceptable value of dv [L T-1 ~> m s-1].
  real, dimension(SZI_(G)), intent(in)    :: vh_tot_0   !< The summed transport with 0 adjustment
                                                        !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G)), intent(in)    :: dvhdv_tot_0 !< The partial derivative of dv_err with
                                                         !! dv at 0 adjustment [H L ~> m2 or kg m-1].
  real, dimension(SZI_(G)), intent(out)   :: dv   !< The barotropic velocity adjustment [L T-1 ~> m s-1].
  real,                     intent(in)    :: dt   !< Time increment [T ~> s].
  type(unit_scale_type),    intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_ppm_cs),  pointer       :: CS   !< This module's control structure.
  integer,                  intent(in)    :: j    !< Spatial index.
  integer,                  intent(in)    :: ish  !< Start of index range.
  integer,                  intent(in)    :: ieh  !< End of index range.
  logical, dimension(SZI_(G)), &
                            intent(in)    :: do_I_in  !< A flag indicating which I values to work on.
  logical,        optional, intent(in)    :: full_precision !< A flag indicating how carefully to
                             !! iterate.  The default is .true. (more accurate).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                  optional, intent(inout) :: vh_3d !< Volume flux through meridional
                             !! faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].
  type(ocean_obc_type), optional, pointer :: OBC !< Open boundaries control structure.
  ! Local variables
  real, dimension(SZI_(G),SZK_(G)) :: &
    vh_aux, &  ! An auxiliary meridional volume flux [H L2 s-1 ~> m3 s-1 or kg s-1].
    dvhdv      ! Partial derivative of vh with v [H m ~> m2 or kg m-1].
  real, dimension(SZI_(G)) :: &
    vh_err, &  ! Difference between vhbt and the summed vh [H L2 T-1 ~> m3 s-1 or kg s-1].
    vh_err_best, & ! The smallest value of vh_err found so far [H L2 T-1 ~> m3 s-1 or kg s-1].
    v_new, &   ! The velocity with the correction added [L T-1 ~> m s-1].
    dvhdv_tot,&! Summed partial derivative of vh with u [H L ~> m2 or kg m-1].
    dv_min, &  ! Min/max limits on dv correction based on CFL limits
    dv_max     ! and previous iterations [L T-1 ~> m s-1].
  real :: dv_prev ! The previous value of dv [L T-1 ~> m s-1].
  real :: ddv    ! The change in dv from the previous iteration [L T-1 ~> m s-1].
  real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
  real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].
  integer :: i, k, nz, itt, max_itts = 20
  logical :: full_prec, domore, do_I(SZI_(G))

  nz = g%ke
  full_prec = .true. ; if (present(full_precision)) full_prec = full_precision

  vh_aux(:,:) = 0.0 ; dvhdv(:,:) = 0.0

  if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh
    vh_aux(i,k) = vh_3d(i,j,k)
  enddo ; enddo ; endif

  do i=ish,ieh
    dv(i) = 0.0 ; do_i(i) = do_i_in(i)
    dv_max(i) = dv_max_cfl(i) ; dv_min(i) = dv_min_cfl(i)
    vh_err(i) = vh_tot_0(i) - vhbt(i) ; dvhdv_tot(i) = dvhdv_tot_0(i)
    vh_err_best(i) = abs(vh_err(i))
  enddo

  do itt=1,max_itts
    if (full_prec) then
      select case (itt)
        case (:1) ; tol_eta = 1e-6 * cs%tol_eta
        case (2)  ; tol_eta = 1e-4 * cs%tol_eta
        case (3)  ; tol_eta = 1e-2 * cs%tol_eta
        case default ; tol_eta = cs%tol_eta
      end select
    else
      tol_eta = cs%tol_eta_aux ; if (itt<=1) tol_eta = 1e-6 * cs%tol_eta_aux
    endif
    tol_vel = cs%tol_vel

    do i=ish,ieh
      if (vh_err(i) > 0.0) then ; dv_max(i) = dv(i)
      elseif (vh_err(i) < 0.0) then ; dv_min(i) = dv(i)
      else ; do_i(i) = .false. ; endif
    enddo
    domore = .false.
    do i=ish,ieh ; if (do_i(i)) then
      if ((dt * min(g%IareaT(i,j),g%IareaT(i,j+1))*abs(vh_err(i)) > tol_eta) .or. &
          (cs%better_iter .and. ((abs(vh_err(i)) > tol_vel * dvhdv_tot(i)) .or. &
                                 (abs(vh_err(i)) > vh_err_best(i))) )) then
        !   Use Newton's method, provided it stays bounded.  Otherwise bisect
        ! the value with the appropriate bound.
        if (full_prec) then
          ddv = -vh_err(i) / dvhdv_tot(i)
          dv_prev = dv(i)
          dv(i) = dv(i) + ddv
          if (abs(ddv) < 1.0e-15*abs(dv(i))) then
            do_i(i) = .false. ! ddv is small enough to quit.
          elseif (ddv > 0.0) then
            if (dv(i) >= dv_max(i)) then
              dv(i) = 0.5*(dv_prev + dv_max(i))
              if (dv_max(i) - dv_prev < 1.0e-15*abs(dv(i))) do_i(i) = .false.
            endif
          else ! dvv(i) < 0.0
            if (dv(i) <= dv_min(i)) then
              dv(i) = 0.5*(dv_prev + dv_min(i))
              if (dv_prev - dv_min(i) < 1.0e-15*abs(dv(i))) do_i(i) = .false.
            endif
          endif
        else
          !   Use Newton's method, provided it stays bounded, just like above.
          dv(i) = dv(i) - vh_err(i) / dvhdv_tot(i)
          if ((dv(i) >= dv_max(i)) .or. (dv(i) <= dv_min(i))) &
            dv(i) = 0.5*(dv_max(i) + dv_min(i))
        endif
        if (do_i(i)) domore = .true.
      else
        do_i(i) = .false.
      endif
    endif ; enddo
    if (.not.domore) exit

    if ((itt < max_itts) .or. present(vh_3d)) then ; do k=1,nz
      do i=ish,ieh ; v_new(i) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
      call merid_flux_layer(v_new, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), &
                            vh_aux(:,k), dvhdv(:,k), visc_rem(:,k), &
                            dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
    enddo ; endif

    if (itt < max_itts) then
      do i=ish,ieh
        vh_err(i) = -vhbt(i) ; dvhdv_tot(i) = 0.0
      enddo
      do k=1,nz ; do i=ish,ieh
        vh_err(i) = vh_err(i) + vh_aux(i,k)
        dvhdv_tot(i) = dvhdv_tot(i) + dvhdv(i,k)
      enddo ; enddo
      do i=ish,ieh
        vh_err_best(i) = min(vh_err_best(i), abs(vh_err(i)))
      enddo
    endif
  enddo ! itt-loop
  ! If there are any faces which have not converged to within the tolerance,
  ! so-be-it, or else use a final upwind correction?
  ! This never seems to happen with 20 iterations as max_itt.

  if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh
    vh_3d(i,j,k) = vh_aux(i,k)
  enddo ; enddo ; endif

end subroutine meridional_flux_adjust

!> Sets of a structure that describes the meridional barotropic volume or mass fluxes as a
!! function of barotropic flow to agree closely with the sum of the layer's transports.
subroutine set_merid_bt_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, &
                             dv_max_CFL, dv_min_CFL, dt, G, US, CS, visc_rem, &
                             visc_rem_max, j, ish, ieh, do_I)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< Meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to calculate fluxes,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the reconstruction,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the reconstruction,
                                                                   !! [H ~> m or kg m-2].
  type(bt_cont_type),                        intent(inout) :: BT_cont !< A structure with elements
                       !! that describe the effective open face areas as a function of barotropic flow.
  real, dimension(SZI_(G)),                  intent(in)    :: vh_tot_0    !< The summed transport
                       !! with 0 adjustment [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G)),                  intent(in)    :: dvhdv_tot_0 !< The partial derivative
                       !! of du_err with dv at 0 adjustment [H L ~> m2 or kg m-1].
  real, dimension(SZI_(G)),                  intent(in)    :: dv_max_CFL !< Maximum acceptable value
                                                                   !!  of dv [L T-1 ~> m s-1].
  real, dimension(SZI_(G)),                  intent(in)    :: dv_min_CFL !< Minimum acceptable value
                                                                   !!  of dv [L T-1 ~> m s-1].
  real,                                      intent(in)    :: dt   !< Time increment [T ~> s].
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_ppm_cs),                   pointer       :: CS   !< This module's control structure.
  real, dimension(SZI_(G),SZK_(G)),          intent(in)    :: visc_rem !< 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.
                       !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G)),                  intent(in)    :: visc_rem_max !< Maximum allowable visc_rem.
  integer,                                   intent(in)    :: j        !< Spatial index.
  integer,                                   intent(in)    :: ish      !< Start of index range.
  integer,                                   intent(in)    :: ieh      !< End of index range.
  logical, dimension(SZI_(G)),               intent(in)    :: do_I     !< A logical flag indicating
                       !! which I values to work on.
  ! Local variables
  real, dimension(SZI_(G)) :: &
    dv0, &        ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1].
    dvL, dvR, &   ! The barotropic velocity increments that give the southerly
                  ! (dvL) and northerly (dvR) test velocities [L T-1 ~> m s-1].
    zeros, &      ! An array of full of 0's.
    dv_cfl, &     ! The velocity increment that corresponds to CFL_min [L T-1 ~> m s-1].
    v_l, v_r, &   ! The southerly (v_L), northerly (v_R), and zero-barotropic
    v_0, &        ! transport (v_0) layer test velocities [L T-1 ~> m s-1].
    dvhdv_l, &    ! The effective layer marginal face areas with the southerly
    dvhdv_r, &    ! (_L), northerly (_R), and zero-barotropic (_0) test
    dvhdv_0, &    ! velocities [H L ~> m2 or kg m-1].
    vh_l, vh_r, & ! The layer transports with the southerly (_L), northerly (_R)
    vh_0, &       ! and zero-barotropic (_0) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
    famt_l, famt_r, & ! The summed effective marginal face areas for the 3
    famt_0, &     ! test velocities [H m ~> m2 or kg m-1].
    vhtot_l, &    ! The summed transport with the southerly (vhtot_L) and
    vhtot_r       ! and northerly (vhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
  real :: FA_0    ! The effective face area with 0 barotropic transport [H L ~> m2 or kg m-1].
  real :: FA_avg  ! The average effective face area [H L ~> m2 or kg m-1], nominally given by
                  ! the realized transport divided by the barotropic velocity.
  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]  This
                       ! limiting is necessary to keep the inverse of visc_rem
                       ! from leading to large CFL numbers.
  real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
                       ! in finding the barotropic velocity that changes the
                       ! flow direction.  This is necessary to keep the inverse
                       ! of visc_rem from leading to large CFL numbers.
  real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
                  ! flow is truly upwind [nondim]
  real :: Idt     ! The inverse of the time step [T-1 ~> s-1].
  logical :: domore
  integer :: i, k, nz

  nz = g%ke ; idt = 1.0 / dt
  min_visc_rem = 0.1 ; cfl_min = 1e-6

 ! Diagnose the zero-transport correction, dv0.
  do i=ish,ieh ; zeros(i) = 0.0 ; enddo
  call meridional_flux_adjust(v, h_in, h_l, h_r, zeros, vh_tot_0, dvhdv_tot_0, dv0, &
                         dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
                         j, ish, ieh, do_i, .true.)

  !   Determine the southerly- and northerly- fluxes.  Choose a sufficiently
  ! negative velocity correction for the northerly-flux, and a sufficiently
  ! positive correction for the southerly-flux.
  domore = .false.
  do i=ish,ieh ; if (do_i(i)) then
    domore = .true.
    dv_cfl(i) = (cfl_min * idt) * g%dyCv(i,j)
    dvr(i) = min(0.0,dv0(i) - dv_cfl(i))
    dvl(i) = max(0.0,dv0(i) + dv_cfl(i))
    famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0
    vhtot_l(i) = 0.0 ; vhtot_r(i) = 0.0
  endif ; enddo

  if (.not.domore) then
    do k=1,nz ; do i=ish,ieh
      bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
      bt_cont%vBT_SS(i,j) = 0.0
      bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
      bt_cont%vBT_NN(i,j) = 0.0
    enddo ; enddo
    return
  endif

  do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then
    visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
    if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
      if (v(i,j,k) + dvr(i)*visc_rem_lim > -dv_cfl(i)*visc_rem(i,k)) &
        dvr(i) = -(v(i,j,k) + dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim
      if (v(i,j,k) + dvl(i)*visc_rem_lim < dv_cfl(i)*visc_rem(i,k)) &
        dvl(i) = -(v(i,j,k) - dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim
    endif
  endif ; enddo ; enddo
  do k=1,nz
    do i=ish,ieh ; if (do_i(i)) then
      v_l(i) = v(i,j,k) + dvl(i) * visc_rem(i,k)
      v_r(i) = v(i,j,k) + dvr(i) * visc_rem(i,k)
      v_0(i) = v(i,j,k) + dv0(i) * visc_rem(i,k)
    endif ; enddo
    call merid_flux_layer(v_0, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_0, dvhdv_0, &
                          visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
    call merid_flux_layer(v_l, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_l, dvhdv_l, &
                          visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
    call merid_flux_layer(v_r, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_r, dvhdv_r, &
                          visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
    do i=ish,ieh ; if (do_i(i)) then
      famt_0(i) = famt_0(i) + dvhdv_0(i)
      famt_l(i) = famt_l(i) + dvhdv_l(i)
      famt_r(i) = famt_r(i) + dvhdv_r(i)
      vhtot_l(i) = vhtot_l(i) + vh_l(i)
      vhtot_r(i) = vhtot_r(i) + vh_r(i)
    endif ; enddo
  enddo
  do i=ish,ieh ; if (do_i(i)) then
    fa_0 = famt_0(i) ; fa_avg = famt_0(i)
    if ((dvl(i) - dv0(i)) /= 0.0) &
      fa_avg = vhtot_l(i) / (dvl(i) - dv0(i))
    if (fa_avg > max(fa_0, famt_l(i))) then ; fa_avg = max(fa_0, famt_l(i))
    elseif (fa_avg < min(fa_0, famt_l(i))) then ; fa_0 = fa_avg ; endif
    bt_cont%FA_v_S0(i,j) = fa_0 ; bt_cont%FA_v_SS(i,j) = famt_l(i)
    if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0) then ; bt_cont%vBT_SS(i,j) = 0.0 ; else
      bt_cont%vBT_SS(i,j) = (1.5 * (dvl(i) - dv0(i))) * &
                   ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))
    endif

    fa_0 = famt_0(i) ; fa_avg = famt_0(i)
    if ((dvr(i) - dv0(i)) /= 0.0) &
      fa_avg = vhtot_r(i) / (dvr(i) - dv0(i))
    if (fa_avg > max(fa_0, famt_r(i))) then ; fa_avg = max(fa_0, famt_r(i))
    elseif (fa_avg < min(fa_0, famt_r(i))) then ; fa_0 = fa_avg ; endif
    bt_cont%FA_v_N0(i,j) = fa_0 ; bt_cont%FA_v_NN(i,j) = famt_r(i)
    if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0) then ; bt_cont%vBT_NN(i,j) = 0.0 ; else
      bt_cont%vBT_NN(i,j) = (1.5 * (dvr(i) - dv0(i))) * &
                   ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
    endif
  else
    bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
    bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
    bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
  endif ; enddo

end subroutine set_merid_bt_cont

!> Calculates left/right edge values for PPM reconstruction.
subroutine ppm_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd, OBC)
  type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_L  !< Left thickness in the reconstruction,
                                                         !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_R  !< Right thickness in the reconstruction,
                                                         !! [H ~> m or kg m-2].
  type(loop_bounds_type),            intent(in)  :: LB   !< Active loop bounds structure.
  real,                              intent(in)  :: h_min !< The minimum thickness
                    !! that can be obtained by a concave parabolic fit.
  logical, optional,                 intent(in)  :: monotonic !< If true, use the
                    !! Colella & Woodward monotonic limiter.
                    !! Otherwise use a simple positive-definite limiter.
  logical, optional,                 intent(in)  :: simple_2nd !< If true, use the
                    !! arithmetic mean thicknesses as the default edge values
                    !! for a simple 2nd order scheme.
  type(ocean_obc_type),    optional, pointer     :: OBC !< Open boundaries control structure.

  ! Local variables with useful mnemonic names.
  real, dimension(SZI_(G),SZJ_(G))  :: slp ! The slopes.
  real, parameter :: oneSixth = 1./6.
  real :: h_ip1, h_im1
  real :: dMx, dMn
  logical :: use_CW84, use_2nd
  character(len=256) :: mesg
  integer :: i, j, isl, iel, jsl, jel, n, stencil
  logical :: local_open_BC
  type(obc_segment_type), pointer :: segment => null()

  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd

  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif

  isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh

  ! This is the stencil of the reconstruction, not the scheme overall.
  stencil = 2 ; if (use_2nd) stencil = 1

  if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
    write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
               & "x-halo that needs to be increased by ",i2,".")') &
               stencil + max(g%isd-isl,iel-g%ied)
    call mom_error(fatal,mesg)
  endif
  if ((jsl < g%jsd) .or. (jel > g%jed)) then
    write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
               & "y-halo that needs to be increased by ",i2,".")') &
               max(g%jsd-jsl,jel-g%jed)
    call mom_error(fatal,mesg)
  endif

  if (use_2nd) then
    do j=jsl,jel ; do i=isl,iel
      h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
      h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
      h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
      h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
    enddo ; enddo
  else
    do j=jsl,jel ; do i=isl-1,iel+1
      if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
        slp(i,j) = 0.0
      else
        ! This uses a simple 2nd order slope.
        slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
        ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
        dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
        dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
        slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
                ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
      endif
    enddo ; enddo

    if (local_open_bc) then
      do n=1, obc%number_of_segments
        segment => obc%segment(n)
        if (.not. segment%on_pe) cycle
        if (segment%direction == obc_direction_e .or. &
            segment%direction == obc_direction_w) then
          i=segment%HI%IsdB
          do j=segment%HI%jsd,segment%HI%jed
            slp(i+1,j) = 0.0
            slp(i,j) = 0.0
          enddo
        endif
      enddo
    endif

    do j=jsl,jel ; do i=isl,iel
      ! Neighboring values should take into account any boundaries.  The 3
      ! following sets of expressions are equivalent.
    ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
    ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
      h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
      h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
      ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
      h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
      h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
    enddo ; enddo
  endif

  if (local_open_bc) then
    do n=1, obc%number_of_segments
      segment => obc%segment(n)
      if (.not. segment%on_pe) cycle
      if (segment%direction == obc_direction_e) then
        i=segment%HI%IsdB
        do j=segment%HI%jsd,segment%HI%jed
          h_l(i+1,j) = h_in(i,j)
          h_r(i+1,j) = h_in(i,j)
          h_l(i,j) = h_in(i,j)
          h_r(i,j) = h_in(i,j)
        enddo
      elseif (segment%direction == obc_direction_w) then
        i=segment%HI%IsdB
        do j=segment%HI%jsd,segment%HI%jed
          h_l(i,j) = h_in(i+1,j)
          h_r(i,j) = h_in(i+1,j)
          h_l(i+1,j) = h_in(i+1,j)
          h_r(i+1,j) = h_in(i+1,j)
        enddo
      endif
    enddo
  endif

  if (use_cw84) then
    call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
  else
    call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
  endif

  return
end subroutine ppm_reconstruction_x

!> Calculates left/right edge values for PPM reconstruction.
subroutine ppm_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd, OBC)
  type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_L  !< Left thickness in the reconstruction,
                                                         !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_R  !< Right thickness in the reconstruction,
                                                         !! [H ~> m or kg m-2].
  type(loop_bounds_type),            intent(in)  :: LB   !< Active loop bounds structure.
  real,                              intent(in)  :: h_min !< The minimum thickness
                    !! that can be obtained by a concave parabolic fit.
  logical, optional,                 intent(in)  :: monotonic !< If true, use the
                    !! Colella & Woodward monotonic limiter.
                    !! Otherwise use a simple positive-definite limiter.
  logical, optional,                 intent(in)  :: simple_2nd !< If true, use the
                    !! arithmetic mean thicknesses as the default edge values
                    !! for a simple 2nd order scheme.
  type(ocean_obc_type),    optional, pointer     :: OBC !< Open boundaries control structure.

  ! Local variables with useful mnemonic names.
  real, dimension(SZI_(G),SZJ_(G))  :: slp ! The slopes.
  real, parameter :: oneSixth = 1./6.
  real :: h_jp1, h_jm1
  real :: dMx, dMn
  logical :: use_CW84, use_2nd
  character(len=256) :: mesg
  integer :: i, j, isl, iel, jsl, jel, n, stencil
  logical :: local_open_BC
  type(obc_segment_type), pointer :: segment => null()

  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd

  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_v_BCs_exist_globally
  endif ; endif

  isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1

  ! This is the stencil of the reconstruction, not the scheme overall.
  stencil = 2 ; if (use_2nd) stencil = 1

  if ((isl < g%isd) .or. (iel > g%ied)) then
    write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
               & "x-halo that needs to be increased by ",i2,".")') &
               max(g%isd-isl,iel-g%ied)
    call mom_error(fatal,mesg)
  endif
  if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
    write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
                 & "y-halo that needs to be increased by ",i2,".")') &
                 stencil + max(g%jsd-jsl,jel-g%jed)
    call mom_error(fatal,mesg)
  endif

  if (use_2nd) then
    do j=jsl,jel ; do i=isl,iel
      h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
      h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
      h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
      h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
    enddo ; enddo
  else
    do j=jsl-1,jel+1 ; do i=isl,iel
      if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
        slp(i,j) = 0.0
      else
        ! This uses a simple 2nd order slope.
        slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
        ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
        dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
        dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
        slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
                ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
      endif
    enddo ; enddo

    if (local_open_bc) then
      do n=1, obc%number_of_segments
        segment => obc%segment(n)
        if (.not. segment%on_pe) cycle
        if (segment%direction == obc_direction_s .or. &
            segment%direction == obc_direction_n) then
          j=segment%HI%JsdB
          do i=segment%HI%isd,segment%HI%ied
            slp(i,j+1) = 0.0
            slp(i,j) = 0.0
          enddo
        endif
      enddo
    endif

    do j=jsl,jel ; do i=isl,iel
      ! Neighboring values should take into account any boundaries.  The 3
      ! following sets of expressions are equivalent.
      h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
      h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
      ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
      h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
      h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
    enddo ; enddo
  endif

  if (local_open_bc) then
    do n=1, obc%number_of_segments
      segment => obc%segment(n)
      if (.not. segment%on_pe) cycle
      if (segment%direction == obc_direction_n) then
        j=segment%HI%JsdB
        do i=segment%HI%isd,segment%HI%ied
          h_l(i,j+1) = h_in(i,j)
          h_r(i,j+1) = h_in(i,j)
          h_l(i,j) = h_in(i,j)
          h_r(i,j) = h_in(i,j)
        enddo
      elseif (segment%direction == obc_direction_s) then
        j=segment%HI%JsdB
        do i=segment%HI%isd,segment%HI%ied
          h_l(i,j) = h_in(i,j+1)
          h_r(i,j) = h_in(i,j+1)
          h_l(i,j+1) = h_in(i,j+1)
          h_r(i,j+1) = h_in(i,j+1)
        enddo
      endif
    enddo
  endif

  if (use_cw84) then
    call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
  else
    call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
  endif

  return
end subroutine ppm_reconstruction_y

!> This subroutine limits the left/right edge values of the PPM reconstruction
!! to give a reconstruction that is positive-definite.  Here this is
!! reinterpreted as giving a constant thickness if the mean thickness is less
!! than h_min, with a minimum of h_min otherwise.
subroutine ppm_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
  type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].
  real,                              intent(in)  :: h_min !< The minimum thickness
                    !! that can be obtained by a concave parabolic fit.
  integer,                           intent(in)  :: iis      !< Start of i index range.
  integer,                           intent(in)  :: iie      !< End of i index range.
  integer,                           intent(in)  :: jis      !< Start of j index range.
  integer,                           intent(in)  :: jie      !< End of j index range.

! Local variables
  real    :: curv, dh, scale
  character(len=256) :: mesg
  integer :: i,j

  do j=jis,jie ; do i=iis,iie
    ! This limiter prevents undershooting minima within the domain with
    ! values less than h_min.
    curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
    if (curv > 0.0) then ! Only minima are limited.
      dh = h_r(i,j) - h_l(i,j)
      if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
        if (h_in(i,j) <= h_min) then
          h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
        elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
          ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
          ! be limited in this case.  0 < scale < 1.
          scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
          h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
          h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
        endif
      endif
    endif
  enddo ; enddo

end subroutine ppm_limit_pos

!> This subroutine limits the left/right edge values of the PPM reconstruction
!! according to the monotonic prescription of Colella and Woodward, 1984.
subroutine ppm_limit_cw84(h_in, h_L, h_R, G, iis, iie, jis, jie)
  type(ocean_grid_type),             intent(in)  :: G     !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in  !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_L !< Left thickness in the reconstruction,
                                                          !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_R !< Right thickness in the reconstruction,
                                                          !! [H ~> m or kg m-2].
  integer,                           intent(in)  :: iis   !< Start of i index range.
  integer,                           intent(in)  :: iie   !< End of i index range.
  integer,                           intent(in)  :: jis   !< Start of j index range.
  integer,                           intent(in)  :: jie   !< End of j index range.

  ! Local variables
  real    :: h_i, RLdiff, RLdiff2, RLmean, FunFac
  character(len=256) :: mesg
  integer :: i,j

  do j=jis,jie ; do i=iis,iie
    ! This limiter monotonizes the parabola following
    ! Colella and Woodward, 1984, Eq. 1.10
    h_i = h_in(i,j)
    if ( ( h_r(i,j) - h_i ) * ( h_i - h_l(i,j) ) <= 0. ) then
      h_l(i,j) = h_i ; h_r(i,j) = h_i
    else
      rldiff = h_r(i,j) - h_l(i,j)            ! Difference of edge values
      rlmean = 0.5 * ( h_r(i,j) + h_l(i,j) )  ! Mean of edge values
      funfac = 6. * rldiff * ( h_i - rlmean ) ! Some funny factor
      rldiff2 = rldiff * rldiff               ! Square of difference
      if ( funfac >  rldiff2 ) h_l(i,j) = 3. * h_i - 2. * h_r(i,j)
      if ( funfac < -rldiff2 ) h_r(i,j) = 3. * h_i - 2. * h_l(i,j)
    endif
  enddo ; enddo

  return
end subroutine ppm_limit_cw84

!> Return the maximum ratio of a/b or maxrat.
function ratio_max(a, b, maxrat) result(ratio)
  real, intent(in) :: a       !< Numerator
  real, intent(in) :: b       !< Denominator
  real, intent(in) :: maxrat  !< Maximum value of ratio.
  real :: ratio               !< Return value.

  if (abs(a) > abs(maxrat*b)) then
    ratio = maxrat
  else
    ratio = a / b
  endif
end function ratio_max

!> Initializes continuity_ppm_cs
subroutine continuity_ppm_init(Time, G, GV, US, param_file, diag, CS)
  type(time_type), target, intent(in)    :: time !< The current model time.
  type(ocean_grid_type),   intent(in)    :: g    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: gv   !< Vertical grid structure.
  type(unit_scale_type),   intent(in)    :: us  !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< A structure indicating
                  !! the open file to parse for model parameter values.
  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to
                  !! regulate diagnostic output.
  type(continuity_ppm_cs), pointer       :: cs   !< Module's control structure.
!> This include declares and sets the variable "version".
#include "version_variable.h"
  real :: tol_eta_m  ! An unscaled version of tol_eta [m].
  character(len=40)  :: mdl = "MOM_continuity_PPM" ! This module's name.

  if (associated(cs)) then
    call mom_error(warning, "continuity_PPM_init called with associated control structure.")
    return
  endif
  allocate(cs)

! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", cs%monotonic, &
                 "If true, CONTINUITY_PPM uses the Colella and Woodward "//&
                 "monotonic limiter.  The default (false) is to use a "//&
                 "simple positive definite limiter.", default=.false.)
  call get_param(param_file, mdl, "SIMPLE_2ND_PPM_CONTINUITY", cs%simple_2nd, &
                 "If true, CONTINUITY_PPM uses a simple 2nd order "//&
                 "(arithmetic mean) interpolation of the edge values. "//&
                 "This may give better PV conservation properties. While "//&
                 "it formally reduces the accuracy of the continuity "//&
                 "solver itself in the strongly advective limit, it does "//&
                 "not reduce the overall order of accuracy of the dynamic "//&
                 "core.", default=.false.)
  call get_param(param_file, mdl, "UPWIND_1ST_CONTINUITY", cs%upwind_1st, &
                 "If true, CONTINUITY_PPM becomes a 1st-order upwind "//&
                 "continuity solver.  This scheme is highly diffusive "//&
                 "but may be useful for debugging or in single-column "//&
                 "mode where its minimal stencil is useful.", default=.false.)
  call get_param(param_file, mdl, "ETA_TOLERANCE", cs%tol_eta, &
                 "The tolerance for the differences between the "//&
                 "barotropic and baroclinic estimates of the sea surface "//&
                 "height due to the fluxes through each face.  The total "//&
                 "tolerance for SSH is 4 times this value.  The default "//&
                 "is 0.5*NK*ANGSTROM, and this should not be set less "//&
                 "than about 10^-15*MAXIMUM_DEPTH.", units="m", scale=gv%m_to_H, &
                 default=0.5*g%ke*gv%Angstrom_m, unscaled=tol_eta_m)

  !### ETA_TOLERANCE_AUX can be obsoleted.
  call get_param(param_file, mdl, "ETA_TOLERANCE_AUX", cs%tol_eta_aux, &
                 "The tolerance for free-surface height discrepancies "//&
                 "between the barotropic solution and the sum of the "//&
                 "layer thicknesses when calculating the auxiliary "//&
                 "corrected velocities. By default, this is the same as "//&
                 "ETA_TOLERANCE, but can be made larger for efficiency.", &
                 units="m", default=tol_eta_m, scale=gv%m_to_H)
  call get_param(param_file, mdl, "VELOCITY_TOLERANCE", cs%tol_vel, &
                 "The tolerance for barotropic velocity discrepancies "//&
                 "between the barotropic solution and  the sum of the "//&
                 "layer thicknesses.", units="m s-1", default=3.0e8, scale=us%m_s_to_L_T)
                 ! The speed of light is the default.

  call get_param(param_file, mdl, "CONT_PPM_AGGRESS_ADJUST", cs%aggress_adjust,&
                 "If true, allow the adjusted velocities to have a "//&
                 "relative CFL change up to 0.5.", default=.false.)
  cs%vol_CFL = cs%aggress_adjust
  call get_param(param_file, mdl, "CONT_PPM_VOLUME_BASED_CFL", cs%vol_CFL, &
                 "If true, use the ratio of the open face lengths to the "//&
                 "tracer cell areas when estimating CFL numbers.  The "//&
                 "default is set by CONT_PPM_AGGRESS_ADJUST.", &
                 default=cs%aggress_adjust, do_not_read=cs%aggress_adjust)
  call get_param(param_file, mdl, "CONTINUITY_CFL_LIMIT", cs%CFL_limit_adjust, &
                 "The maximum CFL of the adjusted velocities.", units="nondim", &
                 default=0.5)
  call get_param(param_file, mdl, "CONT_PPM_BETTER_ITER", cs%better_iter, &
                 "If true, stop corrective iterations using a velocity "//&
                 "based criterion and only stop if the iteration is "//&
                 "better than all predecessors.", default=.true.)
  call get_param(param_file, mdl, "CONT_PPM_USE_VISC_REM_MAX", &
                                 cs%use_visc_rem_max, &
                 "If true, use more appropriate limiting bounds for "//&
                 "corrections in strongly viscous columns.", default=.true.)
  call get_param(param_file, mdl, "CONT_PPM_MARGINAL_FACE_AREAS", cs%marginal_faces, &
                 "If true, use the marginal face areas from the continuity "//&
                 "solver for use as the weights in the barotropic solver. "//&
                 "Otherwise use the transport averaged areas.", default=.true.)

  cs%diag => diag

  id_clock_update = cpu_clock_id('(Ocean continuity update)', grain=clock_routine)
  id_clock_correct = cpu_clock_id('(Ocean continuity correction)', grain=clock_routine)

end subroutine continuity_ppm_init

!> continuity_PPM_stencil returns the continuity solver stencil size
function continuity_ppm_stencil(CS) result(stencil)
  type(continuity_ppm_cs), pointer       :: cs  !< Module's control structure.
  integer ::  stencil !< The continuity solver stencil size with the current settings.

  stencil = 3 ; if (cs%simple_2nd) stencil = 2 ; if (cs%upwind_1st) stencil = 1

end function continuity_ppm_stencil

!> Destructor for continuity_ppm_cs
subroutine continuity_ppm_end(CS)
  type(continuity_ppm_cs), pointer :: cs   !< Module's control structure.
  deallocate(cs)
end subroutine continuity_ppm_end

!> \namespace mom_continuity_ppm
!!
!! This module contains the subroutines that advect layer
!! thickness.  The scheme here uses a Piecewise-Parabolic method with
!! a positive definite limiter.

end module mom_continuity_ppm