mom_continuity_ppm Module Reference

Detailed Description

Solve the layer continuity equation using the PPM method for layer fluxes.

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

Data Types
type	continuity_ppm_cs
	Control structure for mom_continuity_ppm. More...
type	loop_bounds_type
	A container for loop bounds. More...
integer	id_clock_update
	CPU time clock IDs.
integer	id_clock_correct
	CPU time clock IDs.
subroutine, public	continuity_ppm (u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, visc_rem_u, visc_rem_v, u_cor, v_cor, uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont)
	Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme, based on Lin (1994). More...
subroutine	zonal_mass_flux (u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, visc_rem_u, u_cor, uhbt_aux, u_cor_aux, BT_cont)
	Calculates the mass or volume fluxes through the zonal faces, and other related quantities. More...
subroutine	zonal_flux_layer (u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, j, ish, ieh, do_I, vol_CFL, OBC)
	Evaluates the zonal mass or volume fluxes in a layer. More...
subroutine	zonal_face_thickness (u, h, h_L, h_R, h_u, dt, G, LB, vol_CFL, marginal, visc_rem_u, OBC)
	Sets the effective interface thickness at each zonal velocity point. More...
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, OBC)
	Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport. More...
subroutine	set_zonal_bt_cont (u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, du_max_CFL, du_min_CFL, dt, G, US, CS, visc_rem, visc_rem_max, j, ish, ieh, do_I)
	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. More...
subroutine	meridional_mass_flux (v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, visc_rem_v, v_cor, vhbt_aux, v_cor_aux, BT_cont)
	Calculates the mass or volume fluxes through the meridional faces, and other related quantities. More...
subroutine	merid_flux_layer (v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, J, ish, ieh, do_I, vol_CFL, OBC)
	Evaluates the meridional mass or volume fluxes in a layer. More...
subroutine	merid_face_thickness (v, h, h_L, h_R, h_v, dt, G, LB, vol_CFL, marginal, visc_rem_v, OBC)
	Sets the effective interface thickness at each meridional velocity point. More...
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, OBC)
	Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport. More...
subroutine	set_merid_bt_cont (v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, dv_max_CFL, dv_min_CFL, dt, G, US, CS, visc_rem, visc_rem_max, j, ish, ieh, do_I)
	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. More...
subroutine	ppm_reconstruction_x (h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd, OBC)
	Calculates left/right edge values for PPM reconstruction. More...
subroutine	ppm_reconstruction_y (h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd, OBC)
	Calculates left/right edge values for PPM reconstruction. More...
subroutine	ppm_limit_pos (h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
	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. More...
subroutine	ppm_limit_cw84 (h_in, h_L, h_R, G, iis, iie, jis, jie)
	This subroutine limits the left/right edge values of the PPM reconstruction according to the monotonic prescription of Colella and Woodward, 1984. More...
real function	ratio_max (a, b, maxrat)
	Return the maximum ratio of a/b or maxrat. More...
subroutine, public	continuity_ppm_init (Time, G, GV, param_file, diag, CS)
	Initializes continuity_ppm_cs. More...
integer function, public	continuity_ppm_stencil (CS)
	continuity_PPM_stencil returns the continuity solver stencil size More...
subroutine, public	continuity_ppm_end (CS)
	Destructor for continuity_ppm_cs. More...

Function/Subroutine Documentation

continuity_ppm()

subroutine, public mom_continuity_ppm::continuity_ppm	(	real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)	u,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)	v,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	hin,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)	h,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)	uh,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out)	vh,
		real, intent(in)	dt,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(continuity_ppm_cs), pointer	CS,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in), optional	uhbt,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in), optional	vhbt,
		type(ocean_obc_type), optional, pointer	OBC,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in), optional	visc_rem_u,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in), optional	visc_rem_v,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out), optional	u_cor,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out), optional	v_cor,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in), optional	uhbt_aux,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in), optional	vhbt_aux,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out), optional	u_cor_aux,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out), optional	v_cor_aux,
		type(bt_cont_type), optional, pointer	BT_cont
	)

Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme, based on Lin (1994).

Parameters

[in,out]	g	The ocean's grid structure.
	cs	Module's control structure.
[in]	u	Zonal velocity [m s-1].
[in]	v	Meridional velocity [m s-1].
[in]	hin	Initial layer thickness [H ~> m or kg m-2].
[in,out]	h	Final layer thickness [H ~> m or kg m-2].
[out]	uh	Zonal volume flux, u*h*dy [H m2 s-1 ~> m3 s-1 or kg s-1].
[out]	vh	Meridional volume flux, v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	dt	Time increment [s].
[in]	gv	Vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	uhbt	The summed volume flux through zonal faces
[in]	vhbt	The summed volume flux through meridional faces
	obc	Open boundaries control structure.
[in]	visc_rem_u	The fraction of zonal momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]	visc_rem_v	The fraction of meridional momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[out]	u_cor	The zonal velocities that give uhbt as the depth-integrated transport [m s-1].
[out]	v_cor	The meridional velocities that give vhbt as the depth-integrated transport [m s-1].
[in]	uhbt_aux	A second set of summed volume fluxes through zonal faces [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	vhbt_aux	A second set of summed volume fluxes through meridional faces [H m2 s-1 ~> m3 s-1 or kg s-1].
[out]	u_cor_aux	The zonal velocities that give uhbt_aux as the depth-integrated transports [m s-1].
[out]	v_cor_aux	The meridional velocities that give vhbt_aux as the depth-integrated transports [m s-1].
	bt_cont	A structure with elements that describe the effective open face areas as a function of barotropic flow.

Definition at line 79 of file MOM_continuity_PPM.F90.

  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 [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v   !< Meridional velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in)    :: hin !< Initial layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(inout) :: h   !< Final layer thickness [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: uh  !< Zonal volume flux, u*h*dy [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(out)   :: vh  !< Meridional volume flux, v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
  real,                    intent(in)    :: dt  !< Time increment [s].
  type(verticalGrid_type), intent(in)    :: GV  !< Vertical grid structure.
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G)), &
                 optional, intent(in)    :: uhbt !< The summed volume flux through zonal faces
                                                 !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G)), &
                 optional, intent(in)    :: vhbt !< The summed volume flux through meridional faces
                                                 !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  type(ocean_OBC_type),  &
                 optional, pointer       :: OBC !< Open boundaries control structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(in)    :: visc_rem_u
                             !< The fraction of zonal momentum originally
                             !! in a layer that remains after a time-step of viscosity, and the
                             !! fraction of a time-step's worth of a barotropic acceleration that
                             !! a layer experiences after viscosity is applied.
                             !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(in)    :: visc_rem_v
                             !< The fraction of meridional momentum originally
                             !! in a layer that remains after a time-step of viscosity, and the
                             !! fraction of a time-step's worth of a barotropic acceleration that
                             !! a layer experiences after viscosity is applied.
                             !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: u_cor
                             !< The zonal velocities that give uhbt as the depth-integrated transport [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(out)   :: v_cor
                             !< The meridional velocities that give vhbt as the depth-integrated transport [m s-1].
  real, dimension(SZIB_(G),SZJ_(G)),  &
                 optional, intent(in)    :: uhbt_aux
                             !< A second set of summed volume fluxes through zonal faces
                             !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G)),  &
                 optional, intent(in)    :: vhbt_aux
                             !< A second set of summed volume fluxes through meridional faces
                             !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: u_cor_aux
                             !< The zonal velocities that give uhbt_aux as the depth-integrated
                             !! transports [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(out)   :: v_cor_aux
                             !< The meridional velocities that give vhbt_aux as the depth-integrated
                             !! transports [m s-1].
  type(BT_cont_type), optional, pointer  :: BT_cont !< A structure with elements that describe
                             !!  the effective open face areas as a function of barotropic flow.
 
  ! Local variables
  real :: h_min  ! The minimum layer thickness [H ~> m or kg m-2].  h_min could be 0.
  type(loop_bounds_type) :: LB
  integer :: is, ie, js, je, nz, stencil
  integer :: i, j, k
 
  logical :: x_first
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
 
  h_min = gv%Angstrom_H
 
  if (.not.associated(cs)) call mom_error(fatal, &
         "MOM_continuity_PPM: Module must be initialized before it is used.")
  x_first = (mod(g%first_direction,2) == 0)
 
  if (present(visc_rem_u) .neqv. present(visc_rem_v)) call mom_error(fatal, &
      "MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither"// &
      " one must be present in call to continuity_PPM.")
 
  stencil = 3 ; if (cs%simple_2nd) stencil = 2 ; if (cs%upwind_1st) stencil = 1
 
  if (x_first) then
  !    First, advect zonally.
    lb%ish = g%isc ; lb%ieh = g%iec
    lb%jsh = g%jsc-stencil ; lb%jeh = g%jec+stencil
    call zonal_mass_flux(u, hin, uh, dt, g, gv, us, cs, lb, uhbt, obc, visc_rem_u, &
                         u_cor, uhbt_aux, u_cor_aux, bt_cont)
 
    call cpu_clock_begin(id_clock_update)
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
      h(i,j,k) = hin(i,j,k) - dt* g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k))
  !   Uncomment this line to prevent underflow.
  !   if (h(i,j,k) < h_min) h(i,j,k) = h_min
    enddo ; enddo ; enddo
    call cpu_clock_end(id_clock_update)
 
    lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec
 
    !    Now advect meridionally, using the updated thicknesses to determine
    !  the fluxes.
    call meridional_mass_flux(v, h, vh, dt, g, gv, us, cs, lb, vhbt, obc, visc_rem_v, &
                              v_cor, vhbt_aux, v_cor_aux, bt_cont)
 
    call cpu_clock_begin(id_clock_update)
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
      h(i,j,k) = h(i,j,k) - dt*g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k))
  !   This line prevents underflow.
      if (h(i,j,k) < h_min) h(i,j,k) = h_min
    enddo ; enddo ; enddo
    call cpu_clock_end(id_clock_update)
 
  else  ! .not. x_first
  !    First, advect meridionally, so set the loop bounds accordingly.
    lb%ish = g%isc-stencil ; lb%ieh = g%iec+stencil
    lb%jsh = g%jsc ; lb%jeh = g%jec
 
    call meridional_mass_flux(v, hin, vh, dt, g, gv, us, cs, lb, vhbt, obc, visc_rem_v, &
                              v_cor, vhbt_aux, v_cor_aux, bt_cont)
 
    call cpu_clock_begin(id_clock_update)
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
      h(i,j,k) = hin(i,j,k) - dt*g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k))
    enddo ; enddo ; enddo
    call cpu_clock_end(id_clock_update)
 
  !    Now advect zonally, using the updated thicknesses to determine
  !  the fluxes.
    lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec
    call zonal_mass_flux(u, h, uh, dt, g, gv, us, cs, lb, uhbt, obc, visc_rem_u, &
                         u_cor, uhbt_aux, u_cor_aux, bt_cont)
 
    call cpu_clock_begin(id_clock_update)
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=lb%jsh,lb%jeh ; do i=lb%ish,lb%ieh
      h(i,j,k) = h(i,j,k) - dt* g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k))
      ! This line prevents underflow.
      if (h(i,j,k) < h_min) h(i,j,k) = h_min
    enddo ; enddo ; enddo
    call cpu_clock_end(id_clock_update)
 
  endif
 

continuity_ppm_end()

subroutine, public mom_continuity_ppm::continuity_ppm_end

(

type(continuity_ppm_cs), pointer

CS

)

Destructor for continuity_ppm_cs.

Parameters

cs	Module's control structure.

Definition at line 2347 of file MOM_continuity_PPM.F90.

  type(continuity_PPM_CS), pointer :: CS   !< Module's control structure.
  deallocate(cs)

