MOM_continuity_PPM.F90

Go to the documentation of this file.

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

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, in m.
  real :: tol_vel            !< The tolerance for barotropic velocity
                             !! discrepancies between the barotropic solution and
                             !! the sum of the layer thicknesses, in 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, in m.
  real :: cfl_limit_adjust   !< The maximum CFL of the adjusted velocities, ND.
  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, CS, uhbt, vhbt, OBC, &
                          visc_rem_u, visc_rem_v, u_cor, v_cor, &
                          uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont)
  ! In the following documentation, H is used for the units of thickness (usually m or kg m-2.)
  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, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v   !< Meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: hin !< Initial layer thickness, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h   !< Final layer thickness, in H.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out)   :: uh  !< Zonal volume flux,
                                                                  !! u*h*dy, H m2 s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: vh  !< Meridional volume flux,
                                                                  !! v*h*dx, H m2 s-1.
  real,                                      intent(in)    :: dt  !< Time increment in s.
  type(verticalgrid_type),                   intent(in)    :: GV  !< Vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G)),         intent(in),  optional :: uhbt
         !< The summed volume flux through zonal faces, H m2 s-1.
  real, dimension(SZI_(G),SZJB_(G)),         intent(in),  optional :: vhbt
         !< The summed volume flux through meridional faces, H m2 s-1.
  type(ocean_obc_type),                      pointer,     optional :: OBC !< Open boundaries control structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in),  optional :: 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)), intent(in),  optional :: 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)), intent(out), optional :: u_cor
         !< The zonal velocities that give uhbt as the depth-integrated transport, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor
         !< The meridional velocities that give vhbt as the depth-integrated transport, in m s-1.
  real, dimension(SZIB_(G),SZJ_(G)),         intent(in),  optional :: uhbt_aux
         !< A second set of summed volume fluxes through zonal faces, in H m2 s-1.
  real, dimension(SZI_(G),SZJB_(G)),         intent(in),  optional :: vhbt_aux
         !< A second set of summed volume fluxes through meridional faces, in H m2 s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux
         !< The zonal velocities that give uhbt_aux as the depth-integrated transports, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor_aux
         !< The meridional velocities that give vhbt_aux as the depth-integrated transports, in m s-1.
  type(bt_cont_type),                        pointer,     optional :: BT_cont !< A structure with
         !! elements that describe the effective open face areas as a function of barotropic flow.
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_input ! Left and right face thicknesses, in H.
  real :: h_min
  type(loop_bounds_type) :: LB
  integer :: is, ie, js, je, nz, stencil
  integer :: i, j, k

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

  h_min = gv%Angstrom

  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)

  local_flather_ew = .false. ; local_flather_ns = .false.
  if (present(obc)) then ; if (associated(obc)) then ; if (obc%OBC_pe) then
    local_flather_ew = obc%Flather_u_BCs_exist_globally
    local_flather_ns = obc%Flather_v_BCs_exist_globally
    !   If an OBC is being applied, copy the input thicknesses so that the
    ! OBC code works even if hin == h.
    if (local_flather_ew .or.  local_flather_ns) &
      h_input(:,:,:) = hin(:,:,:)
  endif ; endif ; endif

  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, cs, lb, uhbt, obc, visc_rem_u, &
                         u_cor, uhbt_aux, u_cor_aux, bt_cont)

    call cpu_clock_begin(id_clock_update)
!$OMP parallel do default(none) shared(LB,nz,G,uh,hin,dt,h)
    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)

    if (local_flather_ew) then
      do k=1,nz
        do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh+1
          if (obc%segnum_u(i-1,j) /= obc_none) then
            if (obc%segment(obc%segnum_u(i-1,j))%direction == obc_direction_e) &
              h(i,j,k) = h_input(i-1,j,k)
          endif
        enddo
        do i=lb%ish-1,lb%ieh
          if (obc%segnum_u(i,j) /= obc_none) then
            if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) &
              h(i,j,k) = h_input(i+1,j,k)
          endif
        enddo ; enddo
      enddo
    endif
    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, cs, lb, vhbt, obc, visc_rem_v, &
                              v_cor, vhbt_aux, v_cor_aux, bt_cont)

    call cpu_clock_begin(id_clock_update)
!$OMP parallel do default(none) shared(h_min,nz,LB,h,dt,G,vh)
    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)

    if (local_flather_ns) then
      do k=1,nz
        do j=lb%jsh,lb%jeh+1 ; do i=lb%ish-1,lb%ieh+1
          if (obc%segnum_v(i,j-1) /= obc_none) then
            if (obc%segment(obc%segnum_v(i,j-1))%direction == obc_direction_n) &
              h(i,j,k) = h_input(i,j-1,k)
          endif
        enddo ; enddo
        do j=lb%jsh-1,lb%jeh ; do i=lb%ish-1,lb%ieh+1
          if (obc%segnum_v(i,j) /= obc_none) then
            if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) &
              h(i,j,k) = h_input(i,j+1,k)
          endif
        enddo ; enddo
      enddo
    endif
  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, cs, lb, vhbt, obc, visc_rem_v, &
                              v_cor, vhbt_aux, v_cor_aux, bt_cont)

    call cpu_clock_begin(id_clock_update)