continuity_ppm_init()

subroutine, public mom_continuity_ppm::continuity_ppm_init	(	type(time_type), intent(in), target	Time,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(continuity_ppm_cs), pointer	CS
	)

Initializes continuity_ppm_cs.

Parameters

[in]	time	Time increment [s].
[in]	g	The ocean's grid structure.
[in]	gv	Vertical grid structure.
[in]	param_file	A structure indicating the open file to parse for model parameter values.
[in,out]	diag	A structure that is used to regulate diagnostic output.
	cs	Module's control structure.

This include declares and sets the variable "version".

Definition at line 2245 of file MOM_continuity_PPM.F90.

  type(time_type), target, intent(in)    :: Time !< Time increment [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"
  real :: tol_eta_m  ! An unscaled version of tol_eta [m].
  character(len=40)  :: mdl = "MOM_continuity_PPM" ! This module's name.
 
  if (associated(cs)) then
    call mom_error(warning, "continuity_PPM_init called with associated control structure.")
    return
  endif
  allocate(cs)
 
! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", cs%monotonic, &
                 "If true, CONTINUITY_PPM uses the Colella and Woodward "//&
                 "monotonic limiter.  The default (false) is to use a "//&
                 "simple positive definite limiter.", default=.false.)
  call get_param(param_file, mdl, "SIMPLE_2ND_PPM_CONTINUITY", cs%simple_2nd, &
                 "If true, CONTINUITY_PPM uses a simple 2nd order "//&
                 "(arithmetic mean) interpolation of the edge values. "//&
                 "This may give better PV conservation properties. While "//&
                 "it formally reduces the accuracy of the continuity "//&
                 "solver itself in the strongly advective limit, it does "//&
                 "not reduce the overall order of accuracy of the dynamic "//&
                 "core.", default=.false.)
  call get_param(param_file, mdl, "UPWIND_1ST_CONTINUITY", cs%upwind_1st, &
                 "If true, CONTINUITY_PPM becomes a 1st-order upwind "//&
                 "continuity solver.  This scheme is highly diffusive "//&
                 "but may be useful for debugging or in single-column "//&
                 "mode where its minimal stencil is useful.", default=.false.)
  call get_param(param_file, mdl, "ETA_TOLERANCE", cs%tol_eta, &
                 "The tolerance for the differences between the "//&
                 "barotropic and baroclinic estimates of the sea surface "//&
                 "height due to the fluxes through each face.  The total "//&
                 "tolerance for SSH is 4 times this value.  The default "//&
                 "is 0.5*NK*ANGSTROM, and this should not be set less "//&
                 "than about 10^-15*MAXIMUM_DEPTH.", units="m", scale=gv%m_to_H, &
                 default=0.5*g%ke*gv%Angstrom_m, unscaled=tol_eta_m)
 
  call get_param(param_file, mdl, "ETA_TOLERANCE_AUX", cs%tol_eta_aux, &
                 "The tolerance for free-surface height discrepancies "//&
                 "between the barotropic solution and the sum of the "//&
                 "layer thicknesses when calculating the auxiliary "//&
                 "corrected velocities. By default, this is the same as "//&
                 "ETA_TOLERANCE, but can be made larger for efficiency.", &
                 units="m", default=tol_eta_m, scale=gv%m_to_H)
  call get_param(param_file, mdl, "VELOCITY_TOLERANCE", cs%tol_vel, &
                 "The tolerance for barotropic velocity discrepancies "//&
                 "between the barotropic solution and  the sum of the "//&
                 "layer thicknesses.", units="m s-1", default=3.0e8) ! The speed of light is the default.
 
  call get_param(param_file, mdl, "CONT_PPM_AGGRESS_ADJUST", cs%aggress_adjust,&
                 "If true, allow the adjusted velocities to have a "//&
                 "relative CFL change up to 0.5.", default=.false.)
  cs%vol_CFL = cs%aggress_adjust
  call get_param(param_file, mdl, "CONT_PPM_VOLUME_BASED_CFL", cs%vol_CFL, &
                 "If true, use the ratio of the open face lengths to the "//&
                 "tracer cell areas when estimating CFL numbers.  The "//&
                 "default is set by CONT_PPM_AGGRESS_ADJUST.", &
                 default=cs%aggress_adjust, do_not_read=cs%aggress_adjust)
  call get_param(param_file, mdl, "CONTINUITY_CFL_LIMIT", cs%CFL_limit_adjust, &
                 "The maximum CFL of the adjusted velocities.", units="nondim", &
                 default=0.5)
  call get_param(param_file, mdl, "CONT_PPM_BETTER_ITER", cs%better_iter, &
                 "If true, stop corrective iterations using a velocity "//&
                 "based criterion and only stop if the iteration is "//&
                 "better than all predecessors.", default=.true.)
  call get_param(param_file, mdl, "CONT_PPM_USE_VISC_REM_MAX", &
                                 cs%use_visc_rem_max, &
                 "If true, use more appropriate limiting bounds for "//&
                 "corrections in strongly viscous columns.", default=.true.)
  call get_param(param_file, mdl, "CONT_PPM_MARGINAL_FACE_AREAS", cs%marginal_faces, &
                 "If true, use the marginal face areas from the continuity "//&
                 "solver for use as the weights in the barotropic solver. "//&
                 "Otherwise use the transport averaged areas.", default=.true.)
 
  cs%diag => diag
 
  id_clock_update = cpu_clock_id('(Ocean continuity update)', grain=clock_routine)
  id_clock_correct = cpu_clock_id('(Ocean continuity correction)', grain=clock_routine)
 

continuity_ppm_stencil()

integer function, public mom_continuity_ppm::continuity_ppm_stencil

(

type(continuity_ppm_cs), pointer

CS

)

continuity_PPM_stencil returns the continuity solver stencil size

Parameters

cs	Module's control structure.

Returns

The continuity solver stencil size with the current settings.

Definition at line 2338 of file MOM_continuity_PPM.F90.

  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
 

merid_face_thickness()

subroutine mom_continuity_ppm::merid_face_thickness ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_L, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_R, real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)  h_v, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, type(loop_bounds_type), intent(in)  LB, logical, intent(in)  vol_CFL, logical, intent(in)  marginal, real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in), optional  visc_rem_v, type(ocean_obc_type), optional, pointer  OBC  )

private

Sets the effective interface thickness at each meridional velocity point.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	v	Meridional velocity [m s-1].
[in]	h	Layer thickness used to calculate fluxes, [H ~> m or kg m-2].
[in]	h_l	Left thickness in the reconstruction, [H ~> m or kg m-2].
[in]	h_r	Right thickness in the reconstruction, [H ~> m or kg m-2].
[in,out]	h_v	Thickness at meridional faces, [H ~> m or kg m-2].
[in]	dt	Time increment [s].
[in]	lb	Loop bounds structure.
[in]	vol_cfl	If true, rescale the ratio of face areas to the cell areas when estimating the CFL number.
[in]	marginal	If true, report the marginal face thicknesses; otherwise report transport-averaged thicknesses.
[in]	visc_rem_v	Both the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
	obc	Open boundaries control structure.

Definition at line 1444 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< Meridional velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thickness used to calculate fluxes,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the reconstruction,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the reconstruction,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: h_v  !< Thickness at meridional faces,
                                                                   !! [H ~> m or kg m-2].
  real,                                      intent(in)    :: dt   !< Time increment [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)), optional, intent(in) :: visc_rem_v !< Both the fraction
                          !! of the momentum originally in a layer that remains after a time-step of
                          !! viscosity, and the fraction of a time-step's worth of a barotropic
                          !! acceleration that a layer experiences after viscosity is applied.
                          !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  type(ocean_OBC_type),            optional, pointer :: OBC !< Open boundaries control structure.
 
  ! Local variables
  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h_in.
  real :: h_avg  ! The average thickness of a flux [H ~> m or kg m-2].
  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
  logical :: local_open_BC
  integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
 
  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
  do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
    if (v(i,j,k) > 0.0) then
      if (vol_cfl) then ; cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
      else ; cfl = v(i,j,k) * dt * g%IdyT(i,j) ; endif
      curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
      h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
      h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + &
                                3.0*curv_3*(cfl - 1.0))
    elseif (v(i,j,k) < 0.0) then
      if (vol_cfl) then ; cfl = (-v(i,j,k)*dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
      else ; cfl = -v(i,j,k) * dt * g%IdyT(i,j+1) ; endif
      curv_3 = h_l(i,j+1,k) + h_r(i,j+1,k) - 2.0*h(i,j+1,k)
      h_avg = h_l(i,j+1,k) + cfl * (0.5*(h_r(i,j+1,k)-h_l(i,j+1,k)) + curv_3*(cfl - 1.5))
      h_marg = h_l(i,j+1,k) + cfl * ((h_r(i,j+1,k)-h_l(i,j+1,k)) + &
                                    3.0*curv_3*(cfl - 1.0))
    else
      h_avg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
      !   The choice to use the arithmetic mean here is somewhat arbitrariy, but
      ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
      h_marg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
 !    h_marg = (2.0 * h_L(i,j+1,k) * h_R(i,j,k)) / &
 !             (h_L(i,j+1,k) + h_R(i,j,k) + GV%H_subroundoff)
    endif
 
    if (marginal) then ; h_v(i,j,k) = h_marg
    else ; h_v(i,j,k) = h_avg ; endif
  enddo ; enddo ; enddo
 
  if (present(visc_rem_v)) then
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=jsh-1,jeh ; do i=ish,ieh
      h_v(i,j,k) = h_v(i,j,k) * visc_rem_v(i,j,k)
    enddo ; enddo ; enddo
  endif
 
  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
  if (local_open_bc) then
    do n = 1, obc%number_of_segments
      if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
        j = obc%segment(n)%HI%JsdB
        if (obc%segment(n)%direction == obc_direction_n) then
          if (present(visc_rem_v)) then ; do k=1,nz
            do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
              h_v(i,j,k) = h(i,j,k) * visc_rem_v(i,j,k)
            enddo
          enddo ; else ; do k=1,nz
            do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
              h_v(i,j,k) = h(i,j,k)
            enddo
          enddo ; endif
        else
          if (present(visc_rem_v)) then ; do k=1,nz
            do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
              h_v(i,j,k) = h(i,j+1,k) * visc_rem_v(i,j,k)
            enddo
          enddo ; else ; do k=1,nz
            do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
              h_v(i,j,k) = h(i,j+1,k)
            enddo
          enddo ; endif
        endif
      endif
    enddo
  endif
 

merid_flux_layer()

subroutine mom_continuity_ppm::merid_flux_layer ( real, dimension(szi_(g)), intent(in)  v, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_L, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_R, real, dimension( g %isd: g %ied), intent(inout)  vh, real, dimension( g %isd: g %ied), intent(inout)  dvhdv, real, dimension(szi_(g)), intent(in)  visc_rem, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, integer, intent(in)  J, integer, intent(in)  ish, integer, intent(in)  ieh, logical, dimension( g %isd: g %ied), intent(in)  do_I, logical, intent(in)  vol_CFL, type(ocean_obc_type), optional, pointer  OBC  )

private

Evaluates the meridional mass or volume fluxes in a layer.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	v	Meridional velocity [m s-1].
[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).
[in]	h	Layer thickness used to calculate fluxes, [H ~> m or kg m-2].
[in]	h_l	Left thickness in the reconstruction [H ~> m or kg m-2].
[in]	h_r	Right thickness in the reconstruction [H ~> m or kg m-2].
[in,out]	vh	Meridional mass or volume transport [H m2 s-1 ~> m3 s-1 or kg s-1].
[in,out]	dvhdv	Partial derivative of vh with v [H m ~> m2 or kg m-1].
[in]	dt	Time increment [s].
[in]	j	Spatial index.
[in]	ish	Start of index range.
[in]	ieh	End of index range.
[in]	do_i	Which i values to work on.
[in]	vol_cfl	If true, rescale the ratio of face areas to the cell areas when estimating the CFL number.
	obc	Open boundaries control structure.

Definition at line 1364 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),        intent(inout) :: G        !< Ocean's grid structure.
  real, dimension(SZI_(G)),     intent(in)    :: v        !< Meridional velocity [m s-1].
  real, dimension(SZI_(G)),     intent(in)    :: visc_rem !< Both the fraction of the
         !! momentum originally in a layer that remains after a time-step
         !! of viscosity, and the fraction of a time-step's worth of a barotropic
         !! acceleration that a layer experiences after viscosity is applied.
         !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJ_(G)),  intent(in) :: h      !< Layer thickness used to calculate fluxes,
                                                          !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(in) :: h_L    !< Left thickness in the reconstruction
                                                          !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(in) :: h_R    !< Right thickness in the reconstruction
                                                          !! [H ~> m or kg m-2].
  real, dimension(SZI_(G)),     intent(inout) :: vh       !< Meridional mass or volume transport
                                                          !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G)),     intent(inout) :: dvhdv    !< Partial derivative of vh with v
                                                          !! [H m ~> m2 or kg m-1].
  real,                         intent(in)    :: dt       !< Time increment [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.
  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
  ! Local variables
  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h, i.e. [H ~> m or kg m-2].
  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
  integer :: i
  logical :: local_open_BC
 
  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
 
  do i=ish,ieh ; if (do_i(i)) then
    if (v(i) > 0.0) then
      if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
      else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
      curv_3 = h_l(i,j) + h_r(i,j) - 2.0*h(i,j)
      vh(i) = g%dx_Cv(i,j) * v(i) * ( h_r(i,j) + cfl * &
          (0.5*(h_l(i,j) - h_r(i,j)) + curv_3*(cfl - 1.5)) )
      h_marg = h_r(i,j) + cfl * ((h_l(i,j) - h_r(i,j)) + &
                                  3.0*curv_3*(cfl - 1.0))
    elseif (v(i) < 0.0) then
      if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
      else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
      curv_3 = h_l(i,j+1) + h_r(i,j+1) - 2.0*h(i,j+1)
      vh(i) = g%dx_Cv(i,j) * v(i) * ( h_l(i,j+1) + cfl * &
          (0.5*(h_r(i,j+1)-h_l(i,j+1)) + curv_3*(cfl - 1.5)) )
      h_marg = h_l(i,j+1) + cfl * ((h_r(i,j+1)-h_l(i,j+1)) + &
                                    3.0*curv_3*(cfl - 1.0))
    else
      vh(i) = 0.0
      h_marg = 0.5 * (h_l(i,j+1) + h_r(i,j))
    endif
    dvhdv(i) = g%dx_Cv(i,j) * h_marg * visc_rem(i)
  endif ; enddo
 
  if (local_open_bc) then
    do i=ish,ieh ; if (do_i(i)) then
      if (obc%segment(obc%segnum_v(i,j))%open) then
        if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
          vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j)
          dvhdv(i) = g%dx_Cv(i,j) * h(i,j) * visc_rem(i)
        else
          vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j+1)
          dvhdv(i) = g%dx_Cv(i,j) * h(i,j+1) * visc_rem(i)
        endif
      endif
    endif ; enddo
  endif

meridional_flux_adjust()

subroutine mom_continuity_ppm::meridional_flux_adjust ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_in, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_L, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_R, real, dimension( g %isd: g %ied), intent(in), optional  vhbt, real, dimension( g %isd: g %ied), intent(in)  vh_tot_0, real, dimension( g %isd: g %ied), intent(in)  dvhdv_tot_0, real, dimension( g %isd: g %ied), intent(out)  dv, real, dimension( g %isd: g %ied), intent(in)  dv_max_CFL, real, dimension( g %isd: g %ied), intent(in)  dv_min_CFL, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, type(continuity_ppm_cs), pointer  CS, real, dimension(szi_(g),szk_(g)), intent(in)  visc_rem, integer, intent(in)  j, integer, intent(in)  ish, integer, intent(in)  ieh, logical, dimension(szi_(g)), intent(in)  do_I_in, logical, intent(in), optional  full_precision, real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout), optional  vh_3d, type(ocean_obc_type), optional, pointer  OBC  )

private

Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	v	Meridional velocity [m s-1].
[in]	h_in	Layer thickness used to calculate fluxes [H ~> m or kg m-2].
[in]	h_l	Left thickness in the reconstruction [H ~> m or kg m-2].
[in]	h_r	Right thickness in the reconstruction [H ~> m or kg m-2].
[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).
[in]	vhbt	The summed volume flux through meridional faces
[in]	dv_max_cfl	Maximum acceptable value of dv [m s-1].
[in]	dv_min_cfl	Minimum acceptable value of dv [m s-1].
[in]	vh_tot_0	The summed transport with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	dvhdv_tot_0	The partial derivative of dv_err with dv at 0 adjustment [H m ~> m2 or kg m-1].
[out]	dv	The barotropic velocity adjustment [m s-1].
[in]	dt	Time increment [s].
	cs	This module's control structure.
[in]	j	Spatial index.
[in]	ish	Start of index range.
[in]	ieh	End of index range.
[in]	do_i_in	A flag indicating which I values to work on.
[in]	full_precision	A flag indicating how carefully to iterate. The default is .true. (more accurate).
[in,out]	vh_3d	Volume flux through meridional
	obc	Open boundaries control structure.

Definition at line 1552 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),  intent(inout) :: G   !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                         intent(in)    :: v    !< Meridional velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                         intent(in)    :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),&
                         intent(in)    :: h_L  !< Left thickness in the reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                         intent(in)    :: h_R  !< Right thickness in the reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)), intent(in) :: visc_rem
                             !< Both the fraction of the momentum originally
                             !! in a layer that remains after a time-step of viscosity, and the
                             !! fraction of a time-step's worth of a barotropic acceleration that
                             !! a layer experiences after viscosity is applied. Non-dimensional
                             !! between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G)), &
                  optional, intent(in)    :: vhbt !< The summed volume flux through meridional faces
                                                  !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G)), intent(in)    :: dv_max_CFL !< Maximum acceptable value of dv [m s-1].
  real, dimension(SZI_(G)), intent(in)    :: dv_min_CFL !< Minimum acceptable value of dv [m s-1].
  real, dimension(SZI_(G)), intent(in)    :: vh_tot_0   !< The summed transport with 0 adjustment
                                                        !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G)), intent(in)    :: dvhdv_tot_0 !< The partial derivative of dv_err with
                                                         !! dv at 0 adjustment [H m ~> m2 or kg m-1].
  real, dimension(SZI_(G)), intent(out)   :: dv   !< The barotropic velocity adjustment [m s-1].
  real,                     intent(in)    :: dt   !< Time increment [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 flag indicating which I values to work on.
  logical,        optional, intent(in)    :: full_precision !< A flag indicating how carefully to
                             !! iterate.  The default is .true. (more accurate).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                  optional, intent(inout) :: vh_3d !< Volume flux through meridional
                             !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
  type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
  ! Local variables
  real, dimension(SZI_(G),SZK_(G)) :: &
    vh_aux, &  ! An auxiliary meridional volume flux [H m s-1 ~> m2 s-1 or kg m-1 s-1].
    dvhdv      ! Partial derivative of vh with v [H m ~> m2 or kg m-1].
  real, dimension(SZI_(G)) :: &
    vh_err, &  ! Difference between vhbt and the summed vh [H m2 s-1 ~> m3 s-1 or kg s-1].
    vh_err_best, & ! The smallest value of vh_err found so far [H m2 s-1 ~> m3 s-1 or kg s-1].
    v_new, &   ! The velocity with the correction added [m s-1].
    dvhdv_tot,&! Summed partial derivative of vh with u [H m ~> m2 or kg m-1].
    dv_min, &  ! Min/max limits on dv correction based on CFL limits
    dv_max     ! and previous iterations [m s-1].
  real :: dv_prev ! The previous value of dv [m s-1].
  real :: ddv    ! The change in dv from the previous iteration [m s-1].
  real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
  real :: tol_vel ! The tolerance for velocity in the current iteration [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, obc)
    enddo ; endif
 
    if (itt < max_itts) then
      do i=ish,ieh
        vh_err(i) = -vhbt(i) ; dvhdv_tot(i) = 0.0
      enddo
      do k=1,nz ; do i=ish,ieh
        vh_err(i) = vh_err(i) + vh_aux(i,k)
        dvhdv_tot(i) = dvhdv_tot(i) + dvhdv(i,k)
      enddo ; enddo
      do i=ish,ieh
        vh_err_best(i) = min(vh_err_best(i), abs(vh_err(i)))
      enddo
    endif
  enddo ! itt-loop
  ! If there are any faces which have not converged to within the tolerance,
  ! so-be-it, or else use a final upwind correction?
  ! This never seems to happen with 20 iterations as max_itt.
 
  if (present(vh_3d)) then ; do k=1,nz ; do i=ish,ieh
    vh_3d(i,j,k) = vh_aux(i,k)
  enddo ; enddo ; endif
 

meridional_mass_flux()

subroutine mom_continuity_ppm::meridional_mass_flux ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_in, real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)  vh, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(continuity_ppm_cs), pointer  CS, type(loop_bounds_type), intent(in)  LB, real, dimension(szi_(g),szjb_(g)), intent(in), optional  vhbt, type(ocean_obc_type), optional, pointer  OBC, real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in), optional  visc_rem_v, real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out), optional  v_cor, real, dimension(szi_(g),szjb_(g)), intent(in), optional  vhbt_aux, real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out), optional  v_cor_aux, type(bt_cont_type), optional, pointer  BT_cont  )

private

Calculates the mass or volume fluxes through the meridional faces, and other related quantities.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	gv	Ocean's vertical grid structure.
[in]	v	Meridional velocity [m s-1].
[in]	h_in	Layer thickness used to calculate fluxes [H ~> m or kg m-2].
[out]	vh	Volume flux through meridional faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	dt	Time increment [s].
[in]	us	A dimensional unit scaling type
	cs	This module's control structure.
[in]	lb	Loop bounds structure.
	obc	Open boundary condition type specifies whether, where, and what open boundary conditions are used.