!$OMP parallel do default(none) shared(nz,LB,h,hin,dt,G,vh)
    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)

    if (local_flather_ns) then
      do k=1,nz
        do j=lb%jsh,lb%jeh+1 ; do i=lb%ish-1,lb%ieh+1
          if (obc%segnum_v(i,j-1) /= obc_none) then
            if (obc%segment(obc%segnum_v(i,j-1))%direction == obc_direction_n) &
              h(i,j,k) = h_input(i,j-1,k)
          endif
        enddo ; enddo
        do j=lb%jsh-1,lb%jeh ; do i=lb%ish-1,lb%ieh+1
          if (obc%segnum_v(i,j) /= obc_none) then
            if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) &
              h(i,j,k) = h_input(i,j+1,k)
          endif
        enddo ; enddo
      enddo
    endif

  !    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, cs, lb, uhbt, obc, visc_rem_u, &
                         u_cor, uhbt_aux, u_cor_aux, bt_cont)

    call cpu_clock_begin(id_clock_update)
!$OMP parallel do default(none) shared(h_min,nz,LB,h,dt,G,uh)
    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)

    if (local_flather_ew) then
      do k=1,nz
        do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh+1
          if (obc%segnum_u(i-1,j) /= obc_none) then
            if (obc%segment(obc%segnum_u(i-1,j))%direction == obc_direction_e) &
              h(i,j,k) = h_input(i-1,j,k)
          endif
        enddo
        do i=lb%ish-1,lb%ieh
          if (obc%segnum_u(i,j) /= obc_none) then
            if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) &
              h(i,j,k) = h_input(i+1,j,k)
          endif
        enddo ; enddo
      enddo
    endif
  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, CS, LB, uhbt, OBC, &
                           visc_rem_u, u_cor, uhbt_aux, u_cor_aux, 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, in m s-1.
  real,  dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes, in H.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out)   :: uh   !< Volume flux through zonal
                                                                   !! faces = u*h*dy, H m2 s-1.
  real,                                      intent(in)    :: dt   !< Time increment in s.
  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),                      pointer,     optional :: OBC !< Open boundaries control structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in),  optional :: 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)),         intent(in),  optional :: uhbt
         !< The summed volume flux through zonal faces, H m2 s-1.
  real, dimension(SZIB_(G),SZJ_(G)),         intent(in),  optional :: uhbt_aux
         !< A second set of summed volume fluxes through zonal faces, in H m2 s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor
         !< The zonal velocitiess (u with a barotropic correction)
         !! that give uhbt as the depth-integrated transport, m s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux
         !< The zonal velocities (u with a barotropic correction)
         !! that give uhbt_aux as the depth-integrated transports, in m s-1.
  type(bt_cont_type),                        pointer,     optional :: BT_cont !<
         !< A structure with elements that describe the effective
         !! open face areas as a function of barotropic flow.
  ! Local variables
  real, dimension(SZIB_(G),SZK_(G)) :: duhdu ! Partial derivative of uh with u, in H m.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_L, h_R ! Left and right face thicknesses, in H.
  real, dimension(SZIB_(G)) :: &
    du, &      ! Corrective barotropic change in the velocity, in 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, in H m.
    uh_tot_0, & ! Summed transport with no barotropic correction in H m2 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 :: I_vrm   ! 1.0 / visc_rem_max, nondim.
  real :: CFL_dt  ! The maximum CFL ratio of the adjusted velocities divided by
                  ! the time step, in s-1.
  real :: I_dt    ! 1.0 / dt, in s-1.
  real :: du_lim  ! The velocity change that give a relative CFL of 1, in m s-1.
  real :: dx_E, dx_W ! Effective x-grid spacings to the east and west, in m.
  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
  logical :: do_aux, 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

  do_aux = (present(uhbt_aux) .and. present(u_cor_aux))
  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 .or. &
                        obc%Flather_v_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)
  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, cs%monotonic, simple_2nd=cs%simple_2nd)
    endif
    do i=ish-1,ieh ; visc_rem(i,k) = 1.0 ; enddo
  enddo
  if (local_open_bc) then
    do n=1, obc%number_of_segments
      segment => obc%segment(n)
      if (.not. segment%on_pe .or. segment%specified) cycle
      if (segment%direction == obc_direction_e) then
        i=segment%HI%IsdB
        do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
          h_in(i+1,j,k) = h_in(i,j,k)
          h_l(i+1,j,k) = h_in(i,j,k)
          h_r(i+1,j,k) = h_in(i,j,k)
          h_r(i,j,k) = h_in(i,j,k)
        enddo ; enddo
      elseif (segment%direction == obc_direction_w) then
        i=segment%HI%IsdB
        do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
            h_in(i,j,k) = h_in(i+1,j,k)
            h_l(i,j,k) = h_in(i+1,j,k)
            h_r(i,j,k) = h_in(i+1,j,k)
            h_l(i+1,j,k) = h_in(i+1,j,k)
        enddo ; enddo
      endif
    enddo
  endif
  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,G,GV,CS,local_specified_BC,OBC,uhbt,do_aux,set_BT_cont,    &