[in]	visc_rem_v	Both the fraction of the momentum
[in]	vhbt	The summed volume flux through meridional faces [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	vhbt_aux	A second set of summed volume fluxes through meridional faces [H m2 s-1 ~> m3 s-1 or kg s-1].
[out]	v_cor	The meridional velocitiess (v with a barotropic correction) that give vhbt as the depth-integrated transport [m s-1].
[out]	v_cor_aux	The meridional velocities (v with a barotropic correction) that give vhbt_aux as the depth-integrated transports [m s-1].
	bt_cont	A structure with elements that describe the effective open face areas as a function of barotropic flow.

Definition at line 1056 of file MOM_continuity_PPM.F90.

  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 [m s-1].
  real,  dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: vh   !< Volume flux through meridional
                                                                   !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
  real,                                      intent(in)    :: dt   !< Time increment [s].
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_PPM_CS),                   pointer       :: CS   !< This module's control structure.
  type(loop_bounds_type),                    intent(in)    :: LB   !< Loop bounds structure.
  type(ocean_OBC_type),            optional, pointer       :: OBC  !< Open boundary condition type
                                   !! specifies whether, where, and what open boundary conditions are used.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                   optional, intent(in)    :: visc_rem_v !< Both the fraction of the momentum
                                   !! originally in a layer that remains after a time-step of viscosity,
                                   !! and the fraction of a time-step's worth of a barotropic acceleration
                                   !! that a layer experiences after viscosity is applied.  Nondimensional between
                                   !! 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJB_(G)), optional, intent(in)  :: vhbt  !< The summed volume flux through
                                   !< meridional faces [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G)), optional, intent(in)  :: vhbt_aux !< A second set of summed volume fluxes
                                   !! through meridional faces [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                     optional, intent(out) :: v_cor
                                   !< The meridional velocitiess (v with a barotropic correction)
                                   !! that give vhbt as the depth-integrated transport [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                   optional, intent(out)   :: v_cor_aux
                                   !< The meridional velocities (v with a barotropic correction)
                                   !! that give vhbt_aux as the depth-integrated transports [m s-1].
  type(BT_cont_type),              optional, pointer       :: BT_cont !< A structure with elements that describe
                                   !! the effective open face areas as a function of barotropic flow.
  ! Local variables
  real, dimension(SZI_(G),SZK_(G)) :: &
    dvhdv      ! Partial derivative of vh with v [H m ~> m2 or kg m-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    h_L, h_R   ! Left and right face thicknesses [H ~> m or kg m-2].
  real, dimension(SZI_(G)) :: &
    dv, &      ! Corrective barotropic change in the velocity [m s-1].
    dv_min_CFL, & ! Min/max limits on dv correction
    dv_max_CFL, & ! to avoid CFL violations
    dvhdv_tot_0, & ! Summed partial derivative of vh with v [H m ~> m2 or kg m-1].
    vh_tot_0, &   ! Summed transport with no barotropic correction [H m2 s-1 ~> m3 s-1 or kg s-1].
    visc_rem_max  ! The column maximum of visc_rem.
  logical, dimension(SZI_(G)) :: do_I
  real, dimension(SZI_(G)) :: FAvi  ! A list of sums of meridional face areas [H m ~> m2 or kg m-1].
  real :: FA_v    ! A sum of meridional face areas [H m ~> m2 or kg m-1].
  real, dimension(SZI_(G),SZK_(G)) :: &
    visc_rem      ! A 2-D copy of visc_rem_v or an array of 1's.
  real :: I_vrm   ! 1.0 / visc_rem_max, nondim.
  real :: CFL_dt  ! The maximum CFL ratio of the adjusted velocities divided by
                  ! the time step [s-1].
  real :: I_dt    ! 1.0 / dt [s-1].
  real :: dv_lim  ! The velocity change that give a relative CFL of 1 [m s-1].
  real :: dy_N, dy_S ! Effective y-grid spacings to the north and south [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 => null()
 
  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_v_BCs_exist_globally
    local_open_bc = obc%open_v_BCs_exist_globally
  endif ; endif ; endif
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
 
  cfl_dt = cs%CFL_limit_adjust / dt
  i_dt = 1.0 / dt
  if (cs%aggress_adjust) cfl_dt = i_dt
 
  call cpu_clock_begin(id_clock_update)
!$OMP parallel do default(none) shared(nz,ish,ieh,jsh,jeh,h_in,h_L,h_R,G,GV,LB,CS,visc_rem,OBC)
  do k=1,nz
    ! This sets h_L and h_R.
    if (cs%upwind_1st) then
      do j=jsh-1,jeh+1 ; do i=ish,ieh
        h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
      enddo ; enddo
    else
      call ppm_reconstruction_y(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
                                2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
    endif
    do i=ish,ieh ; visc_rem(i,k) = 1.0 ; enddo
  enddo
  call cpu_clock_end(id_clock_update)
 
  call cpu_clock_begin(id_clock_correct)
!$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,v,h_in,h_L,h_R,vh,use_visc_rem, &
!$OMP                                  visc_rem_v,dt,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,FAvi,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, obc)
      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, obc=obc)
 
        if (present(v_cor)) then ; do k=1,nz
          do i=ish,ieh ; v_cor(i,j,k) = v(i,j,k) + dv(i) * visc_rem(i,k) ; enddo
          if (local_specified_bc) then ; do i=ish,ieh
            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., obc=obc)
 
        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, us, 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)) favi(i) = 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) .and. &
                (obc%segment(obc%segnum_v(i,j))%specified)) &
              favi(i) = favi(i) + &
                   obc%segment(obc%segnum_v(i,j))%normal_trans(i,j,k) / &
                   obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
          endif ; enddo ; enddo
          do i=ish,ieh ; if (do_i(i)) then
            bt_cont%FA_v_S0(i,j) = us%m_to_L*favi(i) ; bt_cont%FA_v_N0(i,j) = us%m_to_L*favi(i)
            bt_cont%FA_v_SS(i,j) = us%m_to_L*favi(i) ; bt_cont%FA_v_NN(i,j) = us%m_to_L*favi(i)
            bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
          endif ; enddo
        endif
      endif ! set_BT_cont
 
    endif ! present(vhbt) or do_aux or set_BT_cont
  enddo ! j-loop
 
  if (local_open_bc .and. set_bt_cont) then
    do n = 1, obc%number_of_segments
      if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S) then
        j = obc%segment(n)%HI%JsdB
        if (obc%segment(n)%direction == obc_direction_n) then
          do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
            fa_v = 0.0
            do k=1,nz ; fa_v = fa_v + h_in(i,j,k)*us%m_to_L*g%dx_Cv(i,j) ; enddo
            bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
            bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
            bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
          enddo
        else
          do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
            fa_v = 0.0
            do k=1,nz ; fa_v = fa_v + h_in(i,j+1,k)*us%m_to_L*g%dx_Cv(i,j) ; enddo
            bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
            bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
            bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
          enddo
        endif
      endif
    enddo
  endif
  call cpu_clock_end(id_clock_correct)
 
  if (set_bt_cont) then ; if (allocated(bt_cont%h_v)) then
    if (present(v_cor)) then
      call merid_face_thickness(v_cor, h_in, h_l, h_r, bt_cont%h_v, dt, g, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
    else
      call merid_face_thickness(v, h_in, h_l, h_r, bt_cont%h_v, dt, g, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
    endif
  endif ; endif
 

ppm_limit_cw84()

subroutine mom_continuity_ppm::ppm_limit_cw84 ( real, dimension(szi_(g),szj_(g)), intent(in)  h_in, real, dimension(szi_(g),szj_(g)), intent(inout)  h_L, real, dimension(szi_(g),szj_(g)), intent(inout)  h_R, type(ocean_grid_type), intent(in)  G, integer, intent(in)  iis, integer, intent(in)  iie, integer, intent(in)  jis, integer, intent(in)  jie  )

private

This subroutine limits the left/right edge values of the PPM reconstruction according to the monotonic prescription of Colella and Woodward, 1984.

Parameters

[in]	g	Ocean's grid structure.
[in]	h_in	Layer thickness [H ~> m or kg m-2].
[in,out]	h_l	Left thickness in the reconstruction, [H ~> m or kg m-2].
[in,out]	h_r	Right thickness in the reconstruction, [H ~> m or kg m-2].
[in]	iis	Start of i index range.
[in]	iie	End of i index range.
[in]	jis	Start of j index range.
[in]	jie	End of j index range.

Definition at line 2194 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),             intent(in)  :: G     !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in  !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_L !< Left thickness in the reconstruction,
                                                          !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_R !< Right thickness in the reconstruction,
                                                          !! [H ~> m or kg m-2].
  integer,                           intent(in)  :: iis   !< Start of i index range.
  integer,                           intent(in)  :: iie   !< End of i index range.
  integer,                           intent(in)  :: jis   !< Start of j index range.
  integer,                           intent(in)  :: jie   !< End of j index range.
 
  ! Local variables
  real    :: h_i, RLdiff, RLdiff2, RLmean, FunFac
  character(len=256) :: mesg
  integer :: i,j
 
  do j=jis,jie ; do i=iis,iie
    ! This limiter monotonizes the parabola following
    ! Colella and Woodward, 1984, Eq. 1.10
    h_i = h_in(i,j)
    if ( ( h_r(i,j) - h_i ) * ( h_i - h_l(i,j) ) <= 0. ) then
      h_l(i,j) = h_i ; h_r(i,j) = h_i
    else
      rldiff = h_r(i,j) - h_l(i,j)            ! Difference of edge values
      rlmean = 0.5 * ( h_r(i,j) + h_l(i,j) )  ! Mean of edge values
      funfac = 6. * rldiff * ( h_i - rlmean ) ! Some funny factor
      rldiff2 = rldiff * rldiff               ! Square of difference
      if ( funfac >  rldiff2 ) h_l(i,j) = 3. * h_i - 2. * h_r(i,j)
      if ( funfac < -rldiff2 ) h_r(i,j) = 3. * h_i - 2. * h_l(i,j)
    endif
  enddo ; enddo
 
  return

ppm_limit_pos()

subroutine mom_continuity_ppm::ppm_limit_pos ( real, dimension(szi_(g),szj_(g)), intent(in)  h_in, real, dimension(szi_(g),szj_(g)), intent(inout)  h_L, real, dimension(szi_(g),szj_(g)), intent(inout)  h_R, real, intent(in)  h_min, type(ocean_grid_type), intent(in)  G, integer, intent(in)  iis, integer, intent(in)  iie, integer, intent(in)  jis, integer, intent(in)  jie  )

private

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.

Parameters

[in]	g	Ocean's grid structure.
[in]	h_in	Layer thickness [H ~> m or kg m-2].
[in,out]	h_l	Left thickness in the reconstruction [H ~> m or kg m-2].
[in,out]	h_r	Right thickness in the reconstruction [H ~> m or kg m-2].
[in]	h_min	The minimum thickness that can be obtained by a concave parabolic fit.
[in]	iis	Start of i index range.
[in]	iie	End of i index range.
[in]	jis	Start of j index range.
[in]	jie	End of j index range.

Definition at line 2153 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_L !< Left thickness in the reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: h_R !< Right thickness in the reconstruction [H ~> m or kg m-2].
  real,                              intent(in)  :: h_min !< The minimum thickness
                    !! that can be obtained by a concave parabolic fit.
  integer,                           intent(in)  :: iis      !< Start of i index range.
  integer,                           intent(in)  :: iie      !< End of i index range.
  integer,                           intent(in)  :: jis      !< Start of j index range.
  integer,                           intent(in)  :: jie      !< End of j index range.
 
! Local variables
  real    :: curv, dh, scale
  character(len=256) :: mesg
  integer :: i,j
 
  do j=jis,jie ; do i=iis,iie
    ! This limiter prevents undershooting minima within the domain with
    ! values less than h_min.
    curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
    if (curv > 0.0) then ! Only minima are limited.
      dh = h_r(i,j) - h_l(i,j)
      if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
        if (h_in(i,j) <= h_min) then
          h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
        elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
          ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
          ! be limited in this case.  0 < scale < 1.
          scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
          h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
          h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
        endif
      endif
    endif
  enddo ; enddo
 

ppm_reconstruction_x()

subroutine mom_continuity_ppm::ppm_reconstruction_x ( real, dimension(szi_(g),szj_(g)), intent(in)  h_in, real, dimension(szi_(g),szj_(g)), intent(out)  h_L, real, dimension(szi_(g),szj_(g)), intent(out)  h_R, type(ocean_grid_type), intent(in)  G, type(loop_bounds_type), intent(in)  LB, real, intent(in)  h_min, logical, intent(in), optional  monotonic, logical, intent(in), optional  simple_2nd, type(ocean_obc_type), optional, pointer  OBC  )

private

Calculates left/right edge values for PPM reconstruction.

Parameters

[in]	g	Ocean's grid structure.
[in]	h_in	Layer thickness [H ~> m or kg m-2].
[out]	h_l	Left thickness in the reconstruction, [H ~> m or kg m-2].
[out]	h_r	Right thickness in the reconstruction, [H ~> m or kg m-2].
[in]	lb	Active loop bounds structure.
[in]	h_min	The minimum thickness that can be obtained by a concave parabolic fit.
[in]	monotonic	If true, use the Colella & Woodward monotonic limiter. Otherwise use a simple positive-definite limiter.
[in]	simple_2nd	If true, use the arithmetic mean thicknesses as the default edge values for a simple 2nd order scheme.
	obc	Open boundaries control structure.

Definition at line 1874 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_L  !< Left thickness in the reconstruction,
                                                         !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_R  !< Right thickness in the reconstruction,
                                                         !! [H ~> m or kg m-2].
  type(loop_bounds_type),            intent(in)  :: LB   !< Active loop bounds structure.
  real,                              intent(in)  :: h_min !< The minimum thickness
                    !! that can be obtained by a concave parabolic fit.
  logical, optional,                 intent(in)  :: monotonic !< If true, use the
                    !! Colella & Woodward monotonic limiter.
                    !! Otherwise use a simple positive-definite limiter.
  logical, optional,                 intent(in)  :: simple_2nd !< If true, use the
                    !! arithmetic mean thicknesses as the default edge values
                    !! for a simple 2nd order scheme.
  type(ocean_OBC_type),    optional, pointer     :: OBC !< Open boundaries control structure.
 
  ! Local variables with useful mnemonic names.
  real, dimension(SZI_(G),SZJ_(G))  :: slp ! The slopes.
  real, parameter :: oneSixth = 1./6.
  real :: h_ip1, h_im1
  real :: dMx, dMn
  logical :: use_CW84, use_2nd
  character(len=256) :: mesg
  integer :: i, j, isl, iel, jsl, jel, n, stencil
  logical :: local_open_BC
  type(OBC_segment_type), pointer :: segment => null()
 
  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
 
  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
 
  isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
 
  ! This is the stencil of the reconstruction, not the scheme overall.
  stencil = 2 ; if (use_2nd) stencil = 1
 
  if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
    write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
               & "x-halo that needs to be increased by ",i2,".")') &
               stencil + max(g%isd-isl,iel-g%ied)
    call mom_error(fatal,mesg)
  endif
  if ((jsl < g%jsd) .or. (jel > g%jed)) then
    write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
               & "y-halo that needs to be increased by ",i2,".")') &
               max(g%jsd-jsl,jel-g%jed)
    call mom_error(fatal,mesg)
  endif
 
  if (use_2nd) then
    do j=jsl,jel ; do i=isl,iel
      h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
      h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
      h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
      h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
    enddo ; enddo
  else
    do j=jsl,jel ; do i=isl-1,iel+1
      if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
        slp(i,j) = 0.0
      else
        ! This uses a simple 2nd order slope.
        slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
        ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
        dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
        dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
        slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
                ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
      endif
    enddo ; enddo
 
    if (local_open_bc) then
      do n=1, obc%number_of_segments
        segment => obc%segment(n)
        if (.not. segment%on_pe) cycle
        if (segment%direction == obc_direction_e .or. &
            segment%direction == obc_direction_w) then
          i=segment%HI%IsdB
          do j=segment%HI%jsd,segment%HI%jed
            slp(i+1,j) = 0.0
            slp(i,j) = 0.0
          enddo
        endif
      enddo
    endif
 
    do j=jsl,jel ; do i=isl,iel
      ! Neighboring values should take into account any boundaries.  The 3
      ! following sets of expressions are equivalent.
    ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
    ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
      h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
      h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
      ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
      h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
      h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
    enddo ; enddo
  endif
 
  if (local_open_bc) then
    do n=1, obc%number_of_segments
      segment => obc%segment(n)
      if (.not. segment%on_pe) cycle
      if (segment%direction == obc_direction_e) then
        i=segment%HI%IsdB
        do j=segment%HI%jsd,segment%HI%jed
          h_l(i+1,j) = h_in(i,j)
          h_r(i+1,j) = h_in(i,j)
          h_l(i,j) = h_in(i,j)
          h_r(i,j) = h_in(i,j)
        enddo
      elseif (segment%direction == obc_direction_w) then
        i=segment%HI%IsdB
        do j=segment%HI%jsd,segment%HI%jed
          h_l(i,j) = h_in(i+1,j)
          h_r(i,j) = h_in(i+1,j)
          h_l(i+1,j) = h_in(i+1,j)
          h_r(i+1,j) = h_in(i+1,j)
        enddo
      endif
    enddo
  endif
 
  if (use_cw84) then
    call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
  else
    call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
  endif
 
  return

ppm_reconstruction_y()

subroutine mom_continuity_ppm::ppm_reconstruction_y ( real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_in, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  h_L, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  h_R, type(ocean_grid_type), intent(in)  G, type(loop_bounds_type), intent(in)  LB, real, intent(in)  h_min, logical, intent(in), optional  monotonic, logical, intent(in), optional  simple_2nd, type(ocean_obc_type), optional, pointer  OBC  )

private

Calculates left/right edge values for PPM reconstruction.

Parameters

[in]	g	Ocean's grid structure.
[in]	h_in	Layer thickness [H ~> m or kg m-2].
[out]	h_l	Left thickness in the reconstruction, [H ~> m or kg m-2].
[out]	h_r	Right thickness in the reconstruction, [H ~> m or kg m-2].
[in]	lb	Active loop bounds structure.
[in]	h_min	The minimum thickness that can be obtained by a concave parabolic fit.
[in]	monotonic	If true, use the Colella & Woodward monotonic limiter. Otherwise use a simple positive-definite limiter.
[in]	simple_2nd	If true, use the arithmetic mean thicknesses as the default edge values for a simple 2nd order scheme.
	obc	Open boundaries control structure.

Definition at line 2013 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_L  !< Left thickness in the reconstruction,
                                                         !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_R  !< Right thickness in the reconstruction,
                                                         !! [H ~> m or kg m-2].
  type(loop_bounds_type),            intent(in)  :: LB   !< Active loop bounds structure.
  real,                              intent(in)  :: h_min !< The minimum thickness
                    !! that can be obtained by a concave parabolic fit.
  logical, optional,                 intent(in)  :: monotonic !< If true, use the
                    !! Colella & Woodward monotonic limiter.
                    !! Otherwise use a simple positive-definite limiter.
  logical, optional,                 intent(in)  :: simple_2nd !< If true, use the
                    !! arithmetic mean thicknesses as the default edge values
                    !! for a simple 2nd order scheme.
  type(ocean_OBC_type),    optional, pointer     :: OBC !< Open boundaries control structure.
 
  ! Local variables with useful mnemonic names.
  real, dimension(SZI_(G),SZJ_(G))  :: slp ! The slopes.
  real, parameter :: oneSixth = 1./6.
  real :: h_jp1, h_jm1
  real :: dMx, dMn
  logical :: use_CW84, use_2nd
  character(len=256) :: mesg
  integer :: i, j, isl, iel, jsl, jel, n, stencil
  logical :: local_open_BC
  type(OBC_segment_type), pointer :: segment => null()
 
  use_cw84 = .false. ; if (present(monotonic)) use_cw84 = monotonic
  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
 
  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
 
  isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
 
  ! This is the stencil of the reconstruction, not the scheme overall.
  stencil = 2 ; if (use_2nd) stencil = 1
 
  if ((isl < g%isd) .or. (iel > g%ied)) then
    write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
               & "x-halo that needs to be increased by ",i2,".")') &
               max(g%isd-isl,iel-g%ied)
    call mom_error(fatal,mesg)
  endif
  if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
    write(mesg,'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
                 & "y-halo that needs to be increased by ",i2,".")') &
                 stencil + max(g%jsd-jsl,jel-g%jed)
    call mom_error(fatal,mesg)
  endif
 
  if (use_2nd) then
    do j=jsl,jel ; do i=isl,iel
      h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
      h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
      h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
      h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
    enddo ; enddo
  else
    do j=jsl-1,jel+1 ; do i=isl,iel
      if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
        slp(i,j) = 0.0
      else
        ! This uses a simple 2nd order slope.
        slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
        ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
        dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
        dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
        slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
                ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
      endif
    enddo ; enddo
 
    if (local_open_bc) then
      do n=1, obc%number_of_segments
        segment => obc%segment(n)
        if (.not. segment%on_pe) cycle
        if (segment%direction == obc_direction_s .or. &
            segment%direction == obc_direction_n) then
          j=segment%HI%JsdB
          do i=segment%HI%isd,segment%HI%ied
            slp(i,j+1) = 0.0
            slp(i,j) = 0.0
          enddo
        endif
      enddo
    endif
 
    do j=jsl,jel ; do i=isl,iel
      ! Neighboring values should take into account any boundaries.  The 3
      ! following sets of expressions are equivalent.
      h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
      h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
      ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
      h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
      h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
    enddo ; enddo
  endif
 
  if (local_open_bc) then
    do n=1, obc%number_of_segments
      segment => obc%segment(n)
      if (.not. segment%on_pe) cycle
      if (segment%direction == obc_direction_n) then
        j=segment%HI%JsdB
        do i=segment%HI%isd,segment%HI%ied
          h_l(i,j+1) = h_in(i,j)
          h_r(i,j+1) = h_in(i,j)
          h_l(i,j) = h_in(i,j)
          h_r(i,j) = h_in(i,j)
        enddo
      elseif (segment%direction == obc_direction_s) then
        j=segment%HI%JsdB
        do i=segment%HI%isd,segment%HI%ied
          h_l(i,j) = h_in(i,j+1)
          h_r(i,j) = h_in(i,j+1)
          h_l(i,j+1) = h_in(i,j+1)
          h_r(i,j+1) = h_in(i,j+1)
        enddo
      endif
    enddo
  endif
 
  if (use_cw84) then
    call ppm_limit_cw84(h_in, h_l, h_r, g, isl, iel, jsl, jel)
  else
    call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
  endif
 
  return

ratio_max()

real function mom_continuity_ppm::ratio_max ( real, intent(in)  a, real, intent(in)  b, real, intent(in)  maxrat  )

private

Return the maximum ratio of a/b or maxrat.

Parameters

[in]	a	Numerator
[in]	b	Denominator
[in]	maxrat	Maximum value of ratio.

Returns

Return value.

Definition at line 2231 of file MOM_continuity_PPM.F90.

  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

set_merid_bt_cont()

subroutine mom_continuity_ppm::set_merid_bt_cont ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_in, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_L, real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_R, type(bt_cont_type), intent(inout)  BT_cont, real, dimension(szi_(g)), intent(in)  vh_tot_0, real, dimension(szi_(g)), intent(in)  dvhdv_tot_0, real, dimension(szi_(g)), intent(in)  dv_max_CFL, real, dimension(szi_(g)), intent(in)  dv_min_CFL, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, type(unit_scale_type), intent(in)  US, type(continuity_ppm_cs), pointer  CS, real, dimension( g %isd: g %ied, g %ke), intent(in)  visc_rem, real, dimension(szi_(g)), intent(in)  visc_rem_max, integer, intent(in)  j, integer, intent(in)  ish, integer, intent(in)  ieh, logical, dimension(szi_(g)), intent(in)  do_I  )