!$OMP                                  CFL_dt,I_dt,u_cor,uhbt_aux,u_cor_aux,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, &
!$OMP                                  visc_rem_max, I_vrm, du_lim, dx_E, dx_W, any_simple_OBC ) &
!$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, j, ish, ieh, do_i, cs%vol_CFL)
      if (local_specified_bc) then
        do i=ish-1,ieh
          if (obc%segment(obc%segnum_u(i,j))%specified) &
            uh(i,j,k) = obc%segment(obc%segnum_u(i,j))%normal_trans(i,j,k)
        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. do_aux .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

      ! Up to this point, everything is shared between uhbt and uhbt_aux.

      any_simple_obc = .false.
      if (present(uhbt) .or. do_aux .or. set_bt_cont) then
        if (local_specified_bc .or. local_flather_obc) then ; do i=ish-1,ieh
          ! Avoid reconciling barotropic/baroclinic transports if transport is specified
          is_simple = obc%segment(obc%segnum_u(i,j))%specified
          do_i(i) = .not.(obc%segnum_u(i,j) /= 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, &
                               cs, visc_rem, j, ish, ieh, do_i, .true., uh)

        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
            if (obc%segment(obc%segnum_u(i,j))%specified) &
              u_cor(i,j,k) = obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
          enddo ; endif
        enddo ; endif ! u-corrected

      endif

      if (do_aux) then
        call zonal_flux_adjust(u, h_in, h_l, h_r, uhbt_aux(:,j), uh_tot_0, &
                               duhdu_tot_0, du, du_max_cfl, du_min_cfl, dt, g, &
                               cs, visc_rem, j, ish, ieh, do_i, .false.)

        do k=1,nz
          do i=ish-1,ieh ; u_cor_aux(i,j,k) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
          if (local_specified_bc) then ; do i=ish-1,ieh
            if (obc%segment(obc%segnum_u(i,j))%specified) &
              u_cor_aux(i,j,k) = obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
          enddo ; endif
        enddo
      endif ! do_aux

      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, cs, visc_rem, &
                               visc_rem_max, j, ish, ieh, do_i)
        if (any_simple_obc) then
          do i=ish-1,ieh
            do_i(i) = obc%segment(obc%segnum_u(i,j))%specified
            if (do_i(i)) bt_cont%Fa_u_W0(i,j) = gv%H_subroundoff*g%dy_Cu(i,j)
          enddo
          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) &
              bt_cont%Fa_u_W0(i,j) = bt_cont%Fa_u_W0(i,j) + &
                   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_E0(i,j) = bt_cont%Fa_u_W0(i,j)
            bt_cont%Fa_u_WW(i,j) = bt_cont%Fa_u_W0(i,j)
            bt_cont%Fa_u_EE(i,j) = bt_cont%Fa_u_W0(i,j)
          endif ; enddo
        endif
      endif ! set_BT_cont

    endif ! present(uhbt) or do_aux or set_BT_cont
  enddo ! j-loop
  call cpu_clock_end(id_clock_correct)

  if  (set_bt_cont) then ; if (associated(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, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_u)
    else
      call zonal_face_thickness(u, h_in, h_l, h_r, bt_cont%h_u, dt, g, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_u)
    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, j, &
                            ish, ieh, do_I, vol_CFL)
  type(ocean_grid_type),        intent(inout) :: G        !< Ocean's grid structure.
  real, dimension(SZIB_(G)),    intent(in)    :: u        !< Zonal velocity, in 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, in H.
  real, dimension(SZI_(G)),     intent(in)    :: h_L      !< Left thickness, in H.
  real, dimension(SZI_(G)),     intent(in)    :: h_R      !< Right thickness, in H.
  real, dimension(SZIB_(G)),    intent(inout) :: uh       !< Zonal mass or volume
                                                          !! transport, in H m2 s-1.
  real, dimension(SZIB_(G)),    intent(inout) :: duhdu    !< Partial derivative of uh
                                                          !! with u, in H m.
  real,                         intent(in)    :: dt       !< Time increment in s.
  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.
  ! Local variables
  real :: CFL  ! The CFL number based on the local velocity and grid spacing, ND.
  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, in H.
  integer :: i

  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

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, LB, vol_CFL, &
                                marginal, visc_rem_u)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< Zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thickness used to
                                                                   !! calculate fluxes, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction, in H.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_u  !< Thickness at zonal faces,
                                                                   !! in H.
  real,                                      intent(in)    :: dt   !< Time increment in s.
  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)), intent(in), optional :: visc_rem_u !<
                          !! Both the fraction of the momentum originally in a
                          !! layer that remains after a time-step of viscosity,
                          !! and the fraction of a time-step's worth of a
                          !! barotropic acceleration that a layer experiences
                          !! after viscosity is applied. Non-dimensional between
                          !! 0 (at the bottom) and 1 (far above the bottom).
  ! Local variables
  real :: CFL  ! The CFL number based on the local velocity and grid spacing, ND.
  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, in H.
  real :: h_marg ! The marginal thickness of a flux, in H.
  integer :: i, j, k, ish, ieh, jsh, jeh, nz
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke

!$OMP parallel default(none) shared(ish,ieh,jsh,jeh,nz,u,vol_CFL,dt,G, &
!$OMP                               h_L,h_R,h,h_u,visc_rem_u,marginal) &
!$OMP                       private(CFL,curv_3,h_marg,h_avg)
!$OMP do
  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 do
    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
!$OMP end parallel

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, CS, visc_rem, &
                             j, ish, ieh, do_I_in, full_precision, uh_3d)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< Zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction, in H.
  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),  optional :: uhbt !<
                       !! The summed volume flux through zonal faces, H m2 s-1.
  real, dimension(SZIB_(G)),                 intent(in)    :: du_max_CFL  !< Maximum acceptable
                       !! value of du, in m s-1.
  real, dimension(SZIB_(G)),                 intent(in)    :: du_min_CFL  !< Minimum acceptable
                       !! value of du, in m s-1.
  real, dimension(SZIB_(G)),                 intent(in)    :: uh_tot_0    !<
                       !! The summed transport with 0 adjustment, in H m2 s-1.
  real, dimension(SZIB_(G)),                 intent(in)    :: duhdu_tot_0 !<
                       !! The partial derivative of du_err with du at 0 adjustment, in H m.
  real, dimension(SZIB_(G)),                 intent(out)   :: du !<
                       !! The barotropic velocity adjustment, in m s-1.
  real,                                      intent(in)    :: dt   !< Time increment in s.
  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,                                   intent(in), optional :: full_precision !<
                       !! A flag indicating how carefully to iterate.  The
                       !! default is .true. (more accurate).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout), optional :: uh_3d !<
                       !! Volume flux through zonal faces = u*h*dy, H m2 s-1.
  ! Local variables
  real, dimension(SZIB_(G),SZK_(G)) :: &
    uh_aux, &  ! An auxiliary zonal volume flux, in H m s-1.
    duhdu      ! Partial derivative of uh with u, in H m.
  real, dimension(SZIB_(G)) :: &
    uh_err, &  ! Difference between uhbt and the summed uh, in H m2 s-1.
    uh_err_best, & ! The smallest value of uh_err found so far, in H m2 s-1.
    u_new, &   ! The velocity with the correction added, in m s-1.
    duhdu_tot,&! Summed partial derivative of uh with u, in H m.
    du_min, &  ! Min/max limits on du correction based on CFL limits
    du_max     ! and previous iterations, in m s-1.
  real :: du_prev ! The previous value of du, in m s-1.
  real :: ddu    ! The change in du from the previous iteration, in m s-1.
  real :: tol_eta ! The tolerance for the current iteration, in m.
  real :: tol_vel ! The tolerance for velocity in the current iteration, 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, j, ish, ieh, do_i, cs%vol_CFL)
    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, 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, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction, in H.
  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, in H m2 s-1.
  real, dimension(SZIB_(G)),                 intent(in)    :: duhdu_tot_0 !<
                       !! The partial derivative of du_err with du at 0 adjustment, in H m.
  real, dimension(SZIB_(G)),                 intent(in)    :: du_max_CFL  !< Maximum acceptable
                       !! value of du, in m s-1.
  real, dimension(SZIB_(G)),                 intent(in)    :: du_min_CFL  !< Minimum acceptable
                       !! value of du, in m s-1.
  real,                                      intent(in)    :: dt   !< Time increment in s.
  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, m s-1.
    duL, duR, &   ! The barotropic velocity increments that give the westerly
                  ! (duL) and easterly (duR) test velocities.
    zeros, &      ! An array of full of 0's.
    du_cfl, &     ! The velocity increment that corresponds to CFL_min, in 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, in m s-1.
    fa_marg_l, &  ! The effective layer marginal face areas with the westerly
    fa_marg_r, &  ! (_L), easterly (_R), and zero-barotropic (_0) test
    fa_marg_0, &  ! velocities, in H m.
    uh_l, uh_r, & ! The layer transports with the westerly (_L), easterly (_R),
    uh_0, &       ! and zero-barotropic (_0) test velocities, in H m2 s-1.
    famt_l, famt_r, & ! The summed effective marginal face areas for the 3
    famt_0, &     ! test velocities, in H m.
    uhtot_l, &    ! The summed transport with the westerly (uhtot_L) and
    uhtot_r       ! and easterly (uhtot_R) test velocities, in H m2 s-1.
  real :: FA_0    ! The effective face area with 0 barotropic transport, in m H.
  real :: FA_avg  ! The average effective face area, in m H, nominally given by
                  ! the realized transport divided by the barotropic velocity.
  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem, ND.  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, ND.
  real :: Idt     ! The inverse of the time step, in 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, &
                         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%uBT_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_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, &
                          fa_marg_0, visc_rem(:,k), dt, g, 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, &
                          fa_marg_l, visc_rem(:,k), dt, g, 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, &
                          fa_marg_r, visc_rem(:,k), dt, g, j, ish, ieh, do_i, &
                          cs%vol_CFL)
    do i=ish-1,ieh ; if (do_i(i)) then
      famt_0(i) = famt_0(i) + fa_marg_0(i)
      famt_l(i) = famt_l(i) + fa_marg_l(i)
      famt_r(i) = famt_r(i) + fa_marg_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%uBT_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_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, CS, LB, vhbt, OBC, &
                                visc_rem_v, v_cor, vhbt_aux, v_cor_aux, 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, in m s-1.
  real,  dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes, in H.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: vh   !< Volume flux through meridional
                                                                   !! faces = v*h*dx, H m2 s-1.
  real,                                      intent(in)    :: dt   !< Time increment in s.
  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),                      pointer,     optional :: OBC !<
         !! This open boundary condition type specifies whether, where,
         !! and what open boundary conditions are used.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in),  optional :: visc_rem_v !<
         !! Both the fraction of the momentum originally in a
         !! layer that remains after a time-step of viscosity,
         !! and the fraction of a time-step's worth of a
         !! barotropic acceleration that a layer experiences
         !! after viscosity is applied.  Nondimensional between
         !! 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJB_(G)),         intent(in),  optional :: vhbt !<
         !! The summed volume flux through meridional faces, H m2 s-1.
  real, dimension(SZI_(G),SZJB_(G)),         intent(in),  optional :: vhbt_aux !<
         !! A second set of summed volume fluxes through meridional
         !! faces, in H m2 s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor !<
         !! The meridional velocitiess (v with a barotropic correction)
         !! that give vhbt as the depth-integrated transport, m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor_aux !<
         !! The meridional velocities (v with a barotropic correction)
         !! that give vhbt_aux as the depth-integrated transports, in m s-1.
  type(bt_cont_type),                        pointer,     optional :: BT_cont !<
         !! A structure with elements that describe the effective
  ! Local variables
  real, dimension(SZI_(G),SZK_(G)) :: &
    dvhdv      ! Partial derivative of vh with v, in m2.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    h_L, h_R   ! Left and right face thicknesses, in m.
  real, dimension(SZI_(G)) :: &
    dv, &      ! Corrective barotropic change in the velocity, in 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, in H m.
    vh_tot_0, &   ! Summed transport with no barotropic correction in H m2 s-1.
    visc_rem_max  ! The column maximum of visc_rem.
  logical, dimension(SZI_(G)) :: do_I
  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, in s-1.
  real :: I_dt    ! 1.0 / dt, in s-1.
  real :: dv_lim  ! The velocity change that give a relative CFL of 1, in m s-1.
  real :: dy_N, dy_S ! Effective y-grid spacings to the north and south, in m.
  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
  logical :: do_aux, 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

  do_aux = (present(vhbt_aux) .and. present(v_cor_aux))
  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_u_BCs_exist_globally .or. &
                        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)
  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, cs%monotonic, simple_2nd=cs%simple_2nd)
    endif
    do i=ish,ieh ; visc_rem(i,k) = 1.0 ; enddo
  enddo
  if (local_open_bc) then
    do n=1, obc%number_of_segments
      segment => obc%segment(n)
      if (.not. segment%on_pe .or. segment%specified) cycle
      if (segment%direction == obc_direction_n) then
        j=segment%HI%JsdB
        do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
          h_in(i,j+1,k) = h_in(i,j,k)
          h_l(i,j+1,k) = h_in(i,j,k)
          h_r(i,j+1,k) = h_in(i,j,k)
          h_r(i,j,k) = h_in(i,j,k)
        enddo ; enddo
      elseif (segment%direction == obc_direction_s) then
        j=segment%HI%JsdB
        do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
          h_in(i,j,k) = h_in(i,j+1,k)
          h_l(i,j,k) = h_in(i,j+1,k)
          h_r(i,j,k) = h_in(i,j+1,k)
          h_l(i,j+1,k) = h_in(i,j+1,k)
        enddo ; enddo
      endif
    enddo
  endif
  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,G,GV,CS,local_specified_BC,OBC,vhbt,do_aux, &