private

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.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	v	Meridional velocity [m s-1].
[in]	h_in	Layer thickness used to calculate fluxes, [H ~> m or kg m-2].
[in]	h_l	Left thickness in the reconstruction, [H ~> m or kg m-2].
[in]	h_r	Right thickness in the reconstruction, [H ~> m or kg m-2].
[in,out]	bt_cont	A structure with elements that describe the effective open face areas as a function of barotropic flow.
[in]	vh_tot_0	The summed transport with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	dvhdv_tot_0	The partial derivative of du_err with dv at 0 adjustment [H m ~> m2 or kg m-1].
[in]	dv_max_cfl	Maximum acceptable value of dv [m s-1].
[in]	dv_min_cfl	Minimum acceptable value of dv [m s-1].
[in]	dt	Time increment [s].
[in]	us	A dimensional unit scaling type
	cs	This module's control structure.
[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).
[in]	visc_rem_max	Maximum allowable visc_rem.
[in]	j	Spatial index.
[in]	ish	Start of index range.
[in]	ieh	End of index range.
[in]	do_i	A logical flag indicating which I values to work on.

Definition at line 1714 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< Meridional velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to calculate fluxes,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the reconstruction,
                                                                   !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the reconstruction,
                                                                   !! [H ~> m or kg m-2].
  type(BT_cont_type),                        intent(inout) :: BT_cont !< A structure with elements
                       !! that describe the effective open face areas as a function of barotropic flow.
  real, dimension(SZI_(G)),                  intent(in)    :: vh_tot_0    !< The summed transport
                       !! with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G)),                  intent(in)    :: dvhdv_tot_0 !< The partial derivative
                       !! of du_err with dv at 0 adjustment [H m ~> m2 or kg m-1].
  real, dimension(SZI_(G)),                  intent(in)    :: dv_max_CFL !< Maximum acceptable value of dv [m s-1].
  real, dimension(SZI_(G)),                  intent(in)    :: dv_min_CFL !< Minimum acceptable value of dv [m s-1].
  real,                                      intent(in)    :: dt   !< Time increment [s].
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_PPM_CS),                   pointer       :: CS   !< This module's control structure.
  real, dimension(SZI_(G),SZK_(G)),          intent(in)    :: visc_rem !< Both the fraction of the
                       !! momentum originally in a layer that remains after a time-step
                       !! of viscosity, and the fraction of a time-step's worth of a barotropic
                       !! acceleration that a layer experiences after viscosity is applied.
                       !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G)),                  intent(in)    :: visc_rem_max !< Maximum allowable visc_rem.
  integer,                                   intent(in)    :: j        !< Spatial index.
  integer,                                   intent(in)    :: ish      !< Start of index range.
  integer,                                   intent(in)    :: ieh      !< End of index range.
  logical, dimension(SZI_(G)),               intent(in)    :: do_I     !< A logical flag indicating
                       !! which I values to work on.
  ! Local variables
  real, dimension(SZI_(G)) :: &
    dv0, &        ! The barotropic velocity increment that gives 0 transport [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 [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 [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 [H m ~> m2 or kg m-1].
    vh_l, vh_r, & ! The layer transports with the southerly (_L), northerly (_R)
    vh_0, &       ! and zero-barotropic (_0) test velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
    famt_l, famt_r, & ! The summed effective marginal face areas for the 3
    famt_0, &     ! test velocities [H m ~> m2 or kg m-1].
    vhtot_l, &    ! The summed transport with the southerly (vhtot_L) and
    vhtot_r       ! and northerly (vhtot_R) test velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
  real :: FA_0    ! The effective face area with 0 barotropic transport [H m ~> m2 or kg m-1].
  real :: FA_avg  ! The average effective face area [H m ~> m2 or kg m-1], nominally given by
                  ! the realized transport divided by the barotropic velocity.
  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]  This
                       ! limiting is necessary to keep the inverse of visc_rem
                       ! from leading to large CFL numbers.
  real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
                       ! in finding the barotropic velocity that changes the
                       ! flow direction.  This is necessary to keep the inverse
                       ! of visc_rem from leading to large CFL numbers.
  real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
                  ! flow is truly upwind [nondim]
  real :: Idt     ! The inverse of the time step [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) = us%m_to_L*fa_0 ; bt_cont%FA_v_SS(i,j) = us%m_to_L*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) = us%m_s_to_L_T*(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) = us%m_to_L*fa_0 ; bt_cont%FA_v_NN(i,j) = us%m_to_L*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) = us%m_s_to_L_T*(1.5 * (dvr(i) - dv0(i))) * &
                   ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
    endif
  else
    bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
    bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
    bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
  endif ; enddo
 

set_zonal_bt_cont()

subroutine mom_continuity_ppm::set_zonal_bt_cont ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_in, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_L, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_R, type(bt_cont_type), intent(inout)  BT_cont, real, dimension( g %isdb: g %iedb), intent(in)  uh_tot_0, real, dimension( g %isdb: g %iedb), intent(in)  duhdu_tot_0, real, dimension( g %isdb: g %iedb), intent(in)  du_max_CFL, real, dimension( g %isdb: g %iedb), intent(in)  du_min_CFL, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, type(unit_scale_type), intent(in)  US, type(continuity_ppm_cs), pointer  CS, real, dimension(szib_(g),szk_(g)), intent(in)  visc_rem, real, dimension( g %isdb: g %iedb), intent(in)  visc_rem_max, integer, intent(in)  j, integer, intent(in)  ish, integer, intent(in)  ieh, logical, dimension( g %isdb: g %iedb), intent(in)  do_I  )

private

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.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	u	Zonal velocity [m s-1].
[in]	h_in	Layer thickness used to calculate fluxes [H ~> m or kg m-2].
[in]	h_l	Left thickness in the reconstruction [H ~> m or kg m-2].
[in]	h_r	Right thickness in the reconstruction [H ~> m or kg m-2].
[in,out]	bt_cont	A structure with elements that describe the effective open face areas as a function of barotropic flow.
[in]	uh_tot_0	The summed transport with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	duhdu_tot_0	The partial derivative of du_err with du at 0 adjustment [H m ~> m2 or kg m-1].
[in]	du_max_cfl	Maximum acceptable value of du [m s-1].
[in]	du_min_cfl	Minimum acceptable value of du [m s-1].
[in]	dt	Time increment [s].
[in]	us	A dimensional unit scaling type
	cs	This module's control structure.
[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).
[in]	visc_rem_max	Maximum allowable visc_rem.
[in]	j	Spatial index.
[in]	ish	Start of index range.
[in]	ieh	End of index range.
[in]	do_i	A logical flag indicating which I values to work on.

Definition at line 891 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< Zonal velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  type(BT_cont_type),                        intent(inout) :: BT_cont !< A structure with elements
                       !! that describe the effective open face areas as a function of barotropic flow.
  real, dimension(SZIB_(G)),                 intent(in)    :: uh_tot_0    !< The summed transport
                       !! with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: duhdu_tot_0 !< The partial derivative
                       !! of du_err with du at 0 adjustment [H m ~> m2 or kg m-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: du_max_CFL  !< Maximum acceptable
                       !! value of du [m s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: du_min_CFL  !< Minimum acceptable
                       !! value of du [m s-1].
  real,                                      intent(in)    :: dt   !< Time increment [s].
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_PPM_CS),                   pointer       :: CS   !< This module's control structure.
  real, dimension(SZIB_(G),SZK_(G)),         intent(in)    :: visc_rem !< Both the fraction of the
                       !! momentum originally in a layer that remains after a time-step of viscosity, and
                       !! the fraction of a time-step's worth of a barotropic acceleration that a layer
                       !! experiences after viscosity is applied.
                       !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G)),                 intent(in)    :: visc_rem_max !< Maximum allowable visc_rem.
  integer,                                   intent(in)    :: j        !< Spatial index.
  integer,                                   intent(in)    :: ish      !< Start of index range.
  integer,                                   intent(in)    :: ieh      !< End of index range.
  logical, dimension(SZIB_(G)),              intent(in)    :: do_I     !< A logical flag indicating
                       !! which I values to work on.
  ! Local variables
  real, dimension(SZIB_(G)) :: &
    du0, &        ! The barotropic velocity increment that gives 0 transport [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 [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 [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 [H m ~> m2 or kg m-1].
    uh_l, uh_r, & ! The layer transports with the westerly (_L), easterly (_R),
    uh_0, &       ! and zero-barotropic (_0) test velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
    famt_l, famt_r, & ! The summed effective marginal face areas for the 3
    famt_0, &     ! test velocities [H m ~> m2 or kg m-1].
    uhtot_l, &    ! The summed transport with the westerly (uhtot_L) and
    uhtot_r       ! and easterly (uhtot_R) test velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
  real :: FA_0    ! The effective face area with 0 barotropic transport [m H ~> m2 or kg m].
  real :: FA_avg  ! The average effective face area [m H ~> m2 or kg m], nominally given by
                  ! the realized transport divided by the barotropic velocity.
  real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]  This
                       ! limiting is necessary to keep the inverse of visc_rem
                       ! from leading to large CFL numbers.
  real :: min_visc_rem ! The smallest permitted value for visc_rem that is used
                       ! in finding the barotropic velocity that changes the
                       ! flow direction.  This is necessary to keep the inverse
                       ! of visc_rem from leading to large CFL numbers.
  real :: CFL_min ! A minimal increment in the CFL to try to ensure that the
                  ! flow is truly upwind [nondim]
  real :: Idt     ! The inverse of the time step [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%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
      bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
    enddo ; enddo
    return
  endif
 
  do k=1,nz ; do i=ish-1,ieh ; if (do_i(i)) then
    visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
    if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points.
      if (u(i,j,k) + dur(i)*visc_rem_lim > -du_cfl(i)*visc_rem(i,k)) &
        dur(i) = -(u(i,j,k) + du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
      if (u(i,j,k) + dul(i)*visc_rem_lim < du_cfl(i)*visc_rem(i,k)) &
        dul(i) = -(u(i,j,k) - du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
    endif
  endif ; enddo ; enddo
 
  do k=1,nz
    do i=ish-1,ieh ; if (do_i(i)) then
      u_l(i) = u(i,j,k) + dul(i) * visc_rem(i,k)
      u_r(i) = u(i,j,k) + dur(i) * visc_rem(i,k)
      u_0(i) = u(i,j,k) + du0(i) * visc_rem(i,k)
    endif ; enddo
    call zonal_flux_layer(u_0, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_0, &
                          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) = us%m_to_L*fa_0 ; bt_cont%FA_u_WW(i,j) = us%m_to_L*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) = us%m_s_to_L_T*(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) = us%m_to_L*fa_0 ; bt_cont%FA_u_EE(i,j) = us%m_to_L*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) = us%m_s_to_L_T*(1.5 * (dur(i) - du0(i))) * &
                            ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
    endif
  else
    bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
    bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
    bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
  endif ; enddo
 

zonal_face_thickness()

subroutine mom_continuity_ppm::zonal_face_thickness ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_L, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_R, real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)  h_u, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, type(loop_bounds_type), intent(in)  LB, logical, intent(in)  vol_CFL, logical, intent(in)  marginal, real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in), optional  visc_rem_u, type(ocean_obc_type), optional, pointer  OBC  )

private

Sets the effective interface thickness at each zonal velocity point.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	u	Zonal velocity [m s-1].
[in]	h	Layer thickness used to calculate fluxes [H ~> m or kg m-2].
[in]	h_l	Left thickness in the reconstruction [H ~> m or kg m-2].
[in]	h_r	Right thickness in the reconstruction [H ~> m or kg m-2].
[in,out]	h_u	Thickness at zonal faces [H ~> m or kg m-2].
[in]	dt	Time increment [s].
[in]	lb	Loop bounds structure.
[in]	vol_cfl	If true, rescale the ratio of face areas to the cell areas when estimating the CFL number.
[in]	marginal	If true, report the marginal face thicknesses; otherwise report transport-averaged thicknesses.
[in]	visc_rem_u	Both the fraction of the momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
	obc	Open boundaries control structure.

Definition at line 621 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< Zonal velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thickness used to
                                                                   !! calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_u  !< Thickness at zonal faces [H ~> m or kg m-2].
  real,                                      intent(in)    :: dt   !< Time increment [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)), &
                                   optional, intent(in)    :: visc_rem_u
                          !< Both the fraction of the momentum originally in a layer that remains after
                          !! a time-step of viscosity, and the fraction of a time-step's worth of a
                          !! barotropic acceleration that a layer experiences after viscosity is applied.
                          !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  type(ocean_OBC_type),            optional, pointer       :: OBC !< Open boundaries control structure.
 
  ! Local variables
  real :: CFL  ! The CFL number based on the local velocity and grid spacing [nondim]
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h_in.
  real :: h_avg  ! The average thickness of a flux [H ~> m or kg m-2].
  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
  logical :: local_open_BC
  integer :: i, j, k, ish, ieh, jsh, jeh, nz, n
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
 
  !$OMP parallel do default(shared) private(CFL,curv_3,h_marg,h_avg)
  do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
    if (u(i,j,k) > 0.0) then
      if (vol_cfl) then ; cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
      else ; cfl = u(i,j,k) * dt * g%IdxT(i,j) ; endif
      curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
      h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
      h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + 3.0*curv_3*(cfl - 1.0))
    elseif (u(i,j,k) < 0.0) then
      if (vol_cfl) then ; cfl = (-u(i,j,k)*dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
      else ; cfl = -u(i,j,k) * dt * g%IdxT(i+1,j) ; endif
      curv_3 = h_l(i+1,j,k) + h_r(i+1,j,k) - 2.0*h(i+1,j,k)
      h_avg = h_l(i+1,j,k) + cfl * (0.5*(h_r(i+1,j,k)-h_l(i+1,j,k)) + curv_3*(cfl - 1.5))
      h_marg = h_l(i+1,j,k) + cfl * ((h_r(i+1,j,k)-h_l(i+1,j,k)) + &
                                    3.0*curv_3*(cfl - 1.0))
    else
      h_avg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
      !   The choice to use the arithmetic mean here is somewhat arbitrariy, but
      ! it should be noted that h_L(i+1,j,k) and h_R(i,j,k) are usually the same.
      h_marg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
 !    h_marg = (2.0 * h_L(i+1,j,k) * h_R(i,j,k)) / &
 !             (h_L(i+1,j,k) + h_R(i,j,k) + GV%H_subroundoff)
    endif
 
    if (marginal) then ; h_u(i,j,k) = h_marg
    else ; h_u(i,j,k) = h_avg ; endif
  enddo ; enddo ; enddo
  if (present(visc_rem_u)) then
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=jsh,jeh ; do i=ish-1,ieh
      h_u(i,j,k) = h_u(i,j,k) * visc_rem_u(i,j,k)
    enddo ; enddo ; enddo
  endif
 
  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
  if (local_open_bc) then
    do n = 1, obc%number_of_segments
      if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then
        i = obc%segment(n)%HI%IsdB
        if (obc%segment(n)%direction == obc_direction_e) then
          if (present(visc_rem_u)) then ; do k=1,nz
            do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
              h_u(i,j,k) = h(i,j,k) * visc_rem_u(i,j,k)
            enddo
          enddo ; else ; do k=1,nz
            do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
              h_u(i,j,k) = h(i,j,k)
            enddo
          enddo ; endif
        else
          if (present(visc_rem_u)) then ; do k=1,nz
            do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
              h_u(i,j,k) = h(i+1,j,k) * visc_rem_u(i,j,k)
            enddo
          enddo ; else ; do k=1,nz
            do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
              h_u(i,j,k) = h(i+1,j,k)
            enddo
          enddo ; endif
        endif
      endif
    enddo
  endif
 

zonal_flux_adjust()

subroutine mom_continuity_ppm::zonal_flux_adjust ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_in, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_L, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_R, real, dimension(szib_(g)), intent(in), optional  uhbt, real, dimension(szib_(g)), intent(in)  uh_tot_0, real, dimension(szib_(g)), intent(in)  duhdu_tot_0, real, dimension(szib_(g)), intent(out)  du, real, dimension(szib_(g)), intent(in)  du_max_CFL, real, dimension(szib_(g)), intent(in)  du_min_CFL, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, type(continuity_ppm_cs), pointer  CS, real, dimension( g %isdb: g %iedb, g %ke), intent(in)  visc_rem, integer, intent(in)  j, integer, intent(in)  ish, integer, intent(in)  ieh, logical, dimension( g %isdb: g %iedb), intent(in)  do_I_in, logical, intent(in), optional  full_precision, real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout), optional  uh_3d, type(ocean_obc_type), optional, pointer  OBC  )

private

Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	u	Zonal velocity [m s-1].
[in]	h_in	Layer thickness used to calculate fluxes [H ~> m or kg m-2].
[in]	h_l	Left thickness in the reconstruction [H ~> m or kg m-2].
[in]	h_r	Right thickness in the reconstruction [H ~> m or kg m-2].
[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).
[in]	uhbt	The summed volume flux through zonal faces [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	du_max_cfl	Maximum acceptable value of du [m s-1].
[in]	du_min_cfl	Minimum acceptable value of du [m s-1].
[in]	uh_tot_0	The summed transport with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
[in]	duhdu_tot_0	The partial derivative of du_err with du at 0 adjustment [H m ~> m2 or kg m-1].
[out]	du	The barotropic velocity adjustment [m s-1].
[in]	dt	Time increment [s].
	cs	This module's control structure.
[in]	j	Spatial index.
[in]	ish	Start of index range.
[in]	ieh	End of index range.
[in]	do_i_in	A logical flag indicating which I values to work on.
[in]	full_precision	A flag indicating how carefully to iterate. The default is .true. (more accurate).
[in,out]	uh_3d	Volume flux through zonal faces = u*h*dy [H m2 s-1 ~> m3 s-1 or kg s-1].
	obc	Open boundaries control structure.

Definition at line 728 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),                     intent(inout) :: G    !< Ocean's grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< Zonal velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_in !< Layer thickness used to
                                                                   !! calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_L  !< Left thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h_R  !< Right thickness in the
                                                                   !! reconstruction [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZK_(G)),         intent(in)    :: visc_rem !< Both the fraction of the
                       !! momentum originally in a layer that remains after a time-step of viscosity, and
                       !! the fraction of a time-step's worth of a barotropic acceleration that a layer
                       !! experiences after viscosity is applied.
                       !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G)),       optional, intent(in)    :: uhbt !< The summed volume flux
                       !! through zonal faces [H m2 s-1 ~> m3 s-1 or kg s-1].
 
  real, dimension(SZIB_(G)),                 intent(in)    :: du_max_CFL  !< Maximum acceptable
                       !! value of du [m s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: du_min_CFL  !< Minimum acceptable
                       !! value of du [m s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: uh_tot_0    !< The summed transport
                       !! with 0 adjustment [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G)),                 intent(in)    :: duhdu_tot_0 !< The partial derivative
                       !! of du_err with du at 0 adjustment [H m ~> m2 or kg m-1].
  real, dimension(SZIB_(G)),                 intent(out)   :: du !<
                       !! The barotropic velocity adjustment [m s-1].
  real,                                      intent(in)    :: dt   !< Time increment [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,                         optional, intent(in)    :: full_precision !<
                       !! A flag indicating how carefully to iterate.  The
                       !! default is .true. (more accurate).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), optional, intent(inout) :: uh_3d !<
                       !! Volume flux through zonal faces = u*h*dy [H m2 s-1 ~> m3 s-1 or kg s-1].
  type(ocean_OBC_type),            optional, pointer       :: OBC !< Open boundaries control structure.
  ! Local variables
  real, dimension(SZIB_(G),SZK_(G)) :: &
    uh_aux, &  ! An auxiliary zonal volume flux [H m s-1 ~> m2 s-1 or kg m-1 s-1].
    duhdu      ! Partial derivative of uh with u [H m ~> m2 or kg m-1].
  real, dimension(SZIB_(G)) :: &
    uh_err, &  ! Difference between uhbt and the summed uh [H m2 s-1 ~> m3 s-1 or kg s-1].
    uh_err_best, & ! The smallest value of uh_err found so far [H m2 s-1 ~> m3 s-1 or kg s-1].
    u_new, &   ! The velocity with the correction added [m s-1].
    duhdu_tot,&! Summed partial derivative of uh with u [H m ~> m2 or kg m-1].
    du_min, &  ! Min/max limits on du correction based on CFL limits
    du_max     ! and previous iterations [m s-1].
  real :: du_prev ! The previous value of du [m s-1].
  real :: ddu    ! The change in du from the previous iteration [m s-1].
  real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
  real :: tol_vel ! The tolerance for velocity in the current iteration [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, obc)
    enddo ; endif
 
    if (itt < max_itts) then
      do i=ish-1,ieh
        uh_err(i) = -uhbt(i) ; duhdu_tot(i) = 0.0
      enddo
      do k=1,nz ; do i=ish-1,ieh
        uh_err(i) = uh_err(i) + uh_aux(i,k)
        duhdu_tot(i) = duhdu_tot(i) + duhdu(i,k)
      enddo ; enddo
      do i=ish-1,ieh
        uh_err_best(i) = min(uh_err_best(i), abs(uh_err(i)))
      enddo
    endif
  enddo ! itt-loop
  ! If there are any faces which have not converged to within the tolerance,
  ! so-be-it, or else use a final upwind correction?
  ! This never seems to happen with 20 iterations as max_itt.
 
  if (present(uh_3d)) then ; do k=1,nz ; do i=ish-1,ieh
    uh_3d(i,j,k) = uh_aux(i,k)
  enddo ; enddo ; endif
 

zonal_flux_layer()

subroutine mom_continuity_ppm::zonal_flux_layer ( real, dimension( g %isdb: g %iedb), intent(in)  u, real, dimension(szi_(g)), intent(in)  h, real, dimension(szi_(g)), intent(in)  h_L, real, dimension(szi_(g)), intent(in)  h_R, real, dimension(szib_(g)), intent(inout)  uh, real, dimension(szib_(g)), intent(inout)  duhdu, real, dimension( g %isdb: g %iedb), intent(in)  visc_rem, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, integer, intent(in)  j, integer, intent(in)  ish, integer, intent(in)  ieh, logical, dimension(szib_(g)), intent(in)  do_I, logical, intent(in)  vol_CFL, type(ocean_obc_type), optional, pointer  OBC  )

private

Evaluates the zonal mass or volume fluxes in a layer.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	u	Zonal velocity [m s-1].
[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).
[in]	h	Layer thickness [H ~> m or kg m-2].
[in]	h_l	Left thickness [H ~> m or kg m-2].
[in]	h_r	Right thickness [H ~> m or kg m-2].
[in,out]	uh	Zonal mass or volume transport [H m2 s-1 ~> m3 s-1 or kg s-1].
[in,out]	duhdu	Partial derivative of uh with u [H m ~> m2 or kg m-1].
[in]	dt	Time increment [s].
[in]	j	Spatial index.
[in]	ish	Start of index range.
[in]	ieh	End of index range.
[in]	do_i	Which i values to work on.
[in]	vol_cfl	If true, rescale the ratio of face areas to the cell areas when estimating the CFL number.
	obc	Open boundaries control structure.

Definition at line 545 of file MOM_continuity_PPM.F90.

  type(ocean_grid_type),        intent(inout) :: G        !< Ocean's grid structure.
  real, dimension(SZIB_(G)),    intent(in)    :: u        !< Zonal velocity [m s-1].
  real, dimension(SZIB_(G)),    intent(in)    :: visc_rem !< Both the fraction of the
                        !! momentum originally in a layer that remains after a time-step
                        !! of viscosity, and the fraction of a time-step's worth of a barotropic
                        !! acceleration that a layer experiences after viscosity is applied.
                        !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G)),     intent(in)    :: h        !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G)),     intent(in)    :: h_L      !< Left thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G)),     intent(in)    :: h_R      !< Right thickness [H ~> m or kg m-2].
  real, dimension(SZIB_(G)),    intent(inout) :: uh       !< Zonal mass or volume
                                                          !! transport [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G)),    intent(inout) :: duhdu    !< Partial derivative of uh
                                                          !! with u [H m ~> m2 or kg m-1].
  real,                         intent(in)    :: dt       !< Time increment [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.
  type(ocean_OBC_type), optional, pointer     :: OBC !< Open boundaries control structure.
  ! Local variables
  real :: CFL  ! The CFL number based on the local velocity and grid spacing [nondim]
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h_in.
  real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2].
  integer :: i
  logical :: local_open_BC
 
  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
 
  do i=ish-1,ieh ; if (do_i(i)) then
    ! Set new values of uh and duhdu.
    if (u(i) > 0.0) then
      if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
      else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif
      curv_3 = h_l(i) + h_r(i) - 2.0*h(i)
      uh(i) = g%dy_Cu(i,j) * u(i) * &
          (h_r(i) + cfl * (0.5*(h_l(i) - h_r(i)) + curv_3*(cfl - 1.5)))
      h_marg = h_r(i) + cfl * ((h_l(i) - h_r(i)) + 3.0*curv_3*(cfl - 1.0))
    elseif (u(i) < 0.0) then
      if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
      else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif
      curv_3 = h_l(i+1) + h_r(i+1) - 2.0*h(i+1)
      uh(i) = g%dy_Cu(i,j) * u(i) * &
          (h_l(i+1) + cfl * (0.5*(h_r(i+1)-h_l(i+1)) + curv_3*(cfl - 1.5)))
      h_marg = h_l(i+1) + cfl * ((h_r(i+1)-h_l(i+1)) + 3.0*curv_3*(cfl - 1.0))
    else
      uh(i) = 0.0
      h_marg = 0.5 * (h_l(i+1) + h_r(i))
    endif
    duhdu(i) = g%dy_Cu(i,j) * h_marg * visc_rem(i)
  endif ; enddo
 
  if (local_open_bc) then
    do i=ish-1,ieh ; if (do_i(i)) then
      if (obc%segment(obc%segnum_u(i,j))%open) then
        if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
          uh(i) = g%dy_Cu(i,j) * u(i) * h(i)
          duhdu(i) = g%dy_Cu(i,j) * h(i) * visc_rem(i)
        else
          uh(i) = g%dy_Cu(i,j) * u(i) * h(i+1)
          duhdu(i) = g%dy_Cu(i,j) * h(i+1) * visc_rem(i)
        endif
      endif
    endif ; enddo
  endif