!$OMP                                  set_BT_cont,CFL_dt,I_dt,v_cor,vhbt_aux,          &
!$OMP                                  v_cor_aux,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, &
!$OMP                                  dy_S,any_simple_OBC ) &
!$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, j, ish, ieh, do_i, cs%vol_CFL)
      if (local_specified_bc) then
        do i=ish,ieh
          if (obc%segment(obc%segnum_v(i,j))%specified) &
            vh(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_trans(i,j,k)
        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. do_aux .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

      ! Up to this point, everything is shared between vhbt and vhbt_aux.

      any_simple_obc = .false.
      if (present(vhbt) .or. do_aux .or. set_bt_cont) then
        if (local_specified_bc .or. local_flather_obc) then ; do i=ish,ieh
          ! Avoid reconciling barotropic/baroclinic transports if transport is specified
          is_simple = obc%segment(obc%segnum_v(i,j))%specified
          do_i(i) = .not.(obc%segnum_v(i,j) /= 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, &
                               cs, visc_rem, j, ish, ieh, do_i, .true., vh)

        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
            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)
          enddo ; endif
        enddo ; endif ! v-corrected
      endif

      if (do_aux) then
        call meridional_flux_adjust(v, h_in, h_l, h_r, vhbt_aux(:,j), vh_tot_0, &
                               dvhdv_tot_0, dv, dv_max_cfl, dv_min_cfl, dt, g, &
                               cs, visc_rem, j, ish, ieh, do_i, .false.)

        do k=1,nz
          do i=ish,ieh ; v_cor_aux(i,j,k) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
          if (local_specified_bc) then ; do i=ish,ieh
            if (obc%segment(obc%segnum_v(i,j))%specified) &
              v_cor_aux(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
          enddo ; endif
        enddo
      endif ! do_aux

      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, cs, visc_rem, &
                               visc_rem_max, j, ish, ieh, do_i)
        if (any_simple_obc) then
          do i=ish,ieh
            do_i(i) = (obc%segment(obc%segnum_v(i,j))%specified)
            if (do_i(i)) bt_cont%Fa_v_S0(i,j) = gv%H_subroundoff*g%dx_Cv(i,j)
          enddo
          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) &
              bt_cont%Fa_v_S0(i,j) = bt_cont%Fa_v_S0(i,j) + &
                   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_N0(i,j) = bt_cont%Fa_v_S0(i,j)
            bt_cont%Fa_v_SS(i,j) = bt_cont%Fa_v_S0(i,j)
            bt_cont%Fa_v_NN(i,j) = bt_cont%Fa_v_S0(i,j)
          endif ; enddo
        endif
      endif ! set_BT_cont

    endif ! present(vhbt) or do_aux or set_BT_cont
  enddo ! j-loop
  call cpu_clock_end(id_clock_correct)

  if (set_bt_cont) then ; if (associated(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, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_v)
    else
      call merid_face_thickness(v, h_in, h_l, h_r, bt_cont%h_v, dt, g, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_v)
    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, J, &
                            ish, ieh, do_I, vol_CFL)
  type(ocean_grid_type),        intent(inout) :: G        !< Ocean's grid structure.
  real, dimension(SZI_(G)),     intent(in)    :: v        !< Meridional velocity, in 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, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in) :: h_L    !< Left thickness in the
                                                          !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in) :: h_R    !< Right thickness in the
                                                          !! reconstruction, in H.
  real, dimension(SZI_(G)),     intent(inout) :: vh       !< Meridional mass or volume
                                                          !! transport, in H m2 s-1.
  real, dimension(SZI_(G)),     intent(inout) :: dvhdv    !< Partial derivative of vh
                                                          !! with v, in H m.
  real,                         intent(in)    :: dt       !< Time increment in s.
  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.
  ! Local variables
  real :: CFL ! The CFL number based on the local velocity and grid spacing, ND.
  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, in m.
  integer :: i

  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

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, LB, vol_CFL, &
                                marginal, visc_rem_v)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< Meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thickness used to
                                                                   !! calculate fluxes, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction, in H.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: h_v  !< Thickness at meridional faces,
                                                                   !! in H.
  real,                                      intent(in)    :: dt   !< Time increment in s.
  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(SZI_(G),SZJB_(G),SZK_(G)), intent(in), optional :: visc_rem_v !<
                          !! Both the fraction of the momentum originally in a
                          !! layer that remains after a time-step of viscosity,
                          !! and the fraction of a time-step's worth of a
                          !! barotropic acceleration that a layer experiences
                          !! after viscosity is applied. Non-dimensional between
                          !! 0 (at the bottom) and 1 (far above the bottom).
  ! Local variables
  real :: CFL ! The CFL number based on the local velocity and grid spacing, ND.
  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, in H.
  real :: h_marg ! The marginal thickness of a flux, in H.
  integer :: i, j, k, ish, ieh, jsh, jeh, nz
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke

!$OMP parallel default(none) shared(ish,ieh,jsh,jeh,nz,v,vol_CFL,dt,G, &
!$OMP                               h_L,h_R,h,h_v,visc_rem_v,marginal) &
!$OMP                       private(CFL,curv_3,h_marg,h_avg)
!$OMP do
  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 do
    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
!$OMP end parallel

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, CS, visc_rem, &
                             j, ish, ieh, do_I_in, full_precision, vh_3d)
  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< Meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction, in H.
  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),  optional :: vhbt !<
                       !! The summed volume flux through meridional faces, H m2 s-1.
  real, dimension(SZI_(G)),                  intent(in)    :: dv_max_CFL !< Maximum acceptable value
                       !! of dv, in m s-1.
  real, dimension(SZI_(G)),                  intent(in)    :: dv_min_CFL !< Minimum acceptable value
                       !! of dv, in m s-1.
  real, dimension(SZI_(G)),                  intent(in)    :: vh_tot_0    !<
                       !! The summed transport with 0 adjustment, in H m2 s-1.
  real, dimension(SZI_(G)),                  intent(in)    :: dvhdv_tot_0 !<
                       !! The partial derivative of dv_err with dv at 0 adjustment, in H m.
  real, dimension(SZI_(G)),                  intent(out)   :: dv !<
                       !! The barotropic velocity adjustment, in m s-1.
  real,                                      intent(in)    :: dt   !< Time increment in s.
  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 logical flag indicating which I values to work on.
  logical,                                   intent(in),    optional :: full_precision !<
                       !! full_precision - A flag indicating how carefully to iterate.  The
                       !! default is .true. (more accurate).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout), optional :: vh_3d !<
                       !! Volume flux through meridional faces = v*h*dx, H m2 s-1.
  ! Local variables
  real, dimension(SZI_(G),SZK_(G)) :: &
    vh_aux, &  ! An auxiliary meridional volume flux, in H m s-1.
    dvhdv      ! Partial derivative of vh with v, in H m.
  real, dimension(SZI_(G)) :: &
    vh_err, &  ! Difference between vhbt and the summed vh, in H m2 s-1.
    vh_err_best, & ! The smallest value of vh_err found so far, in H m2 s-1.
    v_new, &   ! The velocity with the correction added, in m s-1.
    dvhdv_tot,&! Summed partial derivative of vh with u, in H m.
    dv_min, &  ! Min/max limits on dv correction based on CFL limits
    dv_max     ! and previous iterations, in m s-1.
  real :: dv_prev ! The previous value of dv, in m s-1.
  real :: ddv    ! The change in dv from the previous iteration, in m s-1.
  real :: tol_eta ! The tolerance for the current iteration, in m.
  real :: tol_vel ! The tolerance for velocity in the current iteration, 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, j, ish, ieh, do_i, cs%vol_CFL)
    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, 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, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction, in H.
  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, in H m2 s-1.
  real, dimension(SZI_(G)),                  intent(in)    :: dvhdv_tot_0 !<
                       !! The partial derivative of du_err with dv at 0 adjustment, in H m.
  real, dimension(SZI_(G)),                  intent(in)    :: dv_max_CFL !< Maximum acceptable value
                       !! of dv, in m s-1.
  real, dimension(SZI_(G)),                  intent(in)    :: dv_min_CFL !< Minimum acceptable value
                       !! of dv, in m s-1.
  real,                                      intent(in)    :: dt   !< Time increment in s.
  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, m s-1.
    dvL, dvR, &   ! The barotropic velocity increments that give the southerly
                  ! (dvL) and northerly (dvR) test velocities.
    zeros, &      ! An array of full of 0's.
    dv_cfl, &     ! The velocity increment that corresponds to CFL_min, in 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, in m s-1.
    fa_marg_l, &  ! The effective layer marginal face areas with the southerly
    fa_marg_r, &  ! (_L), northerly (_R), and zero-barotropic (_0) test
    fa_marg_0, &  ! velocities, in H m.
    vh_l, vh_r, & ! The layer transports with the southerly (_L), northerly (_R)
    vh_0, &       ! and zero-barotropic (_0) test velocities, in H m2 s-1.
    famt_l, famt_r, & ! The summed effective marginal face areas for the 3
    famt_0, &     ! test velocities, in H m.
    vhtot_l, &    ! The summed transport with the southerly (vhtot_L) and
    vhtot_r       ! and northerly (vhtot_R) test velocities, in H m2 s-1.
  real :: FA_0    ! The effective face area with 0 barotropic transport, in m H.
  real :: FA_avg  ! The average effective face area, in m H, nominally given by
                  ! the realized transport divided by the barotropic velocity.
  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem, ND.  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, ND.
  real :: Idt     ! The inverse of the time step, in 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, &
                         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, &
                          fa_marg_0, visc_rem(:,k), dt, g, 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, &
                          fa_marg_l, visc_rem(:,k), dt, g, 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, &
                          fa_marg_r, visc_rem(:,k), dt, g, j, ish, ieh, do_i, &
                          cs%vol_CFL)
    do i=ish,ieh ; if (do_i(i)) then
      famt_0(i) = famt_0(i) + fa_marg_0(i)
      famt_l(i) = famt_l(i) + fa_marg_l(i)
      famt_r(i) = famt_r(i) + fa_marg_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%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
  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)
  type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_L  !< Left thickness in the
                                                         !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_R  !< Right thickness in the
                                                         !! reconstruction, in H.
  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.

  ! 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, stencil

  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
  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

    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 (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)
  type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_L  !< Left thickness in the
                                                         !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_R  !< Right thickness in the
                                                         !! reconstruction, in H.
  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.

  ! 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, stencil

  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
  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

    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 (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, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_L  !< Left thickness in the
                                                         !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_R  !< Right thickness in the
                                                         !! reconstruction, in H.
  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, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_L !< Left thickness in the
                                                                  !! reconstruction, in H.
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_R !< Right thickness in the
                                                                  !! reconstruction, in H.
  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, param_file, diag, CS)
  type(time_type), target, intent(in)    :: Time !< Time increment in s.
  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: GV   !< Vertical grid structure.
  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"
  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 \n"//&
                 "monotonic limiter.  The default (false) is to use a \n"//&
                 "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 \n"//&
                 "(arithmetic mean) interpolation of the edge values. \n"//&
                 "This may give better PV conservation propterties. While \n"//&
                 "it formally reduces the accuracy of the continuity \n"//&
                 "solver itself in the strongly advective limit, it does \n"//&
                 "not reduce the overall order of accuracy of the dynamic \n"//&
                 "core.", default=.false.)
  call get_param(param_file, mdl, "UPWIND_1ST_CONTINUITY", cs%upwind_1st, &
                 "If true, CONTINUITY_PPM becomes a 1st-order upwind \n"//&
                 "continuity solver.  This scheme is highly diffusive \n"//&
                 "but may be useful for debugging or in single-column \n"//&
                 "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 \n"//&
                 "barotropic and baroclinic estimates of the sea surface \n"//&
                 "height due to the fluxes through each face.  The total \n"//&
                 "tolerance for SSH is 4 times this value.  The default \n"//&
                 "is 0.5*NK*ANGSTROM, and this should not be set less x\n"//&
                 "than about 10^-15*MAXIMUM_DEPTH.", units="m", &
                 default=0.5*g%ke*gv%Angstrom_z)

  call get_param(param_file, mdl, "ETA_TOLERANCE_AUX", cs%tol_eta_aux, &
                 "The tolerance for free-surface height discrepancies \n"//&
                 "between the barotropic solution and the sum of the \n"//&
                 "layer thicknesses when calculating the auxiliary \n"//&
                 "corrected velocities. By default, this is the same as \n"//&
                 "ETA_TOLERANCE, but can be made larger for efficiency.", &
                 units="m", default=cs%tol_eta)
  call get_param(param_file, mdl, "VELOCITY_TOLERANCE", cs%tol_vel, &
                 "The tolerance for barotropic velocity discrepancies \n"//&
                 "between the barotropic solution and  the sum of the \n"//&
                 "layer thicknesses.", units="m s-1", default=3.0e8) ! 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 \n"//&
                 "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 \n"//&
                 "tracer cell areas when estimating CFL numbers.  The \n"//&
                 "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 \n"//&
                 "based criterion and only stop if the iteration is \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "solver for use as the weights in the barotropic solver. \n"//&
                 "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)

  cs%tol_eta = cs%tol_eta * gv%m_to_H
  cs%tol_eta_aux = cs%tol_eta_aux * gv%m_to_H

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