zonal_mass_flux()

subroutine mom_continuity_ppm::zonal_mass_flux ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_in, real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)  uh, real, intent(in)  dt, type(ocean_grid_type), intent(inout)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(continuity_ppm_cs), pointer  CS, type(loop_bounds_type), intent(in)  LB, real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in), optional  uhbt, type(ocean_obc_type), optional, pointer  OBC, real, dimension(szib_(g),szj_(g),szk_(g)), intent(in), optional  visc_rem_u, real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out), optional  u_cor, real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in), optional  uhbt_aux, real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out), optional  u_cor_aux, type(bt_cont_type), optional, pointer  BT_cont  )

private

Calculates the mass or volume fluxes through the zonal faces, and other related quantities.

Parameters

[in,out]	g	Ocean's grid structure.
[in]	gv	Ocean's vertical grid structure.
[in]	u	Zonal velocity [m s-1].
[in]	h_in	Layer thickness used to calculate fluxes [H ~> m or kg m-2].
[out]	uh	Volume flux through zonal faces = u*h*dy
[in]	dt	Time increment [s].
[in]	us	A dimensional unit scaling type
	cs	This module's control structure.
[in]	lb	Loop bounds structure.
	obc	Open boundaries control structure.
[in]	visc_rem_u	The fraction of zonal momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that a layer experiences after viscosity is applied. Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
[in]	uhbt	The summed volume flux through zonal faces
[in]	uhbt_aux	A second set of summed volume fluxes through zonal faces [H m2 s-1 ~> m3 s-1 or kg s-1].
[out]	u_cor	The zonal velocitiess (u with a barotropic correction) that give uhbt as the depth-integrated transport, m s-1.
[out]	u_cor_aux	The zonal velocities (u with a barotropic correction) that give uhbt_aux as the depth-integrated transports [m s-1].
	bt_cont	A structure with elements that describe the effective open face areas as a function of barotropic flow.

Definition at line 233 of file MOM_continuity_PPM.F90.

  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 [m s-1].
  real,  dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h_in !< Layer thickness used to calculate fluxes [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: uh   !< Volume flux through zonal faces = u*h*dy
                                                 !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  real,                    intent(in)    :: dt   !< Time increment [s].
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  type(continuity_PPM_CS), pointer       :: CS   !< This module's control structure.
  type(loop_bounds_type),  intent(in)    :: LB   !< Loop bounds structure.
  type(ocean_OBC_type), &
                 optional, pointer       :: OBC  !< Open boundaries control structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(in)    :: visc_rem_u
                     !< The fraction of zonal momentum originally in a layer that remains after a
                     !! time-step of viscosity, and the fraction of a time-step's worth of a barotropic
                     !! acceleration that a layer experiences after viscosity is applied.
                     !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G),SZJ_(G)), &
                 optional, intent(in)    :: uhbt !< The summed volume flux through zonal faces
                                                 !! [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G),SZJ_(G)), &
                 optional, intent(in)    :: uhbt_aux
                     !< A second set of summed volume fluxes through zonal faces [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: u_cor
                     !< The zonal velocitiess (u with a barotropic correction)
                     !! that give uhbt as the depth-integrated transport, m s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: u_cor_aux
                     !< The zonal velocities (u with a barotropic correction)
                     !! that give uhbt_aux as the depth-integrated transports [m s-1].
  type(BT_cont_type), optional, pointer  :: BT_cont !< A structure with elements that describe the
                     !! effective open face areas as a function of barotropic flow.
 
  ! Local variables
  real, dimension(SZIB_(G),SZK_(G)) :: duhdu ! Partial derivative of uh with u [H m ~> m2 or kg m-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_L, h_R ! Left and right face thicknesses [H ~> m or kg m-2].
  real, dimension(SZIB_(G)) :: &
    du, &      ! Corrective barotropic change in the velocity [m s-1].
    du_min_CFL, & ! Min/max limits on du correction
    du_max_CFL, & ! to avoid CFL violations
    duhdu_tot_0, & ! Summed partial derivative of uh with u [H m ~> m2 or kg m-1].
    uh_tot_0, & ! Summed transport with no barotropic correction [H m2 s-1 ~> m3 s-1 or kg s-1].
    visc_rem_max  ! The column maximum of visc_rem.
  logical, dimension(SZIB_(G)) :: do_I
  real, dimension(SZIB_(G),SZK_(G)) :: &
    visc_rem      ! A 2-D copy of visc_rem_u or an array of 1's.
  real, dimension(SZIB_(G)) :: FAuI  ! A list of sums of zonal face areas [H m ~> m2 or kg m-1].
  real :: FA_u    ! A sum of zonal face areas [H m ~> m2 or kg m-1].
  real :: I_vrm   ! 1.0 / visc_rem_max, nondim.
  real :: CFL_dt  ! The maximum CFL ratio of the adjusted velocities divided by
                  ! the time step [s-1].
  real :: I_dt    ! 1.0 / dt [s-1].
  real :: du_lim  ! The velocity change that give a relative CFL of 1 [m s-1].
  real :: dx_E, dx_W ! Effective x-grid spacings to the east and west [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 => null()
 
  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
    local_open_bc = obc%open_u_BCs_exist_globally
  endif ; endif
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
 
  cfl_dt = cs%CFL_limit_adjust / dt
  i_dt = 1.0 / dt
  if (cs%aggress_adjust) cfl_dt = i_dt
 
  call cpu_clock_begin(id_clock_update)
!$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,CS,h_L,h_in,h_R,G,GV,LB,visc_rem,OBC)
  do k=1,nz
    ! This sets h_L and h_R.
    if (cs%upwind_1st) then
      do j=jsh,jeh ; do i=ish-1,ieh+1
        h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
      enddo ; enddo
    else
      call ppm_reconstruction_x(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
                                2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
    endif
    do i=ish-1,ieh ; visc_rem(i,k) = 1.0 ; enddo
  enddo
  call cpu_clock_end(id_clock_update)
 
  call cpu_clock_begin(id_clock_correct)
!$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,u,h_in,h_L,h_R,use_visc_rem,visc_rem_u,  &
!$OMP                                  uh,dt,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,FAuI,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, obc)
      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, obc=obc)
 
        if (present(u_cor)) then ; do k=1,nz
          do i=ish-1,ieh ; u_cor(i,j,k) = u(i,j,k) + du(i) * visc_rem(i,k) ; enddo
          if (local_specified_bc) then ; do i=ish-1,ieh
            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., obc=obc)
 
        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, us, 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)) faui(i) = 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) .and. &
                (obc%segment(obc%segnum_u(i,j))%specified)) &
              faui(i) = faui(i) + obc%segment(obc%segnum_u(i,j))%normal_trans(i,j,k) / &
                                  obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
          endif ; enddo ; enddo
          do i=ish-1,ieh ; if (do_i(i)) then
            bt_cont%FA_u_W0(i,j) = us%m_to_L*faui(i) ; bt_cont%FA_u_E0(i,j) = us%m_to_L*faui(i)
            bt_cont%FA_u_WW(i,j) = us%m_to_L*faui(i) ; bt_cont%FA_u_EE(i,j) = us%m_to_L*faui(i)
            bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
          endif ; enddo
        endif
      endif ! set_BT_cont
 
    endif ! present(uhbt) or do_aux or set_BT_cont
  enddo ! j-loop
  if (local_open_bc .and. set_bt_cont) then
    do n = 1, obc%number_of_segments
      if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W) then
        i = obc%segment(n)%HI%IsdB
        if (obc%segment(n)%direction == obc_direction_e) then
          do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
            fa_u = 0.0
            do k=1,nz ; fa_u = fa_u + h_in(i,j,k)*us%m_to_L*g%dy_Cu(i,j) ; enddo
            bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
            bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
            bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
          enddo
        else
          do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
            fa_u = 0.0
            do k=1,nz ; fa_u = fa_u + h_in(i+1,j,k)*us%m_to_L*g%dy_Cu(i,j) ; enddo
            bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
            bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
            bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
          enddo
        endif
      endif
    enddo
  endif
  call cpu_clock_end(id_clock_correct)
 
  if  (set_bt_cont) then ; if (allocated(bt_cont%h_u)) then
    if (present(u_cor)) then
      call zonal_face_thickness(u_cor, h_in, h_l, h_r, bt_cont%h_u, dt, g, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
    else
      call zonal_face_thickness(u, h_in, h_l, h_r, bt_cont%h_u, dt, g, lb, &
                                cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
    endif
  endif ; endif