MOM_open_boundary.F90

Go to the documentation of this file.

!> Controls where open boundary conditions are applied
module mom_open_boundary
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_coms,                 only : sum_across_pes
use mom_cpu_clock,            only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
use mom_diag_mediator,        only : diag_ctrl, time_type
use mom_domains,              only : pass_var, pass_vector
use mom_domains,              only : to_all, scalar_pair, cgrid_ne
use mom_error_handler,        only : mom_mesg, mom_error, fatal, warning, is_root_pe
use mom_file_parser,          only : get_param, log_version, param_file_type, log_param
use mom_grid,                 only : ocean_grid_type, hor_index_type
use mom_dyn_horgrid,          only : dyn_horgrid_type
use mom_io,                   only : east_face, north_face
use mom_io,                   only : slasher, read_data, field_size, single_file
use mom_io,                   only : vardesc, query_vardesc, var_desc
use mom_restart,              only : register_restart_field, query_initialized, mom_restart_cs
use mom_obsolete_params,      only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char
use mom_string_functions,     only : extract_word, remove_spaces
use mom_tracer_registry,      only : tracer_type, tracer_registry_type, tracer_name_lookup
use time_interp_external_mod, only : init_external_field, time_interp_external
use time_interp_external_mod, only : time_interp_external_init
use mom_remapping,            only : remappingschemesdoc, remappingdefaultscheme, remapping_cs
use mom_remapping,            only : initialize_remapping, remapping_core_h, end_remapping
use mom_regridding,           only : regridding_cs
use mom_unit_scaling,         only : unit_scale_type
use mom_variables,            only : thermo_var_ptrs
use mom_verticalgrid,         only : verticalgrid_type
 
implicit none ; private
 
#include <MOM_memory.h>
 
public open_boundary_apply_normal_flow
public open_boundary_config
public open_boundary_init
public open_boundary_query
public open_boundary_end
public open_boundary_impose_normal_slope
public open_boundary_impose_land_mask
public radiation_open_bdry_conds
public set_tracer_data
public update_obc_segment_data
public open_boundary_test_extern_uv
public open_boundary_test_extern_h
public open_boundary_zero_normal_flow
public register_obc, obc_registry_init
public register_file_obc, file_obc_end
public segment_tracer_registry_init
public segment_tracer_registry_end
public register_segment_tracer
public register_temp_salt_segments
public fill_temp_salt_segments
public open_boundary_register_restarts
public update_segment_tracer_reservoirs
 
integer, parameter, public :: obc_none = 0      !< Indicates the use of no open boundary
integer, parameter, public :: obc_simple = 1    !< Indicates the use of a simple inflow open boundary
integer, parameter, public :: obc_wall = 2      !< Indicates the use of a closed wall
integer, parameter, public :: obc_flather =  3  !< Indicates the use of a Flather open boundary
integer, parameter, public :: obc_radiation = 4 !< Indicates the use of a radiation open boundary
integer, parameter, public :: obc_direction_n = 100 !< Indicates the boundary is an effective northern boundary
integer, parameter, public :: obc_direction_s = 200 !< Indicates the boundary is an effective southern boundary
integer, parameter, public :: obc_direction_e = 300 !< Indicates the boundary is an effective eastern boundary
integer, parameter, public :: obc_direction_w = 400 !< Indicates the boundary is an effective western boundary
integer, parameter         :: max_obc_fields = 100  !< Maximum number of data fields needed for OBC segments
 
!> Open boundary segment data from files (mostly).
type, public :: obc_segment_data_type
  integer :: fid                                !< handle from FMS associated with segment data on disk
  integer :: fid_dz                             !< handle from FMS associated with segment thicknesses on disk
  character(len=8)                :: name       !< a name identifier for the segment data
  real, pointer, dimension(:,:,:) :: buffer_src=>null() !< buffer for segment data located at cell faces
                                                !! and on the original vertical grid
  integer                         :: nk_src     !< Number of vertical levels in the source data
  real, dimension(:,:,:), pointer :: dz_src=>null()     !< vertical grid cell spacing of the incoming segment data [m]
  real, dimension(:,:,:), pointer :: buffer_dst=>null() !< buffer src data remapped to the target vertical grid
  real, dimension(:,:), pointer   :: bt_vel=>null()     !< barotropic velocity [L T-1 ~> m s-1]
  real                            :: value              !< constant value if fid is equal to -1
end type obc_segment_data_type
 
!> Tracer on OBC segment data structure, for putting into a segment tracer registry.
type, public :: obc_segment_tracer_type
  real, dimension(:,:,:), pointer :: t          => null()  !< tracer concentration array
  real                            :: obc_inflow_conc = 0.0 !< tracer concentration for generic inflows
  character(len=32)               :: name                  !< tracer name used for error messages
  type(tracer_type), pointer      :: tr         => null()  !< metadata describing the tracer
  real, dimension(:,:,:), pointer :: tres       => null()  !< tracer reservoir array
  logical                         :: is_initialized        !< reservoir values have been set when True
end type obc_segment_tracer_type
 
!> Registry type for tracers on segments
type, public :: segment_tracer_registry_type
  integer                       :: ntseg = 0         !< number of registered tracer segments
  type(obc_segment_tracer_type) :: tr(max_fields_)   !< array of registered tracers
  logical                       :: locked = .false.  !< New tracers may be registered if locked=.false.
                                                     !! When locked=.true.,no more tracers can be registered.
                                                     !! Not sure who should lock it or when...
end type segment_tracer_registry_type
 
!> Open boundary segment data structure.
type, public :: obc_segment_type
  logical :: flather        !< If true, applies Flather + Chapman radiation of barotropic gravity waves.
  logical :: radiation      !< If true, 1D Orlanksi radiation boundary conditions are applied.
                            !! If False, a gradient condition is applied.
  logical :: radiation_tan  !< If true, 1D Orlanksi radiation boundary conditions are applied to
                            !! tangential flows.
  logical :: radiation_grad !< If true, 1D Orlanksi radiation boundary conditions are applied to
                            !! dudv and dvdx.
  logical :: oblique        !< Oblique waves supported at radiation boundary.
  logical :: oblique_tan    !< If true, 2D radiation boundary conditions are applied to
                            !! tangential flows.
  logical :: oblique_grad   !< If true, 2D radiation boundary conditions are applied to
                            !! dudv and dvdx.
  logical :: nudged         !< Optional supplement to radiation boundary.
  logical :: nudged_tan     !< Optional supplement to nudge tangential velocity.
  logical :: nudged_grad    !< Optional supplement to nudge normal gradient of tangential velocity.
  logical :: specified      !< Boundary normal velocity fixed to external value.
  logical :: specified_tan  !< Boundary tangential velocity fixed to external value.
  logical :: open           !< Boundary is open for continuity solver.
  logical :: gradient       !< Zero gradient at boundary.
  logical :: values_needed  !< Whether or not any external OBC fields are needed.
  logical :: u_values_needed!< Whether or not external u OBC fields are needed.
  logical :: v_values_needed!< Whether or not external v OBC fields are needed.
  logical :: t_values_needed!< Whether or not external T OBC fields are needed.
  logical :: s_values_needed!< Whether or not external S OBC fields are needed.
  logical :: z_values_needed!< Whether or not external zeta OBC fields are needed.
  logical :: g_values_needed!< Whether or not external gradient OBC fields are needed.
  integer :: direction      !< Boundary faces one of the four directions.
  logical :: is_n_or_s      !< True if the OB is facing North or South and exists on this PE.
  logical :: is_e_or_w      !< True if the OB is facing East or West and exists on this PE.
  logical :: is_e_or_w_2    !< True if the OB is facing East or West anywhere.
  type(obc_segment_data_type), pointer, dimension(:) :: field=>null()   !<  OBC data
  integer :: num_fields     !< number of OBC data fields (e.g. u_normal,u_parallel and eta for Flather)
  character(len=32), pointer, dimension(:) :: field_names=>null() !< field names for this segment
  integer :: is_obc         !< i-indices of boundary segment.
  integer :: ie_obc         !< i-indices of boundary segment.
  integer :: js_obc         !< j-indices of boundary segment.
  integer :: je_obc         !< j-indices of boundary segment.
  real :: velocity_nudging_timescale_in  !< Nudging timescale on inflow [T ~> s].
  real :: velocity_nudging_timescale_out !< Nudging timescale on outflow [T ~> s].
  logical :: on_pe          !< true if segment is located in the computational domain
  logical :: temp_segment_data_exists !< true if temperature data arrays are present
  logical :: salt_segment_data_exists !< true if salinity data arrays are present
  real, pointer, dimension(:,:)   :: cg=>null()     !< The external gravity wave speed [L T-1 ~> m s-1]
                                                    !! at OBC-points.
  real, pointer, dimension(:,:)   :: htot=>null()   !< The total column thickness [H ~> m or kg m-2] at OBC-points.
  real, pointer, dimension(:,:,:) :: h=>null()      !< The cell thickness [H ~> m or kg m-2] at OBC-points.
  real, pointer, dimension(:,:,:) :: normal_vel=>null()     !< The layer velocity normal to the OB
                                                            !! segment [L T-1 ~> m s-1].
  real, pointer, dimension(:,:,:) :: tangential_vel=>null() !< The layer velocity tangential to the
                                                            !! OB segment [L T-1 ~> m s-1].
  real, pointer, dimension(:,:,:) :: tangential_grad=>null() !< The gradient of the velocity tangential
                                                            !! to the OB segment [T-1 ~> s-1].
  real, pointer, dimension(:,:,:) :: normal_trans=>null()   !< The layer transport normal to the OB
                                                            !! segment [H L2 T-1 ~> m3 s-1].
  real, pointer, dimension(:,:)   :: normal_vel_bt=>null()  !< The barotropic velocity normal to
                                                            !! the OB segment [L T-1 ~> m s-1].
  real, pointer, dimension(:,:)   :: eta=>null()            !< The sea-surface elevation along the segment [m].
  real, pointer, dimension(:,:,:) :: grad_normal=>null()    !< The gradient of the normal flow along the
                                                            !! segment times the grid spacing [L T-1 ~> m s-1]
  real, pointer, dimension(:,:,:) :: grad_tan=>null()       !< The gradient of the tangential flow along the
                                                            !! segment times the grid spacing [L T-1 ~> m s-1]
  real, pointer, dimension(:,:,:) :: grad_gradient=>null()  !< The gradient of the gradient of tangential flow along
                                                            !! the segment times the grid spacing [T-1 ~> s-1]
  real, pointer, dimension(:,:,:) :: rx_norm_rad=>null()    !< The previous normal phase speed use for EW radiation
                                                            !! OBC, in grid points per timestep [nondim]
  real, pointer, dimension(:,:,:) :: ry_norm_rad=>null()    !< The previous normal phase speed use for NS radiation
                                                            !! OBC, in grid points per timestep [nondim]
  real, pointer, dimension(:,:,:) :: rx_norm_obl=>null()    !< The previous normal radiation coefficient for EW
                                                            !! oblique OBCs [L2 T-2 ~> m2 s-2]
  real, pointer, dimension(:,:,:) :: ry_norm_obl=>null()    !< The previous normal radiation coefficient for NS
                                                            !! oblique OBCs [L2 T-2 ~> m2 s-2]
  real, pointer, dimension(:,:,:) :: cff_normal=>null()     !< The denominator for oblique radiation
                                                            !! for normal velocity [L2 T-2 ~> m2 s-2]
  real, pointer, dimension(:,:,:) :: nudged_normal_vel=>null() !< The layer velocity normal to the OB segment
                                                            !! that values should be nudged towards [L T-1 ~> m s-1].
  real, pointer, dimension(:,:,:) :: nudged_tangential_vel=>null() !< The layer velocity tangential to the OB segment
                                                            !! that values should be nudged towards [L T-1 ~> m s-1].
  real, pointer, dimension(:,:,:) :: nudged_tangential_grad=>null() !< The layer dvdx or dudy towards which nudging
                                                            !! can occur [T-1 ~> s-1].
  type(segment_tracer_registry_type), pointer  :: tr_reg=> null()!< A pointer to the tracer registry for the segment.
  type(hor_index_type) :: hi !< Horizontal index ranges
  real :: tr_invlscale_out                                  !< An effective inverse length scale for restoring
                                                            !! the tracer concentration in a ficticious
                                                            !! reservior towards interior values when flow
                                                            !! is exiting the domain [L-1 ~> m-1]
  real :: tr_invlscale_in                                   !< An effective inverse length scale for restoring
                                                            !! the tracer concentration towards an externally
                                                            !! imposed value when flow is entering [L-1 ~> m-1]
end type obc_segment_type
 
!> Open-boundary data
type, public :: ocean_obc_type
  integer :: number_of_segments = 0                   !< The number of open-boundary segments.
  integer :: ke = 0                                   !< The number of model layers
  logical :: open_u_bcs_exist_globally = .false.      !< True if any zonal velocity points
                                                      !! in the global domain use open BCs.
  logical :: open_v_bcs_exist_globally = .false.      !< True if any meridional velocity points
                                                      !! in the global domain use open BCs.
  logical :: flather_u_bcs_exist_globally = .false.   !< True if any zonal velocity points
                                                      !! in the global domain use Flather BCs.
  logical :: flather_v_bcs_exist_globally = .false.   !< True if any meridional velocity points
                                                      !! in the global domain use Flather BCs.
  logical :: oblique_bcs_exist_globally = .false.     !< True if any velocity points
                                                      !! in the global domain use oblique BCs.
  logical :: nudged_u_bcs_exist_globally = .false.    !< True if any velocity points in the
                                                      !! global domain use nudged BCs.
  logical :: nudged_v_bcs_exist_globally = .false.    !< True if any velocity points in the
                                                      !! global domain use nudged BCs.
  logical :: specified_u_bcs_exist_globally = .false. !< True if any zonal velocity points
                                                      !! in the global domain use specified BCs.
  logical :: specified_v_bcs_exist_globally = .false. !< True if any meridional velocity points
                                                      !! in the global domain use specified BCs.
  logical :: radiation_bcs_exist_globally = .false.   !< True if radiations BCs are in use anywhere.
  logical :: user_bcs_set_globally = .false.          !< True if any OBC_USER_CONFIG is set
                                                      !! for input from user directory.
  logical :: update_obc = .false.                     !< Is OBC data time-dependent
  logical :: needs_io_for_data = .false.              !< Is any i/o needed for OBCs
  logical :: zero_vorticity = .false.                 !< If True, sets relative vorticity to zero on open boundaries.
  logical :: freeslip_vorticity = .false.             !< If True, sets normal gradient of tangential velocity to zero
                                                      !! in the relative vorticity on open boundaries.
  logical :: computed_vorticity = .false.             !< If True, uses external data for tangential velocity
                                                      !! in the relative vorticity on open boundaries.
  logical :: specified_vorticity = .false.            !< If True, uses external data for tangential velocity
                                                      !! gradients in the relative vorticity on open boundaries.
  logical :: zero_strain = .false.                    !< If True, sets strain to zero on open boundaries.
  logical :: freeslip_strain = .false.                !< If True, sets normal gradient of tangential velocity to zero
                                                      !! in the strain on open boundaries.
  logical :: computed_strain = .false.                !< If True, uses external data for tangential velocity to compute
                                                      !! normal gradient in the strain on open boundaries.
  logical :: specified_strain = .false.               !< If True, uses external data for tangential velocity gradients
                                                      !! to compute strain on open boundaries.
  logical :: zero_biharmonic = .false.                !< If True, zeros the Laplacian of flow on open boundaries for
                                                      !! use in the biharmonic viscosity term.
  logical :: brushcutter_mode = .false.               !< If True, read data on supergrid.
  real :: g_earth                                     !< The gravitational acceleration [m s-2].
  logical, pointer, dimension(:) :: &
                   tracer_x_reservoirs_used => null() !< Dimensioned by the number of tracers, set globally,
                                                      !! true for those with x reservoirs (needed for restarts).
  logical, pointer, dimension(:) :: &
                   tracer_y_reservoirs_used => null() !< Dimensioned by the number of tracers, set globally,
                                                      !! true for those with y reservoirs (needed for restarts).
  integer                       :: ntr = 0            !< number of tracers
  ! Properties of the segments used.
  type(obc_segment_type), pointer, dimension(:) :: &
    segment => null()   !< List of segment objects.
  ! Which segment object describes the current point.
  integer, pointer, dimension(:,:) :: &
    segnum_u => null(), &   !< Segment number of u-points.
    segnum_v => null()      !< Segment number of v-points.
 
  ! The following parameters are used in the baroclinic radiation code:
  real :: gamma_uv !< The relative weighting for the baroclinic radiation
                   !! velocities (or speed of characteristics) at the
                   !! new time level (1) or the running mean (0) for velocities.
                   !! Valid values range from 0 to 1, with a default of 0.3.
  real :: rx_max   !< The maximum magnitude of the baroclinic radiation
                   !! velocity (or speed of characteristics) [m s-1].  The
                   !! default value is 10 m s-1.
  logical :: obc_pe !< Is there an open boundary on this tile?
  type(remapping_cs),      pointer :: remap_cs   !< ALE remapping control structure for segments only
  type(obc_registry_type), pointer :: obc_reg => null()  !< Registry type for boundaries
  real, pointer, dimension(:,:,:) :: &
    rx_normal => null(), & !< Array storage for normal phase speed for EW radiation OBCs in units of
                           !! grid points per timestep [nondim]
    ry_normal => null(), & !< Array storage for normal phase speed for NS radiation OBCs in units of
                           !! grid points per timestep [nondim]
    rx_oblique => null(), & !< Array storage for oblique boundary condition restarts [L2 T-2 ~> m2 s-2]
    ry_oblique => null(), & !< Array storage for oblique boundary condition restarts [L2 T-2 ~> m2 s-2]
    cff_normal => null()   !< Array storage for oblique boundary condition restarts [L2 T-2 ~> m2 s-2]
  real, pointer, dimension(:,:,:,:) :: &
    tres_x => null(), & !< Array storage of tracer reservoirs for restarts [conc L ~> conc m]
    tres_y => null()    !< Array storage of tracer reservoirs for restarts [conc L ~> conc m]
  real :: silly_h  !< A silly value of thickness outside of the domain that can be used to test
                   !! the independence of the OBCs to this external data [H ~> m or kg m-2].
  real :: silly_u  !< A silly value of velocity outside of the domain that can be used to test
                   !! the independence of the OBCs to this external data [L T-1 ~> m s-1].
end type ocean_obc_type
 
!> Control structure for open boundaries that read from files.
!! Probably lots to update here.
type, public :: file_obc_cs ; private
  real :: tide_flow = 3.0e6         !< Placeholder for now...
end type file_obc_cs
 
!> Type to carry something (what??) for the OBC registry.
type, public :: obc_struct_type
  character(len=32)               :: name             !< OBC name used for error messages
end type obc_struct_type
 
!> Type to carry basic OBC information needed for updating values.
type, public :: obc_registry_type
  integer               :: nobc = 0          !< number of registered open boundary types.
  type(obc_struct_type) :: ob(max_fields_)   !< array of registered boundary types.
  logical               :: locked = .false.  !< New OBC types may be registered if locked=.false.
                                             !! When locked=.true.,no more boundaries can be registered.
end type obc_registry_type
 
integer :: id_clock_pass !< A CPU time clock
 
character(len=40)  :: mdl = "MOM_open_boundary" !< This module's name.
! This include declares and sets the variable "version".
#include "version_variable.h"
 
contains
 
!> Enables OBC module and reads configuration parameters
!> This routine is called from MOM_initialize_fixed which
!> occurs before the initialization of the vertical coordinate
!> and ALE_init.  Therefore segment data are not fully initialized
!> here. The remainder of the segment data are initialized in a
!> later call to update_open_boundary_data
 
subroutine open_boundary_config(G, US, param_file, OBC)
  type(dyn_horgrid_type),  intent(inout) :: g   !< Ocean grid structure
  type(unit_scale_type),   intent(in)    :: us  !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< Parameter file handle
  type(ocean_obc_type),    pointer       :: obc !< Open boundary control structure
  ! Local variables
  integer :: l ! For looping over segments
  logical :: debug_obc, debug, mask_outside, reentrant_x, reentrant_y
  character(len=15) :: segment_param_str ! The run-time parameter name for each segment
  character(len=100) :: segment_str      ! The contents (rhs) for parameter "segment_param_str"
  character(len=200) :: config1          ! String for OBC_USER_CONFIG
  real               :: lscale_in, lscale_out ! parameters controlling tracer values at the boundaries [L ~> m]
  allocate(obc)
 
  call log_version(param_file, mdl, version, &
                 "Controls where open boundaries are located, what kind of boundary condition "//&
                 "to impose, and what data to apply, if any.")
  call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", obc%number_of_segments, &
                 "The number of open boundary segments.", &
                 default=0)
  call get_param(param_file, mdl, "G_EARTH", obc%g_Earth, &
                 "The gravitational acceleration of the Earth.", &
                 units="m s-2", default = 9.80)
  call get_param(param_file, mdl, "OBC_USER_CONFIG", config1, &
                 "A string that sets how the open boundary conditions are "//&
                 " configured: \n", default="none", do_not_log=.true.)
  call get_param(param_file, mdl, "NK", obc%ke, &
                 "The number of model layers", default=0, do_not_log=.true.)
 
  if (config1 /= "none" .and. config1 /= "dyed_obcs") obc%user_BCs_set_globally = .true.
 
  if (obc%number_of_segments > 0) then
    call get_param(param_file, mdl, "OBC_ZERO_VORTICITY", obc%zero_vorticity, &
         "If true, sets relative vorticity to zero on open boundaries.", &
         default=.false.)
    call get_param(param_file, mdl, "OBC_FREESLIP_VORTICITY", obc%freeslip_vorticity, &
         "If true, sets the normal gradient of tangential velocity to "//&
         "zero in the relative vorticity on open boundaries. This cannot "//&
         "be true if another OBC_XXX_VORTICITY option is True.", default=.true.)
    call get_param(param_file, mdl, "OBC_COMPUTED_VORTICITY", obc%computed_vorticity, &
         "If true, uses the external values of tangential velocity "//&
         "in the relative vorticity on open boundaries. This cannot "//&
         "be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
    call get_param(param_file, mdl, "OBC_SPECIFIED_VORTICITY", obc%specified_vorticity, &
         "If true, uses the external values of tangential velocity "//&
         "in the relative vorticity on open boundaries. This cannot "//&
         "be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
    if ((obc%zero_vorticity .and. obc%freeslip_vorticity) .or.  &
        (obc%zero_vorticity .and. obc%computed_vorticity) .or.  &
        (obc%zero_vorticity .and. obc%specified_vorticity) .or.  &
        (obc%freeslip_vorticity .and. obc%computed_vorticity) .or.  &
        (obc%freeslip_vorticity .and. obc%specified_vorticity) .or.  &
        (obc%computed_vorticity .and. obc%specified_vorticity))  &
         call mom_error(fatal, "MOM_open_boundary.F90, open_boundary_config:\n"//&
         "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"//&
         "and OBC_IMPORTED_VORTICITY can be True at once.")
    call get_param(param_file, mdl, "OBC_ZERO_STRAIN", obc%zero_strain, &
         "If true, sets the strain used in the stress tensor to zero on open boundaries.", &
         default=.false.)
    call get_param(param_file, mdl, "OBC_FREESLIP_STRAIN", obc%freeslip_strain, &
         "If true, sets the normal gradient of tangential velocity to "//&
         "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
         "be true if another OBC_XXX_STRAIN option is True.", default=.true.)
    call get_param(param_file, mdl, "OBC_COMPUTED_STRAIN", obc%computed_strain, &
         "If true, sets the normal gradient of tangential velocity to "//&
         "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
         "be true if another OBC_XXX_STRAIN option is True.", default=.false.)
    call get_param(param_file, mdl, "OBC_SPECIFIED_STRAIN", obc%specified_strain, &
         "If true, sets the normal gradient of tangential velocity to "//&
         "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
         "be true if another OBC_XXX_STRAIN option is True.", default=.false.)
    if ((obc%zero_strain .and. obc%freeslip_strain) .or.  &
        (obc%zero_strain .and. obc%computed_strain) .or.  &
        (obc%zero_strain .and. obc%specified_strain) .or.  &
        (obc%freeslip_strain .and. obc%computed_strain) .or.  &
        (obc%freeslip_strain .and. obc%specified_strain) .or.  &
        (obc%computed_strain .and. obc%specified_strain))  &
         call mom_error(fatal, "MOM_open_boundary.F90, open_boundary_config: \n"//&
         "Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN, OBC_COMPUTED_STRAIN \n"//&
         "and OBC_IMPORTED_STRAIN can be True at once.")
    call get_param(param_file, mdl, "OBC_ZERO_BIHARMONIC", obc%zero_biharmonic, &
         "If true, zeros the Laplacian of flow on open boundaries in the biharmonic "//&
         "viscosity term.", default=.false.)
    call get_param(param_file, mdl, "MASK_OUTSIDE_OBCS", mask_outside, &
         "If true, set the areas outside open boundaries to be land.", &
         default=.false.)
    call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
    call get_param(param_file, mdl, "DEBUG_OBC", debug_obc, default=.false.)
    if (debug_obc .or. debug) &
      call log_param(param_file, mdl, "DEBUG_OBC", debug_obc, &
                 "If true, do additional calls to help debug the performance "//&
                 "of the open boundary condition code.", default=.false., &
                 debuggingparam=.true.)
 
    call get_param(param_file, mdl, "OBC_SILLY_THICK", obc%silly_h, &
                 "A silly value of thicknesses used outside of open boundary "//&
                 "conditions for debugging.", units="m", default=0.0, scale=us%m_to_Z, &
                 do_not_log=.not.debug_obc, debuggingparam=.true.)
    call get_param(param_file, mdl, "OBC_SILLY_VEL", obc%silly_u, &
                 "A silly value of velocities used outside of open boundary "//&
                 "conditions for debugging.", units="m/s", default=0.0, scale=us%m_s_to_L_T, &
                 do_not_log=.not.debug_obc, debuggingparam=.true.)
    reentrant_x = .false.
    call get_param(param_file, mdl, "REENTRANT_X", reentrant_x, default=.true.)
    reentrant_y = .false.
    call get_param(param_file, mdl, "REENTRANT_Y", reentrant_y, default=.false.)
 
    ! Allocate everything
    ! Note the 0-segment is needed when %segnum_u/v(:,:) = 0
    allocate(obc%segment(0:obc%number_of_segments))
    do l=0,obc%number_of_segments
      obc%segment(l)%Flather = .false.
      obc%segment(l)%radiation = .false.
      obc%segment(l)%radiation_tan = .false.
      obc%segment(l)%radiation_grad = .false.
      obc%segment(l)%oblique = .false.
      obc%segment(l)%oblique_tan = .false.
      obc%segment(l)%oblique_grad = .false.
      obc%segment(l)%nudged = .false.
      obc%segment(l)%nudged_tan = .false.
      obc%segment(l)%nudged_grad = .false.
      obc%segment(l)%specified = .false.
      obc%segment(l)%specified_tan = .false.
      obc%segment(l)%open = .false.
      obc%segment(l)%gradient = .false.
      obc%segment(l)%values_needed = .false.
      obc%segment(l)%u_values_needed = .false.
      obc%segment(l)%v_values_needed = .false.
      obc%segment(l)%t_values_needed = .false.
      obc%segment(l)%s_values_needed = .false.
      obc%segment(l)%z_values_needed = .false.
      obc%segment(l)%g_values_needed = .false.
      obc%segment(l)%direction = obc_none
      obc%segment(l)%is_N_or_S = .false.
      obc%segment(l)%is_E_or_W = .false.
      obc%segment(l)%is_E_or_W_2 = .false.
      obc%segment(l)%Velocity_nudging_timescale_in = 0.0
      obc%segment(l)%Velocity_nudging_timescale_out = 0.0
      obc%segment(l)%num_fields = 0
    enddo
    allocate(obc%segnum_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; obc%segnum_u(:,:) = obc_none
    allocate(obc%segnum_v(g%isd:g%ied,g%JsdB:g%JedB)) ; obc%segnum_v(:,:) = obc_none
 
    do l = 1, obc%number_of_segments
      write(segment_param_str(1:15),"('OBC_SEGMENT_',i3.3)") l
      call get_param(param_file, mdl, segment_param_str, segment_str, &
           "Documentation needs to be dynamic?????", &
           fail_if_missing=.true.)
      segment_str = remove_spaces(segment_str)
      if (segment_str(1:2) == 'I=') then
        call setup_u_point_obc(obc, g, us, segment_str, l, param_file, reentrant_y)
      elseif (segment_str(1:2) == 'J=') then
        call setup_v_point_obc(obc, g, us, segment_str, l, param_file, reentrant_x)
      else
        call mom_error(fatal, "MOM_open_boundary.F90, open_boundary_config: "//&
             "Unable to interpret "//segment_param_str//" = "//trim(segment_str))
      endif
    enddo
 
    !    if (open_boundary_query(OBC, needs_ext_seg_data=.true.)) &
    call initialize_segment_data(g, obc, param_file)
 
    if (open_boundary_query(obc, apply_open_obc=.true.)) then
      call get_param(param_file, mdl, "OBC_RADIATION_MAX", obc%rx_max, &
                   "The maximum magnitude of the baroclinic radiation "//&
                   "velocity (or speed of characteristics).  This is only "//&
                   "used if one of the open boundary segments is using Orlanski.", &
                   units="m s-1", default=10.0)
      call get_param(param_file, mdl, "OBC_RAD_VEL_WT", obc%gamma_uv, &
                   "The relative weighting for the baroclinic radiation "//&
                   "velocities (or speed of characteristics) at the new "//&
                   "time level (1) or the running mean (0) for velocities. "//&
                   "Valid values range from 0 to 1. This is only used if "//&
                   "one of the open boundary segments is using Orlanski.", &
                   units="nondim", default=0.3)
    endif
 
    lscale_in = 0.
    lscale_out = 0.
    if (open_boundary_query(obc, apply_open_obc=.true.)) then
      call get_param(param_file, mdl, "OBC_TRACER_RESERVOIR_LENGTH_SCALE_OUT ", lscale_out, &
                 "An effective length scale for restoring the tracer concentration "//&
                 "at the boundaries to externally imposed values when the flow "//&
                 "is exiting the domain.", units="m", default=0.0, scale=us%m_to_L)
 
      call get_param(param_file, mdl, "OBC_TRACER_RESERVOIR_LENGTH_SCALE_IN ", lscale_in, &
                 "An effective length scale for restoring the tracer concentration "//&
                 "at the boundaries to values from the interior when the flow "//&
                 "is entering the domain.", units="m", default=0.0, scale=us%m_to_L)
    endif
 
    if (mask_outside) call mask_outside_obcs(g, us, param_file, obc)
 
    ! All tracers are using the same restoring length scale for now, but we may want to make this
    ! tracer-specific in the future for example, in cases where certain tracers are poorly constrained
    ! by data while others are well constrained - MJH.
    do l = 1, obc%number_of_segments
      obc%segment(l)%Tr_InvLscale_in = 0.0
      if (lscale_in>0.) obc%segment(l)%Tr_InvLscale_in =  1.0/lscale_in
      obc%segment(l)%Tr_InvLscale_out = 0.0
      if (lscale_out>0.) obc%segment(l)%Tr_InvLscale_out =  1.0/lscale_out
    enddo
 
  endif ! OBC%number_of_segments > 0
 
    ! Safety check
  if ((obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally) .and. &
       .not.g%symmetric ) call mom_error(fatal, &
       "MOM_open_boundary, open_boundary_config: "//&
       "Symmetric memory must be used when using Flather OBCs.")
 
  if (.not.(obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
              obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) then
    ! No open boundaries have been requested
    call open_boundary_dealloc(obc)
  else
    ! Need this for ocean_only mode boundary interpolation.
    call time_interp_external_init()
  endif
 
end subroutine open_boundary_config
 
!> Allocate space for reading OBC data from files. It sets up the required vertical
!! remapping. In the process, it does funky stuff with the MPI processes.
subroutine initialize_segment_data(G, OBC, PF)
  use mpp_mod, only : mpp_pe, mpp_set_current_pelist, mpp_get_current_pelist,mpp_npes
 
  type(dyn_horgrid_type), intent(in)    :: G   !< Ocean grid structure
  type(ocean_obc_type),   intent(inout) :: OBC !< Open boundary control structure
  type(param_file_type),  intent(in)    :: PF  !< Parameter file handle
 
  integer :: n,m,num_fields
  character(len=256) :: segstr, filename
  character(len=20)  :: segnam, suffix
  character(len=32)  :: varnam, fieldname
  real               :: value
  character(len=32), dimension(MAX_OBC_FIELDS) :: fields  ! segment field names
  character(len=128) :: inputdir
  type(obc_segment_type), pointer :: segment => null() ! pointer to segment type list
  character(len=32)  :: remappingScheme
  character(len=256) :: mesg    ! Message for error messages.
  logical :: check_reconstruction, check_remapping, force_bounds_in_subcell
  integer, dimension(4) :: siz,siz2
  integer :: is, ie, js, je
  integer :: isd, ied, jsd, jed
  integer :: IsdB, IedB, JsdB, JedB
  integer, dimension(:), allocatable :: saved_pelist
  integer :: current_pe
  integer, dimension(1) :: single_pelist
  !will be able to dynamically switch between sub-sampling refined grid data or model grid
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
 
  ! There is a problem with the order of the OBC initialization
  ! with respect to ALE_init. Currently handling this by copying the
  ! param file so that I can use it later in step_MOM in order to finish
  ! initializing segments on the first step.
 
  call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
 
  call get_param(pf, mdl, "REMAPPING_SCHEME", remappingscheme, &
          "This sets the reconstruction scheme used "//&
          "for vertical remapping for all variables. "//&
          "It can be one of the following schemes: \n"//&
          trim(remappingschemesdoc), default=remappingdefaultscheme,do_not_log=.true.)
  call get_param(pf, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
          "If true, cell-by-cell reconstructions are checked for "//&
          "consistency and if non-monotonicity or an inconsistency is "//&
          "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
  call get_param(pf, mdl, "FATAL_CHECK_REMAPPING", check_remapping, &
          "If true, the results of remapping are checked for "//&
          "conservation and new extrema and if an inconsistency is "//&
          "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
  call get_param(pf, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
          "If true, the values on the intermediate grid used for remapping "//&
          "are forced to be bounded, which might not be the case due to "//&
          "round off.", default=.false.,do_not_log=.true.)
    call get_param(pf, mdl, "BRUSHCUTTER_MODE", obc%brushcutter_mode, &
         "If true, read external OBC data on the supergrid.", &
         default=.false.)
 
  allocate(obc%remap_CS)
  call initialize_remapping(obc%remap_CS, remappingscheme, boundary_extrapolation = .false., &
       check_reconstruction=check_reconstruction, &
       check_remapping=check_remapping, force_bounds_in_subcell=force_bounds_in_subcell)
 
  if (obc%user_BCs_set_globally) return
 
  ! Try this here just for the documentation. It is repeated below.
  do n=1, obc%number_of_segments
    segment => obc%segment(n)
    write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n
    call get_param(pf, mdl, segnam, segstr, 'OBC segment docs')
  enddo
 
  !< temporarily disable communication in order to read segment data independently
 
  allocate(saved_pelist(0:mpp_npes()-1))
  call mpp_get_current_pelist(saved_pelist)
  current_pe = mpp_pe()
  single_pelist(1) = current_pe
  call mpp_set_current_pelist(single_pelist)
 
  do n=1, obc%number_of_segments
    segment => obc%segment(n)
    if (.not. segment%values_needed) cycle
 
    write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n
    write(suffix,"('_segment_',i3.3)") n
    ! needs documentation !!  Yet, unsafe for now, causes grief for
    ! MOM_parameter_docs in circle_obcs on two processes.
!   call get_param(PF, mdl, segnam, segstr, 'xyz')
    ! Clear out any old values
    segstr = ''
    call get_param(pf, mdl, segnam, segstr)
    if (segstr == '') then
      write(mesg,'("No OBC_SEGMENT_XXX_DATA string for OBC segment ",I3)') n
      call mom_error(fatal, mesg)
    endif
 
    call parse_segment_data_str(trim(segstr), fields=fields, num_fields=num_fields)
    if (num_fields == 0) then
      call mom_mesg('initialize_segment_data: num_fields = 0')
      cycle ! cycle to next segment
    endif
 
    allocate(segment%field(num_fields))
    segment%num_fields = num_fields
 
    segment%temp_segment_data_exists=.false.
    segment%salt_segment_data_exists=.false.
!!
! CODE HERE FOR OTHER OPTIONS (CLAMPED, NUDGED,..)
!!
 
    isd = segment%HI%isd ; ied = segment%HI%ied
    jsd = segment%HI%jsd ; jed = segment%HI%jed
    isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
    jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
 
    do m=1,num_fields
      call parse_segment_data_str(trim(segstr), var=trim(fields(m)), value=value, filenam=filename, fieldnam=fieldname)
      if (trim(filename) /= 'none') then
        obc%update_OBC = .true. ! Data is assumed to be time-dependent if we are reading from file
        obc%needs_IO_for_data = .true. ! At least one segment is using I/O for OBC data
!       segment%values_needed = .true. ! Indicates that i/o will be needed for this segment
        segment%field(m)%name = trim(fields(m))
        if (segment%field(m)%name == 'TEMP') then
           segment%temp_segment_data_exists=.true.
           segment%t_values_needed = .false.
        endif
        if (segment%field(m)%name == 'SALT') then
           segment%salt_segment_data_exists=.true.
           segment%s_values_needed = .false.
        endif
        filename = trim(inputdir)//trim(filename)
        fieldname = trim(fieldname)//trim(suffix)
        call field_size(filename,fieldname,siz,no_domain=.true.)
!       if (siz(4) == 1) segment%values_needed = .false.
        if (segment%on_pe) then
          if (obc%brushcutter_mode .and. (modulo(siz(1),2) == 0 .or. modulo(siz(2),2) == 0)) then
            call mom_error(fatal,'segment data are not on the supergrid')
          endif
          siz2(1)=1
 
          if (siz(1)>1) then
            if (obc%brushcutter_mode) then
              siz2(1)=(siz(1)-1)/2
            else
              siz2(1)=siz(1)
            endif
          endif
          siz2(2)=1
          if (siz(2)>1) then
            if (obc%brushcutter_mode) then
              siz2(2)=(siz(2)-1)/2
            else
              siz2(2)=siz(2)
            endif
          endif
          siz2(3)=siz(3)
 
          if (segment%is_E_or_W) then
            if (segment%field(m)%name == 'V') then
              allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
              segment%v_values_needed = .false.
            else if (segment%field(m)%name == 'DVDX') then
              allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
              segment%g_values_needed = .false.
            else
              allocate(segment%field(m)%buffer_src(isdb:iedb,jsd:jed,siz2(3)))
              if (segment%field(m)%name == 'U') then
                segment%u_values_needed = .false.
              else if (segment%field(m)%name == 'SSH') then
                segment%z_values_needed = .false.
              else if (segment%field(m)%name == 'TEMP') then
                segment%t_values_needed = .false.
              else if (segment%field(m)%name == 'SALT') then
                segment%s_values_needed = .false.
              endif
            endif
          else
            if (segment%field(m)%name == 'U') then
              allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
              segment%u_values_needed = .false.
            else if (segment%field(m)%name == 'DUDY') then
              allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
              segment%g_values_needed = .false.
            else
              allocate(segment%field(m)%buffer_src(isd:ied,jsdb:jedb,siz2(3)))
              if (segment%field(m)%name == 'V') then
                segment%v_values_needed = .false.
              else if (segment%field(m)%name == 'SSH') then
                segment%z_values_needed = .false.
              else if (segment%field(m)%name == 'TEMP') then
                segment%t_values_needed = .false.
              else if (segment%field(m)%name == 'SALT') then
                segment%s_values_needed = .false.
              endif
            endif
          endif
          segment%field(m)%buffer_src(:,:,:)=0.0
          segment%field(m)%fid = init_external_field(trim(filename),&
               trim(fieldname),ignore_axis_atts=.true.,threading=single_file)
          if (siz(3) > 1) then
            fieldname = 'dz_'//trim(fieldname)
            call field_size(filename,fieldname,siz,no_domain=.true.)
            if (segment%is_E_or_W) then
              if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
                allocate(segment%field(m)%dz_src(isdb:iedb,jsdb:jedb,siz(3)))
              else
                allocate(segment%field(m)%dz_src(isdb:iedb,jsd:jed,siz(3)))
              endif
            else
              if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
                allocate(segment%field(m)%dz_src(isdb:iedb,jsdb:jedb,siz(3)))
              else
                allocate(segment%field(m)%dz_src(isd:ied,jsdb:jedb,siz(3)))
              endif
            endif
            segment%field(m)%dz_src(:,:,:)=0.0
            segment%field(m)%nk_src=siz(3)
            segment%field(m)%fid_dz = init_external_field(trim(filename),trim(fieldname),&
                       ignore_axis_atts=.true.,threading=single_file)
          else
            segment%field(m)%nk_src=1
          endif
        endif
      else
        segment%field(m)%fid = -1
        segment%field(m)%value = value
        segment%field(m)%name = trim(fields(m))
        if (segment%field(m)%name == 'U') then
          segment%u_values_needed = .false.
        elseif (segment%field(m)%name == 'V') then
          segment%v_values_needed = .false.
        elseif (segment%field(m)%name == 'SSH') then
          segment%z_values_needed = .false.
        elseif (segment%field(m)%name == 'TEMP') then
          segment%t_values_needed = .false.
        elseif (segment%field(m)%name == 'SALT') then
          segment%s_values_needed = .false.
        elseif (segment%field(m)%name == 'DVDX' .or. segment%field(m)%name == 'DUDY') then
          segment%g_values_needed = .false.
        endif
      endif
    enddo
    if (segment%u_values_needed .or. segment%v_values_needed .or. &
        segment%t_values_needed .or. segment%s_values_needed .or. &
        segment%z_values_needed .or. segment%g_values_needed) then
      write(mesg,'("Values needed for OBC segment ",I3)') n
      call mom_error(fatal, mesg)
    endif
  enddo
 
  call mpp_set_current_pelist(saved_pelist)
 
end subroutine initialize_segment_data
 
!> Define indices for segment and store in hor_index_type
!> using global segment bounds corresponding to q-points
subroutine setup_segment_indices(G, seg, Is_obc, Ie_obc, Js_obc, Je_obc)
  type(dyn_horgrid_type), intent(in) :: G !< grid type
  type(obc_segment_type), intent(inout) :: seg  !< Open boundary segment
  integer, intent(in) :: Is_obc !< Q-point global i-index of start of segment
  integer, intent(in) :: Ie_obc !< Q-point global i-index of end of segment
  integer, intent(in) :: Js_obc !< Q-point global j-index of start of segment
  integer, intent(in) :: Je_obc !< Q-point global j-index of end of segment
  ! Local variables
  integer :: Isg,Ieg,Jsg,Jeg
 
  ! Isg, Ieg will be I*_obc in global space
  if (ie_obc<is_obc) then
    isg=ie_obc;ieg=is_obc
  else
    isg=is_obc;ieg=ie_obc
  endif
  if (je_obc<js_obc) then
    jsg=je_obc;jeg=js_obc
  else
    jsg=js_obc;jeg=je_obc
  endif
 
  ! Global space I*_obc but sorted
  seg%HI%IsgB = isg ; seg%HI%IegB = ieg
  seg%HI%isg = isg+1 ; seg%HI%ieg = ieg
  seg%HI%JsgB = jsg ; seg%HI%JegB = jeg
  seg%HI%jsg = jsg+1 ; seg%HI%Jeg = jeg
 
  ! Move into local index space
  isg = isg - g%idg_offset
  jsg = jsg - g%jdg_offset
  ieg = ieg - g%idg_offset
  jeg = jeg - g%jdg_offset
 
  ! This is the i-extent of the segment on this PE.
  ! The values are nonsense if the segment is not on this PE.
  seg%HI%IsdB = min( max(isg, g%HI%IsdB), g%HI%IedB)
  seg%HI%IedB = min( max(ieg, g%HI%IsdB), g%HI%IedB)
  seg%HI%isd = min( max(isg+1, g%HI%isd), g%HI%ied)
  seg%HI%ied = min( max(ieg, g%HI%isd), g%HI%ied)
  seg%HI%IscB = min( max(isg, g%HI%IscB), g%HI%IecB)
  seg%HI%IecB = min( max(ieg, g%HI%IscB), g%HI%IecB)
  seg%HI%isc = min( max(isg+1, g%HI%isc), g%HI%iec)
  seg%HI%iec = min( max(ieg, g%HI%isc), g%HI%iec)
 
  ! This is the j-extent of the segment on this PE.
  ! The values are nonsense if the segment is not on this PE.
  seg%HI%JsdB = min( max(jsg, g%HI%JsdB), g%HI%JedB)
  seg%HI%JedB = min( max(jeg, g%HI%JsdB), g%HI%JedB)
  seg%HI%jsd = min( max(jsg+1, g%HI%jsd), g%HI%jed)
  seg%HI%jed = min( max(jeg, g%HI%jsd), g%HI%jed)
  seg%HI%JscB = min( max(jsg, g%HI%JscB), g%HI%JecB)
  seg%HI%JecB = min( max(jeg, g%HI%JscB), g%HI%JecB)
  seg%HI%jsc = min( max(jsg+1, g%HI%jsc), g%HI%jec)
  seg%HI%jec = min( max(jeg, g%HI%jsc), g%HI%jec)
 
end subroutine setup_segment_indices
 
!> Parse an OBC_SEGMENT_%%% string starting with "I=" and configure placement and type of OBC accordingly
subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y)
  type(ocean_obc_type),    pointer    :: OBC !< Open boundary control structure
  type(dyn_horgrid_type),  intent(in) :: G   !< Ocean grid structure
  type(unit_scale_type),   intent(in) :: US  !< A dimensional unit scaling type
  character(len=*),        intent(in) :: segment_str !< A string in form of "I=%,J=%:%,string"
  integer,                 intent(in) :: l_seg !< which segment is this?
  type(param_file_type), intent(in)   :: PF  !< Parameter file handle
  logical, intent(in)                 :: reentrant_y !< is the domain reentrant in y?
  ! Local variables
  integer :: I_obc, Js_obc, Je_obc ! Position of segment in global index space
  integer :: j, a_loop
  character(len=32) :: action_str(8)
  character(len=128) :: segment_param_str
  real, allocatable, dimension(:)  :: tnudge
  ! This returns the global indices for the segment
  call parse_segment_str(g%ieg, g%jeg, segment_str, i_obc, js_obc, je_obc, action_str, reentrant_y)
 
  call setup_segment_indices(g, obc%segment(l_seg),i_obc,i_obc,js_obc,je_obc)
 
  i_obc = i_obc - g%idg_offset ! Convert to local tile indices on this tile
  js_obc = js_obc - g%jdg_offset ! Convert to local tile indices on this tile
  je_obc = je_obc - g%jdg_offset ! Convert to local tile indices on this tile
 
  if (je_obc>js_obc) then
    obc%segment(l_seg)%direction = obc_direction_e
  elseif (je_obc<js_obc) then
    obc%segment(l_seg)%direction = obc_direction_w
    j=js_obc;js_obc=je_obc;je_obc=j
  endif
 
  obc%segment(l_seg)%on_pe = .false.
 
  do a_loop = 1,8 ! up to 8 options available
    if (len_trim(action_str(a_loop)) == 0) then
      cycle
    elseif (trim(action_str(a_loop)) == 'FLATHER') then
      obc%segment(l_seg)%Flather = .true.
      obc%segment(l_seg)%open = .true.
      obc%Flather_u_BCs_exist_globally = .true.
      obc%open_u_BCs_exist_globally = .true.
      obc%segment%z_values_needed = .true.
      obc%segment%u_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'ORLANSKI') then
      obc%segment(l_seg)%radiation = .true.
      obc%segment(l_seg)%open = .true.
      obc%open_u_BCs_exist_globally = .true.
      obc%radiation_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'ORLANSKI_TAN') then
      obc%segment(l_seg)%radiation = .true.
      obc%segment(l_seg)%radiation_tan = .true.
      obc%radiation_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'ORLANSKI_GRAD') then
      obc%segment(l_seg)%radiation = .true.
      obc%segment(l_seg)%radiation_grad = .true.
    elseif (trim(action_str(a_loop)) == 'OBLIQUE') then
      obc%segment(l_seg)%oblique = .true.
      obc%segment(l_seg)%open = .true.
      obc%oblique_BCs_exist_globally = .true.
      obc%open_u_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'OBLIQUE_TAN') then
      obc%segment(l_seg)%oblique = .true.
      obc%segment(l_seg)%oblique_tan = .true.
      obc%oblique_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'OBLIQUE_GRAD') then
      obc%segment(l_seg)%oblique = .true.
      obc%segment(l_seg)%oblique_grad = .true.
    elseif (trim(action_str(a_loop)) == 'NUDGED') then
      obc%segment(l_seg)%nudged = .true.
      obc%nudged_u_BCs_exist_globally = .true.
      obc%segment%u_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'NUDGED_TAN') then
      obc%segment(l_seg)%nudged_tan = .true.
      obc%nudged_u_BCs_exist_globally = .true.
      obc%segment%v_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'NUDGED_GRAD') then
      obc%segment(l_seg)%nudged_grad = .true.
      obc%segment%g_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'GRADIENT') then
      obc%segment(l_seg)%gradient = .true.
      obc%segment(l_seg)%open = .true.
      obc%open_u_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'SIMPLE') then
      obc%segment(l_seg)%specified = .true.
      obc%specified_u_BCs_exist_globally = .true. ! This avoids deallocation
      obc%segment%u_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then
      obc%segment(l_seg)%specified_tan = .true.
    else
      call mom_error(fatal, "MOM_open_boundary.F90, setup_u_point_obc: "//&
                     "String '"//trim(action_str(a_loop))//"' not understood.")
    endif
    if (obc%segment(l_seg)%nudged .or. obc%segment(l_seg)%nudged_tan) then
      write(segment_param_str(1:43),"('OBC_SEGMENT_',i3.3,'_VELOCITY_NUDGING_TIMESCALES')") l_seg
      allocate(tnudge(2))
      call get_param(pf, mdl, segment_param_str(1:43), tnudge, &
                     "Timescales in days for nudging along a segment, "//&
                     "for inflow, then outflow. Setting both to zero should "//&
                     "behave like SIMPLE obcs for the baroclinic velocities.", &
                     fail_if_missing=.true., default=0., units="days", scale=86400.0*us%s_to_T)
      obc%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1)
      obc%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2)
      deallocate(tnudge)
    endif
 
  enddo ! a_loop
 
  obc%segment(l_seg)%is_E_or_W_2 = .true.
 
  if (i_obc<=g%HI%IsdB+1 .or. i_obc>=g%HI%IedB-1) return ! Boundary is not on tile
  if (je_obc<=g%HI%JsdB .or. js_obc>=g%HI%JedB) return ! Segment is not on tile
 
  obc%segment(l_seg)%on_pe = .true.
  obc%segment(l_seg)%is_E_or_W = .true.
 
  do j=g%HI%jsd, g%HI%jed
    if (j>js_obc .and. j<=je_obc) then
      obc%segnum_u(i_obc,j) = l_seg
    endif
  enddo
  obc%segment(l_seg)%Is_obc = i_obc
  obc%segment(l_seg)%Ie_obc = i_obc
  obc%segment(l_seg)%Js_obc = js_obc
  obc%segment(l_seg)%Je_obc = je_obc
  call allocate_obc_segment_data(obc, obc%segment(l_seg))
 
  if (obc%segment(l_seg)%oblique .and.  obc%segment(l_seg)%radiation) &
         call mom_error(fatal, "MOM_open_boundary.F90, setup_u_point_obc: \n"//&
         "Orlanski and Oblique OBC options cannot be used together on one segment.")
 
  if (obc%segment(l_seg)%u_values_needed .or. obc%segment(l_seg)%v_values_needed .or. &
      obc%segment(l_seg)%t_values_needed .or. obc%segment(l_seg)%s_values_needed .or. &
      obc%segment(l_seg)%z_values_needed .or. obc%segment(l_seg)%g_values_needed) &
    obc%segment(l_seg)%values_needed = .true.
end subroutine setup_u_point_obc
 
!> Parse an OBC_SEGMENT_%%% string starting with "J=" and configure placement and type of OBC accordingly
subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x)
  type(ocean_obc_type),    pointer    :: OBC !< Open boundary control structure
  type(dyn_horgrid_type),  intent(in) :: G   !< Ocean grid structure
  type(unit_scale_type),   intent(in) :: US  !< A dimensional unit scaling type
  character(len=*),        intent(in) :: segment_str !< A string in form of "J=%,I=%:%,string"
  integer,                 intent(in) :: l_seg !< which segment is this?
  type(param_file_type),   intent(in) :: PF  !< Parameter file handle
  logical, intent(in)                 :: reentrant_x !< is the domain reentrant in x?
  ! Local variables
  integer :: J_obc, Is_obc, Ie_obc ! Position of segment in global index space
  integer :: i, a_loop
  character(len=32) :: action_str(8)
  character(len=128) :: segment_param_str
  real, allocatable, dimension(:)  :: tnudge
 
  ! This returns the global indices for the segment
  call parse_segment_str(g%ieg, g%jeg, segment_str, j_obc, is_obc, ie_obc, action_str, reentrant_x)
 
  call setup_segment_indices(g, obc%segment(l_seg),is_obc,ie_obc,j_obc,j_obc)
 
  j_obc = j_obc - g%jdg_offset ! Convert to local tile indices on this tile
  is_obc = is_obc - g%idg_offset ! Convert to local tile indices on this tile
  ie_obc = ie_obc - g%idg_offset ! Convert to local tile indices on this tile
 
  if (ie_obc>is_obc) then
     obc%segment(l_seg)%direction = obc_direction_s
  elseif (ie_obc<is_obc) then
     obc%segment(l_seg)%direction = obc_direction_n
     i=is_obc;is_obc=ie_obc;ie_obc=i
  endif
 
  obc%segment(l_seg)%on_pe = .false.
 
  do a_loop = 1,8
    if (len_trim(action_str(a_loop)) == 0) then
      cycle
    elseif (trim(action_str(a_loop)) == 'FLATHER') then
      obc%segment(l_seg)%Flather = .true.
      obc%segment(l_seg)%open = .true.
      obc%Flather_v_BCs_exist_globally = .true.
      obc%open_v_BCs_exist_globally = .true.
      obc%segment%z_values_needed = .true.
      obc%segment%v_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'ORLANSKI') then
      obc%segment(l_seg)%radiation = .true.
      obc%segment(l_seg)%open = .true.
      obc%open_v_BCs_exist_globally = .true.
      obc%radiation_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'ORLANSKI_TAN') then
      obc%segment(l_seg)%radiation = .true.
      obc%segment(l_seg)%radiation_tan = .true.
      obc%radiation_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'ORLANSKI_GRAD') then
      obc%segment(l_seg)%radiation = .true.
      obc%segment(l_seg)%radiation_grad = .true.
    elseif (trim(action_str(a_loop)) == 'OBLIQUE') then
      obc%segment(l_seg)%oblique = .true.
      obc%segment(l_seg)%open = .true.
      obc%oblique_BCs_exist_globally = .true.
      obc%open_v_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'OBLIQUE_TAN') then
      obc%segment(l_seg)%oblique = .true.
      obc%segment(l_seg)%oblique_tan = .true.
      obc%oblique_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'OBLIQUE_GRAD') then
      obc%segment(l_seg)%oblique = .true.
      obc%segment(l_seg)%oblique_grad = .true.
    elseif (trim(action_str(a_loop)) == 'NUDGED') then
      obc%segment(l_seg)%nudged = .true.
      obc%nudged_v_BCs_exist_globally = .true.
      obc%segment%v_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'NUDGED_TAN') then
      obc%segment(l_seg)%nudged_tan = .true.
      obc%nudged_v_BCs_exist_globally = .true.
      obc%segment%u_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'NUDGED_GRAD') then
      obc%segment(l_seg)%nudged_grad = .true.
      obc%segment%g_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'GRADIENT') then
      obc%segment(l_seg)%gradient = .true.
      obc%segment(l_seg)%open = .true.
      obc%open_v_BCs_exist_globally = .true.
    elseif (trim(action_str(a_loop)) == 'SIMPLE') then
      obc%segment(l_seg)%specified = .true.
      obc%specified_v_BCs_exist_globally = .true. ! This avoids deallocation
      obc%segment%v_values_needed = .true.
    elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then
      obc%segment(l_seg)%specified_tan = .true.
    else
      call mom_error(fatal, "MOM_open_boundary.F90, setup_v_point_obc: "//&
                     "String '"//trim(action_str(a_loop))//"' not understood.")
    endif
    if (obc%segment(l_seg)%nudged .or. obc%segment(l_seg)%nudged_tan) then
      write(segment_param_str(1:43),"('OBC_SEGMENT_',i3.3,'_VELOCITY_NUDGING_TIMESCALES')") l_seg
      allocate(tnudge(2))
      call get_param(pf, mdl, segment_param_str(1:43), tnudge, &
                     "Timescales in days for nudging along a segment, "//&
                     "for inflow, then outflow. Setting both to zero should "//&
                     "behave like SIMPLE obcs for the baroclinic velocities.", &
                     fail_if_missing=.true., default=0., units="days", scale=86400.0*us%s_to_T)
      obc%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1)
      obc%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2)
      deallocate(tnudge)
    endif
 
  enddo ! a_loop
 
  if (j_obc<=g%HI%JsdB+1 .or. j_obc>=g%HI%JedB-1) return ! Boundary is not on tile
  if (ie_obc<=g%HI%IsdB .or. is_obc>=g%HI%IedB) return ! Segment is not on tile
 
  obc%segment(l_seg)%on_pe = .true.
  obc%segment(l_seg)%is_N_or_S = .true.
 
  do i=g%HI%isd, g%HI%ied
    if (i>is_obc .and. i<=ie_obc) then
      obc%segnum_v(i,j_obc) = l_seg
    endif
  enddo
  obc%segment(l_seg)%Is_obc = is_obc
  obc%segment(l_seg)%Ie_obc = ie_obc
  obc%segment(l_seg)%Js_obc = j_obc
  obc%segment(l_seg)%Je_obc = j_obc
  call allocate_obc_segment_data(obc, obc%segment(l_seg))
 
  if (obc%segment(l_seg)%oblique .and.  obc%segment(l_seg)%radiation) &
         call mom_error(fatal, "MOM_open_boundary.F90, setup_v_point_obc: \n"//&
         "Orlanski and Oblique OBC options cannot be used together on one segment.")
 
  if (obc%segment(l_seg)%u_values_needed .or. obc%segment(l_seg)%v_values_needed .or. &
      obc%segment(l_seg)%t_values_needed .or. obc%segment(l_seg)%s_values_needed .or. &
      obc%segment(l_seg)%z_values_needed .or. obc%segment(l_seg)%g_values_needed) &
    obc%segment(l_seg)%values_needed = .true.
end subroutine setup_v_point_obc
 
!> Parse an OBC_SEGMENT_%%% string
subroutine parse_segment_str(ni_global, nj_global, segment_str, l, m, n, action_str, reentrant)
  integer,          intent(in)  :: ni_global !< Number of h-points in zonal direction
  integer,          intent(in)  :: nj_global !< Number of h-points in meridional direction
  character(len=*), intent(in)  :: segment_str !< A string in form of "I=l,J=m:n,string" or "J=l,I=m,n,string"
  integer,          intent(out) :: l !< The value of I=l, if segment_str begins with I=l, or the value of J=l
  integer,          intent(out) :: m !< The value of J=m, if segment_str begins with I=, or the value of I=m
  integer,          intent(out) :: n !< The value of J=n, if segment_str begins with I=, or the value of I=n
  character(len=*), intent(out) :: action_str(:) !< The "string" part of segment_str
  logical,          intent(in)  :: reentrant !< is domain reentrant in relevant direction?
  ! Local variables
  character(len=24) :: word1, word2, m_word, n_word !< Words delineated by commas in a string in form of
                                                    !! "I=%,J=%:%,string"
  integer :: l_max !< Either ni_global or nj_global, depending on whether segment_str begins with "I=" or "J="
  integer :: mn_max !< Either nj_global or ni_global, depending on whether segment_str begins with "I=" or "J="
  integer :: j
  integer, parameter :: halo = 10
 
  ! Process first word which will started with either 'I=' or 'J='
  word1 = extract_word(segment_str,',',1)
  word2 = extract_word(segment_str,',',2)
  if (word1(1:2)=='I=') then
    l_max = ni_global
    mn_max = nj_global
    if (.not. (word2(1:2)=='J=')) call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
                     "Second word of string '"//trim(segment_str)//"' must start with 'J='.")
  elseif (word1(1:2)=='J=') then ! Note that the file_parser uniformaly expands "=" to " = "
    l_max = nj_global
    mn_max = ni_global
    if (.not. (word2(1:2)=='I=')) call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
                     "Second word of string '"//trim(segment_str)//"' must start with 'I='.")
  else
    call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str"//&
                   "String '"//segment_str//"' must start with 'I=' or 'J='.")
  endif
 
  ! Read l
  l = interpret_int_expr( word1(3:24), l_max )
  if (l<0 .or. l>l_max) then
    call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
                   "First value from string '"//trim(segment_str)//"' is outside of the physical domain.")
  endif
 
  ! Read m
  m_word = extract_word(word2(3:24),':',1)
  m = interpret_int_expr( m_word, mn_max )
  if (reentrant) then
    if (m<-halo .or. m>mn_max+halo) then
      call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
                     "Beginning of range in string '"//trim(segment_str)//"' is outside of the physical domain.")
    endif
  else
    if (m<-1 .or. m>mn_max+1) then
      call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
                     "Beginning of range in string '"//trim(segment_str)//"' is outside of the physical domain.")
    endif
  endif
 
  ! Read n
  n_word = extract_word(word2(3:24),':',2)
  n = interpret_int_expr( n_word, mn_max )
  if (reentrant) then
    if (n<-halo .or. n>mn_max+halo) then
      call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
                     "End of range in string '"//trim(segment_str)//"' is outside of the physical domain.")
    endif
  else
    if (n<-1 .or. n>mn_max+1) then
      call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
                     "End of range in string '"//trim(segment_str)//"' is outside of the physical domain.")
    endif
  endif
 
  if (abs(n-m)==0) then
    call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
                   "Range in string '"//trim(segment_str)//"' must span one cell.")
  endif
 
  ! Type of open boundary condition
  do j = 1, size(action_str)
    action_str(j) = extract_word(segment_str,',',2+j)
  enddo
 
  contains
 
  ! Returns integer value interpreted from string in form of %I, N or N+-%I
  integer function interpret_int_expr(string, imax)
    character(len=*), intent(in) :: string !< Integer in form or %I, N or N-%I
    integer,          intent(in) :: imax !< Value to replace 'N' with
    ! Local variables
    integer slen
 
    slen = len_trim(string)
    if (slen==0) call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str"//&
                                "Parsed string was empty!")
    if (len_trim(string)==1 .and. string(1:1)=='N') then
      interpret_int_expr = imax
    elseif (string(1:1)=='N') then
      if (string(2:2)=='+') then
        read(string(3:slen),*,err=911) interpret_int_expr
        interpret_int_expr = imax + interpret_int_expr
      elseif (string(2:2)=='-') then
        read(string(3:slen),*,err=911) interpret_int_expr
        interpret_int_expr = imax - interpret_int_expr
      endif
    else
      read(string(1:slen),*,err=911) interpret_int_expr
    endif
    return
    911 call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str"//&
                       "Problem reading value from string '"//trim(string)//"'.")
  end function interpret_int_expr
end subroutine parse_segment_str
 
!> Parse an OBC_SEGMENT_%%%_DATA string
 subroutine parse_segment_data_str(segment_str, var, value, filenam, fieldnam, fields, num_fields, debug )
   character(len=*),           intent(in)   :: segment_str !< A string in form of
                                                          !! "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..."
   character(len=*), optional, intent(in)   :: var        !< The name of the variable for which parameters are needed
   character(len=*), optional, intent(out)  :: filenam    !< The name of the input file if using "file" method
   character(len=*), optional, intent(out)  :: fieldnam   !< The name of the variable in the input file if using
                                                          !! "file" method
   real,             optional, intent(out)  :: value      !< A constant value if using the "value" method
   character(len=*), dimension(MAX_OBC_FIELDS), &
                     optional, intent(out)  :: fields     !< List of fieldnames for each segment
   integer, optional, intent(out)           :: num_fields !< The number of fields in the segment data
   logical, optional, intent(in)            :: debug      !< If present and true, write verbose debugging messages
   ! Local variables
   character(len=128) :: word1, word2, word3, method
   integer :: lword, nfields, n, m
   logical :: continue,dbg
   character(len=32), dimension(MAX_OBC_FIELDS) :: flds
 
   nfields=0
   continue=.true.
   dbg=.false.
   if (PRESENT(debug)) dbg=debug
 
   do while (continue)
      word1 = extract_word(segment_str,',',nfields+1)
      if (trim(word1) == '') exit
      nfields=nfields+1
      word2 = extract_word(word1,'=',1)
      flds(nfields) = trim(word2)
   enddo
 
   if (PRESENT(fields)) then
     do n=1,nfields
       fields(n) = flds(n)
     enddo
   endif
 
   if (PRESENT(num_fields)) then
      num_fields=nfields
      return
   endif
 
   m=0
   if (PRESENT(var)) then
     do n=1,nfields
       if (trim(var)==trim(flds(n))) then
          m=n
          exit
       endif
     enddo
     if (m==0) then
        call abort()
     endif
 
    ! Process first word which will start with the fieldname
     word3 = extract_word(segment_str,',',m)
     word1 = extract_word(word3,':',1)
!     if (trim(word1) == '') exit
     word2 = extract_word(word1,'=',1)
     if (trim(word2) == trim(var)) then
        method=trim(extract_word(word1,'=',2))
        lword=len_trim(method)
        if (method(lword-3:lword) == 'file') then
           ! raise an error id filename/fieldname not in argument list
           word1 = extract_word(word3,':',2)
           filenam = extract_word(word1,'(',1)
           fieldnam = extract_word(word1,'(',2)
           lword=len_trim(fieldnam)
           fieldnam = fieldnam(1:lword-1)  ! remove trailing parenth
           value=-999.
        elseif (method(lword-4:lword) == 'value') then
           filenam = 'none'
           fieldnam = 'none'
           word1 = extract_word(word3,':',2)
           lword=len_trim(word1)
           read(word1(1:lword),*,end=986,err=987) value
        endif
      endif
    endif
 
   return
 986 call mom_error(fatal,'End of record while parsing segment data specification! '//trim(segment_str))
 987 call mom_error(fatal,'Error while parsing segment data specification! '//trim(segment_str))
 
 end subroutine parse_segment_data_str
 
 
!> Parse all the OBC_SEGMENT_%%%_DATA strings again
!! to see which need tracer reservoirs (all pes need to know).
 subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature)
  type(ocean_obc_type),   intent(inout) :: OBC !< Open boundary control structure
  type(param_file_type),  intent(in)    :: PF  !< Parameter file handle
  logical,                intent(in) :: use_temperature !< If true, T and S are used
 
  ! Local variables
  integer :: n,m,num_fields
  character(len=256) :: segstr, filename
  character(len=20)  :: segnam, suffix
  character(len=32)  :: varnam, fieldname
  real               :: value
  character(len=32), dimension(MAX_OBC_FIELDS) :: fields  ! segment field names
  type(obc_segment_type), pointer :: segment => null() ! pointer to segment type list
  character(len=256) :: mesg    ! Message for error messages.
 
  do n=1, obc%number_of_segments
    segment => obc%segment(n)
    write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n
    write(suffix,"('_segment_',i3.3)") n
    ! Clear out any old values
    segstr = ''
    call get_param(pf, mdl, segnam, segstr)
    if (segstr == '') cycle
 
    call parse_segment_data_str(trim(segstr), fields=fields, num_fields=num_fields)
    if (num_fields == 0) cycle
 
    ! At this point, just search for TEMP and SALT as tracers 1 and 2.
    do m=1,num_fields
      call parse_segment_data_str(trim(segstr), var=trim(fields(m)), value=value, filenam=filename, fieldnam=fieldname)
      if (trim(filename) /= 'none') then
        if (fields(m) == 'TEMP') then
          if (segment%is_E_or_W_2) then
            obc%tracer_x_reservoirs_used(1) = .true.
          else
            obc%tracer_y_reservoirs_used(1) = .true.
          endif
        endif
        if (fields(m) == 'SALT') then
          if (segment%is_E_or_W_2) then
            obc%tracer_x_reservoirs_used(2) = .true.
          else
            obc%tracer_y_reservoirs_used(2) = .true.
          endif
        endif
      endif
    enddo
    ! Alternately, set first two to true if use_temperature is true
    if (use_temperature) then
      if (segment%is_E_or_W_2) then
        obc%tracer_x_reservoirs_used(1) = .true.
        obc%tracer_x_reservoirs_used(2) = .true.
      else
        obc%tracer_y_reservoirs_used(1) = .true.
        obc%tracer_y_reservoirs_used(2) = .true.
      endif
    endif
  enddo
 
  return
 
end subroutine parse_for_tracer_reservoirs
 
!> Parse an OBC_SEGMENT_%%%_PARAMS string
 subroutine parse_segment_param_real(segment_str, var, param_value, debug )
   character(len=*),  intent(in)  :: segment_str !< A string in form of
                                                 !! "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..."
   character(len=*),  intent(in)  :: var         !< The name of the variable for which parameters are needed
   real,              intent(out) :: param_value !< The value of the parameter
   logical, optional, intent(in)  :: debug       !< If present and true, write verbose debugging messages
   ! Local variables
   character(len=128) :: word1, word2, word3, method
   integer :: lword, nfields, n, m
   logical :: continue,dbg
   character(len=32), dimension(MAX_OBC_FIELDS) :: flds
 
   nfields=0
   continue=.true.
   dbg=.false.
   if (PRESENT(debug)) dbg=debug
 
   do while (continue)
      word1 = extract_word(segment_str,',',nfields+1)
      if (trim(word1) == '') exit
      nfields=nfields+1
      word2 = extract_word(word1,'=',1)
      flds(nfields) = trim(word2)
   enddo
 
   ! if (PRESENT(fields)) then
   !   do n=1,nfields
   !     fields(n) = flds(n)
   !   enddo
   ! endif
 
   ! if (PRESENT(num_fields)) then
   !    num_fields=nfields
   !    return
   ! endif
 
   m=0
!   if (PRESENT(var)) then
     do n=1,nfields
       if (trim(var)==trim(flds(n))) then
          m=n
          exit
       endif
     enddo
     if (m==0) then
        call abort()
     endif
 
    ! Process first word which will start with the fieldname
     word3 = extract_word(segment_str,',',m)
!     word1 = extract_word(word3,':',1)
!     if (trim(word1) == '') exit
     word2 = extract_word(word1,'=',1)
     if (trim(word2) == trim(var)) then
        method=trim(extract_word(word1,'=',2))
        lword=len_trim(method)
        read(method(1:lword),*,err=987) param_value
        ! if (method(lword-3:lword) == 'file') then
        !    ! raise an error id filename/fieldname not in argument list
        !    word1 = extract_word(word3,':',2)
        !    filenam = extract_word(word1,'(',1)
        !    fieldnam = extract_word(word1,'(',2)
        !    lword=len_trim(fieldnam)
        !    fieldnam = fieldnam(1:lword-1)  ! remove trailing parenth
        !    value=-999.
        ! elseif (method(lword-4:lword) == 'value') then
        !    filenam = 'none'
        !    fieldnam = 'none'
        !    word1 = extract_word(word3,':',2)
        !    lword=len_trim(word1)
        !    read(word1(1:lword),*,end=986,err=987) value
        ! endif
      endif
!    endif
 
   return
 986 call mom_error(fatal,'End of record while parsing segment data specification! '//trim(segment_str))
 987 call mom_error(fatal,'Error while parsing segment parameter specification! '//trim(segment_str))
 
 end subroutine parse_segment_param_real
 
!> Initialize open boundary control structure and do any necessary rescaling of OBC
!! fields that have been read from a restart file.
subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CSp)
  type(ocean_grid_type),   intent(in) :: g   !< Ocean grid structure
  type(verticalgrid_type), intent(in) :: gv  !< Container for vertical grid information
  type(unit_scale_type),   intent(in) :: us  !< A dimensional unit scaling type
  type(param_file_type),   intent(in) :: param_file !< Parameter file handle
  type(ocean_obc_type),    pointer    :: obc !< Open boundary control structure
  type(mom_restart_cs),    pointer    :: restart_csp !< Restart structure, data intent(inout)
 
  ! Local variables
  real :: vel2_rescale ! A rescaling factor for squared velocities from the representation in
                       ! a restart file to the internal representation in this run.
  integer :: i, j, k, isd, ied, jsd, jed, nz
  integer :: isdb, iedb, jsdb, jedb
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nz = gv%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
 
  if (.not.associated(obc)) return
 
  id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=clock_routine)
 
  ! The rx_normal and ry_normal arrays used with radiation OBCs are currently in units of grid
  ! points per timestep, but if this were to be corrected to [L T-1 ~> m s-1] or [T-1 ~> s-1] to
  ! permit timesteps to change between calls to the OBC code, the following would be needed:
!  if ( OBC%radiation_BCs_exist_globally .and. (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. &
!       ((US%m_to_L * US%s_to_T_restart) /= (US%m_to_L_restart * US%s_to_T)) ) then
!    vel_rescale = (US%m_to_L * US%s_to_T_restart) /  (US%m_to_L_restart * US%s_to_T)
!    if (query_initialized(OBC%rx_normal, "rx_normal", restart_CSp)) then
!      do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB
!        OBC%rx_normal(I,j,k) = vel_rescale * OBC%rx_normal(I,j,k)
!      enddo ; enddo ; enddo
!    endif
!    if (query_initialized(OBC%ry_normal, "ry_normal", restart_CSp)) then
!      do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied
!        OBC%ry_normal(i,J,k) = vel_rescale * OBC%ry_normal(i,J,k)
!      enddo ; enddo ; enddo
!    endif
!  endif
 
  ! The oblique boundary condition terms have units of [L2 T-2 ~> m2 s-2] and may need to be rescaled.
  if ( obc%oblique_BCs_exist_globally .and. (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
       ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) ) then
    vel2_rescale = (us%m_to_L * us%s_to_T_restart)**2 /  (us%m_to_L_restart * us%s_to_T)**2
    if (query_initialized(obc%rx_oblique, "rx_oblique", restart_csp)) then
      do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
        obc%rx_oblique(i,j,k) = vel2_rescale * obc%rx_oblique(i,j,k)
      enddo ; enddo ; enddo
    endif
    if (query_initialized(obc%ry_oblique, "ry_oblique", restart_csp)) then
      do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
        obc%ry_oblique(i,j,k) = vel2_rescale * obc%ry_oblique(i,j,k)
      enddo ; enddo ; enddo
    endif
    if (query_initialized(obc%cff_normal, "cff_normal", restart_csp)) then
      do k=1,nz ; do j=jsdb,jedb ; do i=isdb,iedb
        obc%cff_normal(i,j,k) = vel2_rescale * obc%cff_normal(i,j,k)
      enddo ; enddo ; enddo
    endif
   endif
 
end subroutine open_boundary_init
 
logical function open_boundary_query(OBC, apply_open_OBC, apply_specified_OBC, apply_Flather_OBC, &
                                     apply_nudged_OBC, needs_ext_seg_data)
  type(ocean_obc_type), pointer    :: obc !< Open boundary control structure
  logical, optional,    intent(in) :: apply_open_obc      !< Returns True if open_*_BCs_exist_globally is true
  logical, optional,    intent(in) :: apply_specified_obc !< Returns True if specified_*_BCs_exist_globally is true
  logical, optional,    intent(in) :: apply_flather_obc   !< Returns True if Flather_*_BCs_exist_globally is true
  logical, optional,    intent(in) :: apply_nudged_obc    !< Returns True if nudged_*_BCs_exist_globally is true
  logical, optional,    intent(in) :: needs_ext_seg_data  !< Returns True if external segment data needed
  open_boundary_query = .false.
  if (.not. associated(obc)) return
  if (present(apply_open_obc)) open_boundary_query = obc%open_u_BCs_exist_globally .or. &
                                                     obc%open_v_BCs_exist_globally
  if (present(apply_specified_obc)) open_boundary_query = obc%specified_u_BCs_exist_globally .or. &
                                                          obc%specified_v_BCs_exist_globally
  if (present(apply_flather_obc)) open_boundary_query = obc%Flather_u_BCs_exist_globally .or. &
                                                        obc%Flather_v_BCs_exist_globally
  if (present(apply_nudged_obc)) open_boundary_query = obc%nudged_u_BCs_exist_globally .or. &
                                                       obc%nudged_v_BCs_exist_globally
  if (present(needs_ext_seg_data)) open_boundary_query = obc%needs_IO_for_data
 
end function open_boundary_query
 
!> Deallocate open boundary data
subroutine open_boundary_dealloc(OBC)
  type(ocean_obc_type), pointer :: OBC !< Open boundary control structure
  type(obc_segment_type), pointer :: segment => null()
  integer :: n
 
  if (.not. associated(obc)) return
 
  do n=1, obc%number_of_segments
    segment => obc%segment(n)
    call deallocate_obc_segment_data(obc, segment)
  enddo
  if (associated(obc%segment)) deallocate(obc%segment)
  if (associated(obc%segnum_u)) deallocate(obc%segnum_u)
  if (associated(obc%segnum_v)) deallocate(obc%segnum_v)
  if (associated(obc%rx_normal)) deallocate(obc%rx_normal)
  if (associated(obc%ry_normal)) deallocate(obc%ry_normal)
  if (associated(obc%rx_oblique)) deallocate(obc%rx_oblique)
  if (associated(obc%ry_oblique)) deallocate(obc%ry_oblique)
  if (associated(obc%cff_normal)) deallocate(obc%cff_normal)
  if (associated(obc%tres_x)) deallocate(obc%tres_x)
  if (associated(obc%tres_y)) deallocate(obc%tres_y)
  deallocate(obc)
end subroutine open_boundary_dealloc
 
!> Close open boundary data
subroutine open_boundary_end(OBC)
  type(ocean_obc_type), pointer :: obc !< Open boundary control structure
  call open_boundary_dealloc(obc)
end subroutine open_boundary_end
 
!> Sets the slope of bathymetry normal to an open bounndary to zero.
subroutine open_boundary_impose_normal_slope(OBC, G, depth)
  type(ocean_obc_type),             pointer       :: obc !< Open boundary control structure
  type(dyn_horgrid_type),           intent(in)    :: g !< Ocean grid structure
  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points
  ! Local variables
  integer :: i, j, n
  type(obc_segment_type), pointer :: segment => null()
 
  if (.not.associated(obc)) return
 
  if (.not.(obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
              obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
    return
 
  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
        depth(i+1,j) = depth(i,j)
      enddo
    elseif (segment%direction == obc_direction_w) then
      i=segment%HI%IsdB
      do j=segment%HI%jsd,segment%HI%jed
        depth(i,j) = depth(i+1,j)
      enddo
    elseif (segment%direction == obc_direction_n) then
      j=segment%HI%JsdB
      do i=segment%HI%isd,segment%HI%ied
        depth(i,j+1) = depth(i,j)
      enddo
    elseif (segment%direction == obc_direction_s) then
      j=segment%HI%JsdB
      do i=segment%HI%isd,segment%HI%ied
        depth(i,j) = depth(i,j+1)
      enddo
    endif
  enddo
 
end subroutine open_boundary_impose_normal_slope
 
!> Reconcile masks and open boundaries, deallocate OBC on PEs where it is not needed.
!! Also adjust u- and v-point cell area on specified open boundaries and mask all
!! points outside open boundaries.
subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US)
  type(ocean_obc_type),              pointer       :: obc !< Open boundary control structure
  type(dyn_horgrid_type),            intent(inout) :: g   !< Ocean grid structure
  type(unit_scale_type),             intent(in)    :: us  !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: areacu !< Area of a u-cell [L2 ~> m2]
  real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: areacv !< Area of a u-cell [L2 ~> m2]
  ! Local variables
  integer :: i, j, n
  type(obc_segment_type), pointer :: segment => null()
  logical :: any_u, any_v
 
  if (.not.associated(obc)) return
 
  do n=1,obc%number_of_segments
    segment=>obc%segment(n)
    if (.not. segment%on_pe) cycle
    if (segment%is_E_or_W) then
      ! Sweep along u-segments and delete the OBC for blocked points.
      ! Also, mask all points outside.
      i=segment%HI%IsdB
      do j=segment%HI%jsd,segment%HI%jed
        if (g%mask2dCu(i,j) == 0) obc%segnum_u(i,j) = obc_none
        if (segment%direction == obc_direction_w) then
          g%mask2dT(i,j) = 0
        else
          g%mask2dT(i+1,j) = 0
        endif
      enddo
      do j=segment%HI%JsdB+1,segment%HI%JedB-1
        if (segment%direction == obc_direction_w) then
          g%mask2dCv(i,j) = 0
        else
          g%mask2dCv(i+1,j) = 0
        endif
      enddo
    else
      ! Sweep along v-segments and delete the OBC for blocked points.
      j=segment%HI%JsdB
      do i=segment%HI%isd,segment%HI%ied
        if (g%mask2dCv(i,j) == 0) obc%segnum_v(i,j) = obc_none
        if (segment%direction == obc_direction_s) then
          g%mask2dT(i,j) = 0
        else
          g%mask2dT(i,j+1) = 0
        endif
      enddo
      do i=segment%HI%IsdB+1,segment%HI%IedB-1
        if (segment%direction == obc_direction_s) then
          g%mask2dCu(i,j) = 0
        else
          g%mask2dCu(i,j+1) = 0
        endif
      enddo
    endif
  enddo
 
  do n=1,obc%number_of_segments
    segment=>obc%segment(n)
    if (.not. segment%on_pe .or. .not. segment%specified) cycle
    if (segment%is_E_or_W) then
      ! Sweep along u-segments and for %specified BC points reset the u-point area which was masked out
      i=segment%HI%IsdB
      do j=segment%HI%jsd,segment%HI%jed
        if (segment%direction == obc_direction_e) then
          areacu(i,j) = g%areaT(i,j)   ! Both of these are in [L2]
        else   ! West
          areacu(i,j) = g%areaT(i+1,j) ! Both of these are in [L2]
        endif
      enddo
    else
      ! Sweep along v-segments and for %specified BC points reset the v-point area which was masked out
      j=segment%HI%JsdB
      do i=segment%HI%isd,segment%HI%ied
        if (segment%direction == obc_direction_s) then
          areacv(i,j) = g%areaT(i,j+1) ! Both of these are in [L2]
        else      ! North
          areacu(i,j) = g%areaT(i,j)   ! Both of these are in [L2]
        endif
      enddo
    endif
  enddo
 
  ! G%mask2du will be open wherever bathymetry allows it.
  ! Bathymetry outside of the open boundary was adjusted to match
  ! the bathymetry inside so these points will be open unless the
  ! bathymetry inside the boundary was too shallow and flagged as land.
  any_u = .false.
  any_v = .false.
  do n=1,obc%number_of_segments
    segment=>obc%segment(n)
    if (.not. segment%on_pe) cycle
    if (segment%is_E_or_W) then
      i=segment%HI%IsdB
      do j=segment%HI%jsd,segment%HI%jed
        if (obc%segnum_u(i,j) /= obc_none) any_u = .true.
      enddo
    else
      j=segment%HI%JsdB
      do i=segment%HI%isd,segment%HI%ied
        if (obc%segnum_v(i,j) /= obc_none) any_v = .true.
      enddo
    endif
  enddo
 
  obc%OBC_pe = .true.
  if (.not.(any_u .or. any_v)) obc%OBC_pe = .false.
 
end subroutine open_boundary_impose_land_mask
 
!> Make sure the OBC tracer reservoirs are initialized.
subroutine setup_obc_tracer_reservoirs(G, OBC)
  type(ocean_grid_type),      intent(inout) :: G          !< Ocean grid structure
  type(ocean_obc_type),       pointer       :: OBC !< Open boundary control structure
  ! Local variables
  type(obc_segment_type), pointer :: segment => null()
  integer :: i, j, k, m, n
 
  do n=1,obc%number_of_segments
    segment=>obc%segment(n)
    if (associated(segment%tr_Reg)) then
      if (segment%is_E_or_W) then
        i = segment%HI%IsdB
        do m=1,obc%ntr
          if (associated(segment%tr_Reg%Tr(m)%tres)) then
            do k=1,g%ke
              do j=segment%HI%jsd,segment%HI%jed
                obc%tres_x(i,j,k,m) = segment%tr_Reg%Tr(m)%t(i,j,k)
              enddo
            enddo
          endif
        enddo
      else
        j = segment%HI%JsdB
        do m=1,obc%ntr
          if (associated(segment%tr_Reg%Tr(m)%tres)) then
            do k=1,g%ke
              do i=segment%HI%isd,segment%HI%ied
                obc%tres_y(i,j,k,m) = segment%tr_Reg%Tr(m)%t(i,j,k)
              enddo
            enddo
          endif
        enddo
      endif
    endif
  enddo
 
end subroutine setup_obc_tracer_reservoirs
 
!> Apply radiation conditions to 3D  u,v at open boundaries
subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, US, dt)
  type(ocean_grid_type),                     intent(inout) :: g !< Ocean grid structure
  type(ocean_obc_type),                      pointer       :: obc !< Open boundary control structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new !< On exit, new u values on open boundaries
                                                                    !! On entry, the old time-level v but including
                                                                    !! barotropic accelerations [L T-1 ~> m s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u_old !< Original unadjusted u [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new !< On exit, new v values on open boundaries.
                                                                    !! On entry, the old time-level v but including
                                                                    !! barotropic accelerations [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v_old !< Original unadjusted v  [L T-1 ~> m s-1]
  type(unit_scale_type),                     intent(in)    :: us    !< A dimensional unit scaling type
  real,                                      intent(in)    :: dt    !< Appropriate timestep [T ~> s]
  ! Local variables
  real :: dhdt, dhdx, dhdy  ! One-point differences in time or space [L T-1 ~> m s-1]
  real :: gamma_u, gamma_2  ! Fractional weightings of new values [nondim]
  real :: tau            ! A local nudging timescale [T ~> s]
  real :: rx_max, ry_max ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2]
  real :: rx_new, rx_avg ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2]
  real :: ry_new, ry_avg ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2]
  real :: cff_new, cff_avg ! denominator in oblique [L2 T-2 ~> m2 s-2]
  real, allocatable, dimension(:,:,:) :: &
    rx_tang_rad, & ! The phase speed at u-points for tangential oblique OBCs
                   ! in units of grid points per timestep [nondim]
    ry_tang_rad, & ! The phase speed at v-points for tangential oblique OBCs
                   ! in units of grid points per timestep [nondim]
    rx_tang_obl, & ! The x-coefficient for tangential oblique OBCs [L2 T-2 ~> m2 s-2]
    ry_tang_obl, & ! The y-coefficient for tangential oblique OBCs [L2 T-2 ~> m2 s-2]
    cff_tangential ! The denominator for tangential oblique OBCs [L2 T-2 ~> m2 s-2]
  real :: eps      ! A small velocity squared [L2 T-2 ~> m2 s-2]
  type(obc_segment_type), pointer :: segment => null()
  integer :: i, j, k, is, ie, js, je, m, nz, n
  integer :: is_obc, ie_obc, js_obc, je_obc
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
 
  if (.not.associated(obc)) return
 
  if (.not.(obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
    return
 
  eps = 1.0e-20*us%m_s_to_L_T**2
 
  !! Copy previously calculated phase velocity from global arrays into segments
  !! This is terribly inefficient and temporary solution for continuity across restarts
  !! and needs to be revisited in the future.
  if (obc%gamma_uv < 1.0) then
    do n=1,obc%number_of_segments
      segment=>obc%segment(n)
      if (.not. segment%on_pe) cycle
      if (segment%is_E_or_W .and. segment%radiation) then
        do k=1,g%ke
          i=segment%HI%IsdB
          do j=segment%HI%jsd,segment%HI%jed
            segment%rx_norm_rad(i,j,k) = obc%rx_normal(i,j,k)
          enddo
        enddo
      elseif (segment%is_N_or_S .and. segment%radiation) then
        do k=1,g%ke
          j=segment%HI%JsdB
          do i=segment%HI%isd,segment%HI%ied
            segment%ry_norm_rad(i,j,k) = obc%ry_normal(i,j,k)
          enddo
        enddo
      endif
      if (segment%is_E_or_W .and. segment%oblique) then
        do k=1,g%ke
          i=segment%HI%IsdB
          do j=segment%HI%jsd,segment%HI%jed
            segment%rx_norm_obl(i,j,k) = obc%rx_oblique(i,j,k)
            segment%ry_norm_obl(i,j,k) = obc%ry_oblique(i,j,k)
            segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
          enddo
        enddo
      elseif (segment%is_N_or_S .and. segment%oblique) then
        do k=1,g%ke
          j=segment%HI%JsdB
          do i=segment%HI%isd,segment%HI%ied
            segment%rx_norm_obl(i,j,k) = obc%rx_oblique(i,j,k)
            segment%ry_norm_obl(i,j,k) = obc%ry_oblique(i,j,k)
            segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
          enddo
        enddo
      endif
    enddo
  endif
 
  ! Now tracers (if any)
  do n=1,obc%number_of_segments
    segment=>obc%segment(n)
    if (associated(segment%tr_Reg)) then
      if (segment%is_E_or_W) then
        i = segment%HI%IsdB
        do m=1,obc%ntr
          if (associated(segment%tr_Reg%Tr(m)%tres)) then
            do k=1,g%ke
              do j=segment%HI%jsd,segment%HI%jed
                segment%tr_Reg%Tr(m)%tres(i,j,k) = obc%tres_x(i,j,k,m)
              enddo
            enddo
          endif
        enddo
      else
        j = segment%HI%JsdB
        do m=1,obc%ntr
          if (associated(segment%tr_Reg%Tr(m)%tres)) then
            do k=1,g%ke
              do i=segment%HI%isd,segment%HI%ied
                segment%tr_Reg%Tr(m)%tres(i,j,k) = obc%tres_y(i,j,k,m)
              enddo
            enddo
          endif
        enddo
      endif
    endif
  enddo
 
  gamma_u = obc%gamma_uv
  rx_max = obc%rx_max ; ry_max = obc%rx_max
  do n=1,obc%number_of_segments
     segment=>obc%segment(n)
     if (.not. segment%on_pe) cycle
     if (segment%oblique) call gradient_at_q_points(g, segment, u_new(:,:,:), v_new(:,:,:))
     if (segment%direction == obc_direction_e) then
       i=segment%HI%IsdB
       if (i<g%HI%IscB) cycle
       do k=1,nz ;  do j=segment%HI%jsd,segment%HI%jed
         if (segment%radiation) then
           dhdt = (u_old(i-1,j,k) - u_new(i-1,j,k)) !old-new
           dhdx = (u_new(i-1,j,k) - u_new(i-2,j,k)) !in new time backward sasha for I-1
           rx_new = 0.0
           if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max) ! outward phase speed
           if (gamma_u < 1.0) then
             rx_avg = (1.0-gamma_u)*segment%rx_norm_rad(i,j,k) + gamma_u*rx_new
           else
             rx_avg = rx_new
           endif
           segment%rx_norm_rad(i,j,k) = rx_avg
           ! The new boundary value is interpolated between future interior
           ! value, u_new(I-1) and past boundary value but with barotropic
           ! accelerations, u_new(I).
           segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) / (1.0+rx_avg)
           ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
           ! implemented as a work-around to limitations in restart capability
           if (gamma_u < 1.0) then
             obc%rx_normal(i,j,k) = segment%rx_norm_rad(i,j,k)
           endif
         elseif (segment%oblique) then
           dhdt = (u_old(i-1,j,k) - u_new(i-1,j,k)) !old-new
           dhdx = (u_new(i-1,j,k) - u_new(i-2,j,k)) !in new time backward sasha for I-1
           if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0) then
             dhdy = segment%grad_normal(j-1,1,k)
           elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0) then
             dhdy = 0.0
           else
             dhdy = segment%grad_normal(j,1,k)
           endif
           if (dhdt*dhdx < 0.0) dhdt = 0.0
           rx_new = dhdt*dhdx
           cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
           ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
           if (gamma_u < 1.0) then
             rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
             ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
             cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
           else
             rx_avg = rx_new
             ry_avg = ry_new
             cff_avg = cff_new
           endif
           segment%rx_norm_obl(i,j,k) = rx_avg
           segment%ry_norm_obl(i,j,k) = ry_avg
           segment%cff_normal(i,j,k) = cff_avg
           segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) - &
                              (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + &
                               min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
                            (cff_avg + rx_avg)
           if (gamma_u < 1.0) then
             ! Copy restart fields into 3-d arrays. This is an inefficient and temporary
             ! implementation as a work-around to limitations in restart capability
             obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
             obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
             obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
           endif
         elseif (segment%gradient) then
           segment%normal_vel(i,j,k) = u_new(i-1,j,k)
         endif
         if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then
           ! dhdt gets set to 0 on inflow in oblique case
           if (dhdt*dhdx <= 0.0) then
             tau = segment%Velocity_nudging_timescale_in
           else
             tau = segment%Velocity_nudging_timescale_out
           endif
           gamma_2 = dt / (tau + dt)
           segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_normal_vel(i,j,k)
         endif
       enddo ; enddo
       if (segment%radiation_tan .or. segment%radiation_grad) then
         i=segment%HI%IsdB
         allocate(rx_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         do k=1,nz
           if (gamma_u < 1.0) then
             rx_tang_rad(i,segment%HI%JsdB,k) = segment%rx_norm_rad(i,segment%HI%jsd,k)
             rx_tang_rad(i,segment%HI%JedB,k) = segment%rx_norm_rad(i,segment%HI%jed,k)
             do j=segment%HI%JsdB+1,segment%HI%JedB-1
               rx_tang_rad(i,j,k) = 0.5*(segment%rx_norm_rad(i,j,k) + segment%rx_norm_rad(i,j+1,k))
             enddo
           else
             do j=segment%HI%JsdB,segment%HI%JedB
               dhdt = v_old(i,j,k)-v_new(i,j,k)   !old-new
               dhdx = v_new(i,j,k)-v_new(i-1,j,k) !in new time backward sasha for I-1
               rx_tang_rad(i,j,k) = 0.0
               if (dhdt*dhdx > 0.0) rx_tang_rad(i,j,k) = min( (dhdt/dhdx), rx_max) ! outward phase speed
             enddo
           endif
         enddo
         if (segment%radiation_tan) then
           do k=1,nz ;  do j=segment%HI%JsdB,segment%HI%JedB
             rx_avg = rx_tang_rad(i,j,k)
             segment%tangential_vel(i,j,k) = (v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) / (1.0+rx_avg)
           enddo ; enddo
         endif
         if (segment%nudged_tan) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (rx_tang_rad(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_vel(i,j,k)
           enddo ; enddo
         endif
         if (segment%radiation_grad) then
           js_obc = max(segment%HI%JsdB,g%jsd+1)
           je_obc = min(segment%HI%JedB,g%jed-1)
           do k=1,nz ; do j=js_obc,je_obc
             rx_avg = rx_tang_rad(i,j,k)
!            if (G%mask2dCu(I-1,j) > 0.0 .and. G%mask2dCu(I-1,j+1) > 0.0) then
!              rx_avg = 0.5*(u_new(I-1,j,k) + u_new(I-1,j+1,k)) * dt * G%IdxBu(I-1,J)
!            elseif (G%mask2dCu(I-1,j) > 0.0) then
!              rx_avg = u_new(I-1,j,k) * dt * G%IdxBu(I-1,J)
!            elseif (G%mask2dCu(I-1,j+1) > 0.0) then
!              rx_avg = u_new(I-1,j+1,k) * dt * G%IdxBu(I-1,J)
!            else
!              rx_avg = 0.0
!            endif
             segment%tangential_grad(i,j,k) = ((v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
                               rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) / (1.0+rx_avg)
           enddo ; enddo
         endif
         if (segment%nudged_grad) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (rx_tang_rad(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_grad(i,j,k)
           enddo ; enddo
         endif
         deallocate(rx_tang_rad)
       endif
       if (segment%oblique_tan .or. segment%oblique_grad) then
         i=segment%HI%IsdB
         allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         do k=1,nz
           if (gamma_u < 1.0) then
             rx_tang_obl(i,segment%HI%JsdB,k) = segment%rx_norm_obl(i,segment%HI%jsd,k)
             rx_tang_obl(i,segment%HI%JedB,k) = segment%rx_norm_obl(i,segment%HI%jed,k)
             ry_tang_obl(i,segment%HI%JsdB,k) = segment%ry_norm_obl(i,segment%HI%jsd,k)
             ry_tang_obl(i,segment%HI%JedB,k) = segment%ry_norm_obl(i,segment%HI%jed,k)
             cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
             cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
             do j=segment%HI%JsdB+1,segment%HI%JedB-1
               rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i,j+1,k))
               ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i,j+1,k))
               cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
             enddo
           else
             do j=segment%HI%JsdB,segment%HI%JedB
               dhdt = v_old(i,j,k)-v_new(i,j,k)   !old-new
               dhdx = v_new(i,j,k)-v_new(i-1,j,k) !in new time backward sasha for I-1
               if (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) > 0.0) then
                 dhdy = segment%grad_tan(j,1,k)
               elseif (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) == 0.0) then
                 dhdy = 0.0
               else
                 dhdy = segment%grad_tan(j+1,1,k)
               endif
               if (dhdt*dhdx < 0.0) dhdt = 0.0
               rx_new = dhdt*dhdx
               cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
               ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
               rx_tang_obl(i,j,k) = rx_new
               ry_tang_obl(i,j,k) = ry_new
               cff_tangential(i,j,k) = cff_new
             enddo
           endif
         enddo
         if (segment%oblique_tan) then
           do k=1,nz ;  do j=segment%HI%JsdB,segment%HI%JedB
             rx_avg = rx_tang_obl(i,j,k)
             ry_avg = ry_tang_obl(i,j,k)
             cff_avg = cff_tangential(i,j,k)
             segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) - &
                                              (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + &
                                               min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
                                             (cff_avg + rx_avg)
           enddo ; enddo
         endif
         if (segment%nudged_tan) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (rx_tang_obl(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_vel(i,j,k)
           enddo ; enddo
         endif
         if (segment%oblique_grad) then
           js_obc = max(segment%HI%JsdB,g%jsd+1)
           je_obc = min(segment%HI%JedB,g%jed-1)
           do k=1,nz ;  do j=segment%HI%JsdB+1,segment%HI%JedB-1
             rx_avg = rx_tang_obl(i,j,k)
             ry_avg = ry_tang_obl(i,j,k)
             cff_avg = cff_tangential(i,j,k)
             segment%tangential_grad(i,j,k) =  &
                 ((cff_avg*(v_new(i,j,k)  - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
                   rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) - &
                  (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + &
                   min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k)) ) / &
                 (cff_avg + rx_avg)
           enddo ; enddo
         endif
         if (segment%nudged_grad) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (rx_tang_obl(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_grad(i,j,k)
           enddo ; enddo
         endif
         deallocate(rx_tang_obl)
         deallocate(ry_tang_obl)
         deallocate(cff_tangential)
       endif
     endif
 
     if (segment%direction == obc_direction_w) then
       i=segment%HI%IsdB
       if (i>g%HI%IecB) cycle
       do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
         if (segment%radiation) then
           dhdt = (u_old(i+1,j,k) - u_new(i+1,j,k)) !old-new
           dhdx = (u_new(i+1,j,k) - u_new(i+2,j,k)) !in new time forward sasha for I+1
           rx_new = 0.0
           if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
           if (gamma_u < 1.0) then
             rx_avg = (1.0-gamma_u)*segment%rx_norm_rad(i,j,k) + gamma_u*rx_new
           else
             rx_avg = rx_new
           endif
           segment%rx_norm_rad(i,j,k) = rx_avg
           ! The new boundary value is interpolated between future interior
           ! value, u_new(I+1) and past boundary value but with barotropic
           ! accelerations, u_new(I).
           segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) / (1.0+rx_avg)
           if (gamma_u < 1.0) then
             ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
             ! implemented as a work-around to limitations in restart capability
             obc%rx_normal(i,j,k) = segment%rx_norm_rad(i,j,k)
           endif
         elseif (segment%oblique) then
           dhdt = (u_old(i+1,j,k) - u_new(i+1,j,k)) !old-new
           dhdx = (u_new(i+1,j,k) - u_new(i+2,j,k)) !in new time forward sasha for I+1
           if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0) then
             dhdy = segment%grad_normal(j-1,1,k)
           elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0) then
             dhdy = 0.0
           else
             dhdy = segment%grad_normal(j,1,k)
           endif
           if (dhdt*dhdx < 0.0) dhdt = 0.0
 
           rx_new = dhdt*dhdx
           cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
           ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
           if (gamma_u < 1.0) then
             rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
             ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
             cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
           else
             rx_avg = rx_new
             ry_avg = ry_new
             cff_avg = cff_new
           endif
           segment%rx_norm_obl(i,j,k) = rx_avg
           segment%ry_norm_obl(i,j,k) = ry_avg
           segment%cff_normal(i,j,k) = cff_avg
           segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) - &
                                        (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + &
                                         min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
                                       (cff_avg + rx_avg)
           if (gamma_u < 1.0) then
             ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
             ! implemented as a work-around to limitations in restart capability
             obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
             obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
             obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
           endif
         elseif (segment%gradient) then
           segment%normal_vel(i,j,k) = u_new(i+1,j,k)
         endif
         if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then
           ! dhdt gets set to 0. on inflow in oblique case
           if (dhdt*dhdx <= 0.0) then
             tau = segment%Velocity_nudging_timescale_in
           else
             tau = segment%Velocity_nudging_timescale_out
           endif
           gamma_2 = dt / (tau + dt)
           segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_normal_vel(i,j,k)
         endif
       enddo ; enddo
       if (segment%radiation_tan .or. segment%radiation_grad) then
         i=segment%HI%IsdB
         allocate(rx_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         do k=1,nz
           if (gamma_u < 1.0) then
             rx_tang_rad(i,segment%HI%JsdB,k) = segment%rx_norm_rad(i,segment%HI%jsd,k)
             rx_tang_rad(i,segment%HI%JedB,k) = segment%rx_norm_rad(i,segment%HI%jed,k)
             do j=segment%HI%JsdB+1,segment%HI%JedB-1
               rx_tang_rad(i,j,k) = 0.5*(segment%rx_norm_rad(i,j,k) + segment%rx_norm_rad(i,j+1,k))
             enddo
           else
             do j=segment%HI%JsdB,segment%HI%JedB
               dhdt = v_old(i+1,j,k)-v_new(i+1,j,k)   !old-new
               dhdx = v_new(i+1,j,k)-v_new(i+2,j,k) !in new time backward sasha for I-1
               rx_tang_rad(i,j,k) = 0.0
               if (dhdt*dhdx > 0.0) rx_tang_rad(i,j,k) = min( (dhdt/dhdx), rx_max) ! outward phase speed
             enddo
           endif
         enddo
         if (segment%radiation_tan) then
           do k=1,nz ;  do j=segment%HI%JsdB,segment%HI%JedB
             rx_avg = rx_tang_rad(i,j,k)
             segment%tangential_vel(i,j,k) = (v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) / (1.0+rx_avg)
           enddo ; enddo
         endif
         if (segment%nudged_tan) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (rx_tang_rad(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_vel(i,j,k)
           enddo ; enddo
         endif
         if (segment%radiation_grad) then
           js_obc = max(segment%HI%JsdB,g%jsd+1)
           je_obc = min(segment%HI%JedB,g%jed-1)
           do k=1,nz ;  do j=js_obc,je_obc
             rx_avg = rx_tang_rad(i,j,k)
!            if (G%mask2dCu(I+1,j) > 0.0 .and. G%mask2dCu(I+1,j+1) > 0.0) then
!              rx_avg = 0.5*(u_new(I+1,j,k) + u_new(I+1,j+1,k)) * dt * G%IdxBu(I+1,J)
!            elseif (G%mask2dCu(I+1,j) > 0.0) then
!              rx_avg = u_new(I+1,j,k) * dt * G%IdxBu(I+1,J)
!            elseif (G%mask2dCu(I+1,j+1) > 0.0) then
!              rx_avg = u_new(I+1,j+1,k) * dt * G%IdxBu(I+1,J)
!            else
!              rx_avg = 0.0
!            endif
             segment%tangential_grad(i,j,k) = ((v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
                               rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) / (1.0+rx_avg)
           enddo ; enddo
         endif
         if (segment%nudged_grad) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (rx_tang_rad(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_grad(i,j,k)
           enddo ; enddo
         endif
         deallocate(rx_tang_rad)
       endif
       if (segment%oblique_tan .or. segment%oblique_grad) then
         i=segment%HI%IsdB
         allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         do k=1,nz
           if (gamma_u < 1.0) then
             rx_tang_obl(i,segment%HI%JsdB,k) = segment%rx_norm_obl(i,segment%HI%jsd,k)
             rx_tang_obl(i,segment%HI%JedB,k) = segment%rx_norm_obl(i,segment%HI%jed,k)
             ry_tang_obl(i,segment%HI%JsdB,k) = segment%ry_norm_obl(i,segment%HI%jsd,k)
             ry_tang_obl(i,segment%HI%JedB,k) = segment%ry_norm_obl(i,segment%HI%jed,k)
             cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
             cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
             do j=segment%HI%JsdB+1,segment%HI%JedB-1
               rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i,j+1,k))
               ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i,j+1,k))
               cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
             enddo
           else
             do j=segment%HI%JsdB,segment%HI%JedB
               dhdt = v_old(i+1,j,k)-v_new(i+1,j,k)   !old-new
               dhdx = v_new(i+1,j,k)-v_new(i+2,j,k) !in new time backward sasha for I-1
               if (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) > 0.0) then
                 dhdy = segment%grad_tan(j,1,k)
               elseif (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) == 0.0) then
                 dhdy = 0.0
               else
                 dhdy = segment%grad_tan(j+1,1,k)
               endif
               if (dhdt*dhdx < 0.0) dhdt = 0.0
               rx_new = dhdt*dhdx
               cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
               ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
               rx_tang_obl(i,j,k) = rx_new
               ry_tang_obl(i,j,k) = ry_new
               cff_tangential(i,j,k) = cff_new
             enddo
           endif
         enddo
         if (segment%oblique_tan) then
           do k=1,nz ;  do j=segment%HI%JsdB,segment%HI%JedB
             rx_avg = rx_tang_obl(i,j,k)
             ry_avg = ry_tang_obl(i,j,k)
             cff_avg = cff_tangential(i,j,k)
             segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) - &
                                              (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + &
                                               min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
                                             (cff_avg + rx_avg)
           enddo ; enddo
         endif
         if (segment%nudged_tan) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (rx_tang_obl(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_vel(i,j,k)
           enddo ; enddo
         endif
         if (segment%oblique_grad) then
           js_obc = max(segment%HI%JsdB,g%jsd+1)
           je_obc = min(segment%HI%JedB,g%jed-1)
           do k=1,nz ;  do j=segment%HI%JsdB+1,segment%HI%JedB-1
             rx_avg = rx_tang_obl(i,j,k)
             ry_avg = ry_tang_obl(i,j,k)
             cff_avg = cff_tangential(i,j,k)
             segment%tangential_grad(i,j,k) = &
                 ((cff_avg*(v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
                    rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) - &
                  (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + &
                   min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k))) / &
                 (cff_avg + rx_avg)
           enddo ; enddo
         endif
         if (segment%nudged_grad) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (rx_tang_obl(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_grad(i,j,k)
           enddo ; enddo
         endif
         deallocate(rx_tang_obl)
         deallocate(ry_tang_obl)
         deallocate(cff_tangential)
       endif
     endif
 
     if (segment%direction == obc_direction_n) then
       j=segment%HI%JsdB
       if (j<g%HI%JscB) cycle
       do k=1,nz ;  do i=segment%HI%isd,segment%HI%ied
         if (segment%radiation) then
           dhdt = (v_old(i,j-1,k) - v_new(i,j-1,k)) !old-new
           dhdy = (v_new(i,j-1,k) - v_new(i,j-2,k)) !in new time backward sasha for J-1
           ry_new = 0.0
           if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
           if (gamma_u < 1.0) then
             ry_avg = (1.0-gamma_u)*segment%ry_norm_rad(i,j,k) + gamma_u*ry_new
           else
             ry_avg = ry_new
           endif
           segment%ry_norm_rad(i,j,k) = ry_avg
           ! The new boundary value is interpolated between future interior
           ! value, v_new(J-1) and past boundary value but with barotropic
           ! accelerations, v_new(J).
           segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) / (1.0+ry_avg)
           if (gamma_u < 1.0) then
             ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
             ! implemented as a work-around to limitations in restart capability
             obc%ry_normal(i,j,k) = segment%ry_norm_rad(i,j,k)
           endif
         elseif (segment%oblique) then
           dhdt = (v_old(i,j-1,k) - v_new(i,j-1,k)) !old-new
           dhdy = (v_new(i,j-1,k) - v_new(i,j-2,k)) !in new time backward sasha for J-1
           if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0) then
             dhdx = segment%grad_normal(i-1,1,k)
           elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0) then
             dhdx = 0.0
           else
             dhdx = segment%grad_normal(i,1,k)
           endif
           if (dhdt*dhdy < 0.0) dhdt = 0.0
           ry_new = dhdt*dhdy
           cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
           rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
           if (gamma_u < 1.0) then
             rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
             ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
             cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
           else
             rx_avg = rx_new
             ry_avg = ry_new
             cff_avg = cff_new
           endif
           segment%rx_norm_obl(i,j,k) = rx_avg
           segment%ry_norm_obl(i,j,k) = ry_avg
           segment%cff_normal(i,j,k) = cff_avg
           segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) - &
                                        (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) +&
                                         min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
                                       (cff_avg + ry_avg)
           if (gamma_u < 1.0) then
             ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
             ! implemented as a work-around to limitations in restart capability
             obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
             obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
             obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
           endif
         elseif (segment%gradient) then
           segment%normal_vel(i,j,k) = v_new(i,j-1,k)
         endif
         if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then
           ! dhdt gets set to 0 on inflow in oblique case
           if (dhdt*dhdy <= 0.0) then
             tau = segment%Velocity_nudging_timescale_in
           else
             tau = segment%Velocity_nudging_timescale_out
           endif
           gamma_2 = dt / (tau + dt)
           segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_normal_vel(i,j,k)
         endif
       enddo ; enddo
       if (segment%radiation_tan .or. segment%radiation_grad) then
         j=segment%HI%JsdB
         allocate(ry_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         do k=1,nz
           if (gamma_u < 1.0) then
             ry_tang_rad(segment%HI%IsdB,j,k) = segment%ry_norm_rad(segment%HI%isd,j,k)
             ry_tang_rad(segment%HI%IedB,j,k) = segment%ry_norm_rad(segment%HI%ied,j,k)
             do i=segment%HI%IsdB+1,segment%HI%IedB-1
               ry_tang_rad(i,j,k) = 0.5*(segment%ry_norm_rad(i,j,k) + segment%ry_norm_rad(i+1,j,k))
             enddo
           else
             do i=segment%HI%IsdB,segment%HI%IedB
               dhdt = u_old(i,j-1,k)-u_new(i,j-1,k)   !old-new
               dhdy = u_new(i,j-1,k)-u_new(i,j-2,k) !in new time backward sasha for I-1
               ry_tang_rad(i,j,k) = 0.0
               if (dhdt*dhdy > 0.0) ry_tang_rad(i,j,k) = min( (dhdt/dhdy), rx_max) ! outward phase speed
             enddo
           endif
         enddo
         if (segment%radiation_tan) then
           do k=1,nz ;  do i=segment%HI%IsdB,segment%HI%IedB
             ry_avg = ry_tang_rad(i,j,k)
             segment%tangential_vel(i,j,k) = (u_new(i,j,k) + ry_avg*u_new(i,j-1,k)) / (1.0+ry_avg)
           enddo ; enddo
         endif
         if (segment%nudged_tan) then
           do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (ry_tang_rad(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_vel(i,j,k)
           enddo ; enddo
         endif
         if (segment%radiation_grad) then
           is_obc = max(segment%HI%IsdB,g%isd+1)
           ie_obc = min(segment%HI%IedB,g%ied-1)
           do k=1,nz ;  do i=is_obc,ie_obc
             ry_avg = ry_tang_rad(i,j,k)
!            if (G%mask2dCv(i,J-1) > 0.0 .and. G%mask2dCv(i+1,J-1) > 0.0) then
!              ry_avg = 0.5*(v_new(i,J-1,k) + v_new(i+1,J-1,k) * dt * G%IdyBu(I,J-1))
!            elseif (G%mask2dCv(i,J-1) > 0.0) then
!              ry_avg = v_new(i,J-1,k) * dt *G%IdyBu(I,J-1)
!            elseif (G%mask2dCv(i+1,J-1) > 0.0) then
!              ry_avg = v_new(i+1,J-1,k) * dt *G%IdyBu(I,J-1)
!            else
!              ry_avg = 0.0
!            endif
             segment%tangential_grad(i,j,k) = ((u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
                               ry_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) / (1.0+ry_avg)
           enddo ; enddo
         endif
         if (segment%nudged_grad) then
           do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (ry_tang_rad(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_grad(i,j,k)
           enddo ; enddo
         endif
         deallocate(ry_tang_rad)
       endif
       if (segment%oblique_tan .or. segment%oblique_grad) then
         j=segment%HI%JsdB
         allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         do k=1,nz
           if (gamma_u < 1.0) then
             rx_tang_obl(segment%HI%IsdB,j,k) = segment%rx_norm_obl(segment%HI%isd,j,k)
             rx_tang_obl(segment%HI%IedB,j,k) = segment%rx_norm_obl(segment%HI%ied,j,k)
             ry_tang_obl(segment%HI%IsdB,j,k) = segment%ry_norm_obl(segment%HI%isd,j,k)
             ry_tang_obl(segment%HI%IedB,j,k) = segment%ry_norm_obl(segment%HI%ied,j,k)
             cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
             cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
             do i=segment%HI%IsdB+1,segment%HI%IedB-1
               rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i+1,j,k))
               ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i+1,j,k))
               cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
             enddo
           else
             do i=segment%HI%IsdB,segment%HI%IedB
               dhdt = u_old(i,j,k)-u_new(i,j,k)   !old-new
               dhdy = u_new(i,j,k)-u_new(i,j-1,k) !in new time backward sasha for I-1
               if (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) > 0.0) then
                 dhdx = segment%grad_tan(i,1,k)
               elseif (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) == 0.0) then
                 dhdx = 0.0
               else
                 dhdx = segment%grad_tan(i+1,1,k)
               endif
               if (dhdt*dhdy < 0.0) dhdt = 0.0
               ry_new = dhdt*dhdy
               cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
               rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
               rx_tang_obl(i,j,k) = rx_new
               ry_tang_obl(i,j,k) = ry_new
               cff_tangential(i,j,k) = cff_new
             enddo
           endif
         enddo
         if (segment%oblique_tan) then
           do k=1,nz ;  do i=segment%HI%IsdB,segment%HI%IedB
             rx_avg = rx_tang_obl(i,j,k)
             ry_avg = ry_tang_obl(i,j,k)
             cff_avg = cff_tangential(i,j,k)
             segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + ry_avg*u_new(i,j-1,k)) - &
                                              (max(rx_avg,0.0)*segment%grad_tan(i,2,k) + &
                                               min(rx_avg,0.0)*segment%grad_tan(i+1,2,k))) / &
                                             (cff_avg + ry_avg)
           enddo ; enddo
         endif
         if (segment%nudged_tan) then
           do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (ry_tang_obl(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_vel(i,j,k)
           enddo ; enddo
         endif
         if (segment%oblique_grad) then
           is_obc = max(segment%HI%IsdB,g%isd+1)
           ie_obc = min(segment%HI%IedB,g%ied-1)
           do k=1,nz ;  do i=segment%HI%IsdB+1,segment%HI%IedB-1
             rx_avg = rx_tang_obl(i,j,k)
             ry_avg = ry_tang_obl(i,j,k)
             cff_avg = cff_tangential(i,j,k)
             segment%tangential_grad(i,j,k) =  &
                 ((cff_avg*(u_new(i,j,k)   - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
                    ry_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) - &
                                  (max(rx_avg,0.0)*segment%grad_gradient(i,2,k) + &
                                   min(rx_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
                                 (cff_avg + ry_avg)
           enddo ; enddo
         endif
         if (segment%nudged_grad) then
           do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (ry_tang_obl(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_grad(i,j,k)
           enddo ; enddo
         endif
         deallocate(rx_tang_obl)
         deallocate(ry_tang_obl)
         deallocate(cff_tangential)
       endif
     endif
 
     if (segment%direction == obc_direction_s) then
       j=segment%HI%JsdB
       if (j>g%HI%JecB) cycle
       do k=1,nz ;  do i=segment%HI%isd,segment%HI%ied
         if (segment%radiation) then
           dhdt = (v_old(i,j+1,k) - v_new(i,j+1,k)) !old-new
           dhdy = (v_new(i,j+1,k) - v_new(i,j+2,k)) !in new time backward sasha for J-1
           ry_new = 0.0
           if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
           if (gamma_u < 1.0) then
             ry_avg = (1.0-gamma_u)*segment%ry_norm_rad(i,j,k) + gamma_u*ry_new
           else
             ry_avg = ry_new
           endif
           segment%ry_norm_rad(i,j,k) = ry_avg
           ! The new boundary value is interpolated between future interior
           ! value, v_new(J+1) and past boundary value but with barotropic
           ! accelerations, v_new(J).
           segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) / (1.0+ry_avg)
           if (gamma_u < 1.0) then
             ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
             ! implemented as a work-around to limitations in restart capability
             obc%ry_normal(i,j,k) = segment%ry_norm_rad(i,j,k)
           endif
         elseif (segment%oblique) then
           dhdt = (v_old(i,j+1,k) - v_new(i,j+1,k)) !old-new
           dhdy = (v_new(i,j+1,k) - v_new(i,j+2,k)) !in new time backward sasha for J-1
           if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0) then
             dhdx = segment%grad_normal(i-1,1,k)
           elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0) then
             dhdx = 0.0
           else
             dhdx = segment%grad_normal(i,1,k)
           endif
           if (dhdt*dhdy < 0.0) dhdt = 0.0
 
           ry_new = dhdt*dhdy
           cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
           rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
           if (gamma_u < 1.0) then
             rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
             ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
             cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
           else
             rx_avg = rx_new
             ry_avg = ry_new
             cff_avg = cff_new
           endif
           segment%rx_norm_obl(i,j,k) = rx_avg
           segment%ry_norm_obl(i,j,k) = ry_avg
           segment%cff_normal(i,j,k) = cff_avg
           segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) - &
                                        (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) + &
                                         min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
                                       (cff_avg + ry_avg)
           if (gamma_u < 1.0) then
             ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
             ! implemented as a work-around to limitations in restart capability
             obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
             obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
             obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
           endif
         elseif (segment%gradient) then
           segment%normal_vel(i,j,k) = v_new(i,j+1,k)
         endif
         if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then
           ! dhdt gets set to 0 on inflow in oblique case
           if (dhdt*dhdy <= 0.0) then
             tau = segment%Velocity_nudging_timescale_in
           else
             tau = segment%Velocity_nudging_timescale_out
           endif
           gamma_2 = dt / (tau + dt)
           segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_normal_vel(i,j,k)
         endif
       enddo ; enddo
       if (segment%radiation_tan .or. segment%radiation_grad) then
         j=segment%HI%JsdB
         allocate(ry_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         do k=1,nz
           if (gamma_u < 1.0) then
             ry_tang_rad(segment%HI%IsdB,j,k) = segment%ry_norm_rad(segment%HI%isd,j,k)
             ry_tang_rad(segment%HI%IedB,j,k) = segment%ry_norm_rad(segment%HI%ied,j,k)
             do i=segment%HI%IsdB+1,segment%HI%IedB-1
               ry_tang_rad(i,j,k) = 0.5*(segment%ry_norm_rad(i,j,k) + segment%ry_norm_rad(i+1,j,k))
             enddo
           else
             do i=segment%HI%IsdB,segment%HI%IedB
               dhdt = u_old(i,j+1,k)-u_new(i,j+1,k)   !old-new
               dhdy = u_new(i,j+1,k)-u_new(i,j+2,k) !in new time backward sasha for I-1
               ry_tang_rad(i,j,k) = 0.0
               if (dhdt*dhdy > 0.0) ry_tang_rad(i,j,k) = min( (dhdt/dhdy), rx_max) ! outward phase speed
             enddo
           endif
         enddo
         if (segment%radiation_tan) then
           do k=1,nz ;  do i=segment%HI%IsdB,segment%HI%IedB
             ry_avg = ry_tang_rad(i,j,k)
             segment%tangential_vel(i,j,k) = (u_new(i,j+1,k) + ry_avg*u_new(i,j+2,k)) / (1.0+ry_avg)
           enddo ; enddo
         endif
         if (segment%nudged_tan) then
           do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (ry_tang_rad(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_vel(i,j,k)
           enddo ; enddo
         endif
         if (segment%radiation_grad) then
           is_obc = max(segment%HI%IsdB,g%isd+1)
           ie_obc = min(segment%HI%IedB,g%ied-1)
           do k=1,nz ;  do i=is_obc,ie_obc
             ry_avg = ry_tang_rad(i,j,k)
!            if (G%mask2dCv(i,J+1) > 0.0 .and. G%mask2dCv(i+1,J+1) > 0.0) then
!              ry_avg = 0.5*(v_new(i,J+1,k) + v_new(i+1,J+1,k)) * dt * G%IdyBu(I,J+1)
!            elseif (G%mask2dCv(i,J+1) > 0.0) then
!              ry_avg = v_new(i,J+1,k) * dt * G%IdyBu(I,J+1)
!            elseif (G%mask2dCv(i+1,J+1) > 0.0) then
!              ry_avg = v_new(i+1,J+1,k) * dt * G%IdyBu(I,J+1)
!            else
!              ry_avg = 0.0
!            endif
             segment%tangential_grad(i,j,k) = ((u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
                               ry_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) / (1.0+ry_avg)
           enddo ; enddo
         endif
         if (segment%nudged_grad) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (ry_tang_rad(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_grad(i,j,k)
           enddo ; enddo
         endif
         deallocate(ry_tang_rad)
       endif
       if (segment%oblique_tan .or. segment%oblique_grad) then
         j=segment%HI%JsdB
         allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
         do k=1,nz
           if (gamma_u < 1.0) then
             rx_tang_obl(segment%HI%IsdB,j,k) = segment%rx_norm_obl(segment%HI%isd,j,k)
             rx_tang_obl(segment%HI%IedB,j,k) = segment%rx_norm_obl(segment%HI%ied,j,k)
             ry_tang_obl(segment%HI%IsdB,j,k) = segment%ry_norm_obl(segment%HI%isd,j,k)
             ry_tang_obl(segment%HI%IedB,j,k) = segment%ry_norm_obl(segment%HI%ied,j,k)
             cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
             cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
             do i=segment%HI%IsdB+1,segment%HI%IedB-1
               rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i+1,j,k))
               ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i+1,j,k))
               cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
             enddo
           else
             do i=segment%HI%IsdB,segment%HI%IedB
               dhdt = u_old(i,j+1,k)-u_new(i,j+1,k)   !old-new
               dhdy = u_new(i,j+1,k)-u_new(i,j+2,k) !in new time backward sasha for I-1
               if (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) > 0.0) then
                 dhdx = segment%grad_tan(i,1,k)
               elseif (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) == 0.0) then
                 dhdx = 0.0
               else
                 dhdx = segment%grad_tan(i+1,1,k)
               endif
               if (dhdt*dhdy < 0.0) dhdt = 0.0
               ry_new = dhdt*dhdy
               cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
               rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
               rx_tang_obl(i,j,k) = rx_new
               ry_tang_obl(i,j,k) = ry_new
               cff_tangential(i,j,k) = cff_new
             enddo
           endif
         enddo
         if (segment%oblique_tan) then
           do k=1,nz ;  do i=segment%HI%IsdB,segment%HI%IedB
             rx_avg = rx_tang_obl(i,j,k)
             ry_avg = ry_tang_obl(i,j,k)
             cff_avg = cff_tangential(i,j,k)
             segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j+1,k) + ry_avg*u_new(i,j+2,k)) - &
                                              (max(rx_avg,0.0)*segment%grad_tan(i,2,k) + &
                                               min(rx_avg,0.0)*segment%grad_tan(i+1,2,k)) ) / &
                                             (cff_avg + ry_avg)
           enddo ; enddo
         endif
         if (segment%nudged_tan) then
           do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (ry_tang_obl(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
                                             gamma_2 * segment%nudged_tangential_vel(i,j,k)
           enddo ; enddo
         endif
         if (segment%oblique_grad) then
           is_obc = max(segment%HI%IsdB,g%isd+1)
           ie_obc = min(segment%HI%IedB,g%ied-1)
           do k=1,nz ;  do i=segment%HI%IsdB+1,segment%HI%IedB-1
             rx_avg = rx_tang_obl(i,j,k)
             ry_avg = ry_tang_obl(i,j,k)
             cff_avg = cff_tangential(i,j,k)
             segment%tangential_grad(i,j,k) = &
                 ((cff_avg*(u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
                    ry_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) - &
                                  (max(rx_avg,0.0)*segment%grad_gradient(i,2,k) + &
                                   min(rx_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
                                 (cff_avg + ry_avg)
           enddo ; enddo
         endif
         if (segment%nudged_grad) then
           do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
             ! dhdt gets set to 0 on inflow in oblique case
             if (ry_tang_obl(i,j,k) <= 0.0) then
               tau = segment%Velocity_nudging_timescale_in
             else
               tau = segment%Velocity_nudging_timescale_out
             endif
             gamma_2 = dt / (tau + dt)
             segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
                                 gamma_2 * segment%nudged_tangential_grad(i,j,k)
           enddo ; enddo
         endif
         deallocate(rx_tang_obl)
         deallocate(ry_tang_obl)
         deallocate(cff_tangential)
       endif
     endif
  enddo
 
  ! Actually update u_new, v_new
  call open_boundary_apply_normal_flow(obc, g, u_new, v_new)
 
  call pass_vector(u_new, v_new, g%Domain, clock=id_clock_pass)
 
end subroutine radiation_open_bdry_conds
 
!> Applies OBC values stored in segments to 3d u,v fields
subroutine open_boundary_apply_normal_flow(OBC, G, u, v)
  ! Arguments
  type(ocean_obc_type),                      pointer       :: obc !< Open boundary control structure
  type(ocean_grid_type),                     intent(inout) :: g   !< Ocean grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u   !< u field to update on open
                                                                  !! boundaries [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v   !< v field to update on open
                                                                  !! boundaries [L T-1 ~> m s-1]
  ! Local variables
  integer :: i, j, k, n
  type(obc_segment_type), pointer :: segment => null()
 
  if (.not.associated(obc)) return ! Bail out if OBC is not available
 
  do n=1,obc%number_of_segments
    segment => obc%segment(n)
    if (.not. segment%on_pe) then
      cycle
    elseif (segment%radiation .or. segment%oblique .or. segment%gradient) then
      if (segment%is_E_or_W) then
        i=segment%HI%IsdB
        do k=1,g%ke ;  do j=segment%HI%jsd,segment%HI%jed
          u(i,j,k) = segment%normal_vel(i,j,k)
        enddo ; enddo
      elseif (segment%is_N_or_S) then
        j=segment%HI%JsdB
        do k=1,g%ke ;  do i=segment%HI%isd,segment%HI%ied
          v(i,j,k) = segment%normal_vel(i,j,k)
        enddo ; enddo
      endif
    endif
  enddo
 
end subroutine open_boundary_apply_normal_flow
 
!> Applies zero values to 3d u,v fields on OBC segments
subroutine open_boundary_zero_normal_flow(OBC, G, u, v)
  ! Arguments
  type(ocean_obc_type),                      pointer       :: obc !< Open boundary control structure
  type(ocean_grid_type),                     intent(inout) :: g   !< Ocean grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u   !< u field to update on open boundaries
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v   !< v field to update on open boundaries
  ! Local variables
  integer :: i, j, k, n
  type(obc_segment_type), pointer :: segment => null()
 
  if (.not.associated(obc)) return ! Bail out if OBC is not available
 
  do n=1,obc%number_of_segments
    segment => obc%segment(n)
    if (.not. segment%on_pe) then
      cycle
    elseif (segment%is_E_or_W) then
      i=segment%HI%IsdB
      do k=1,g%ke ;  do j=segment%HI%jsd,segment%HI%jed
        u(i,j,k) = 0.
      enddo ; enddo
    elseif (segment%is_N_or_S) then
      j=segment%HI%JsdB
      do k=1,g%ke ;  do i=segment%HI%isd,segment%HI%ied
        v(i,j,k) = 0.
      enddo ; enddo
    endif
  enddo
 
end subroutine open_boundary_zero_normal_flow
 
!> Calculate the tangential gradient of the normal flow at the boundary q-points.
subroutine gradient_at_q_points(G, segment, uvel, vvel)
  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
  type(obc_segment_type), pointer :: segment !< OBC segment structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: uvel !< zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: vvel !< meridional velocity [L T-1 ~> m s-1]
  integer :: i,j,k
 
  if (.not. segment%on_pe) return
 
  if (segment%is_E_or_W) then
    if (segment%direction == obc_direction_e) then
      i=segment%HI%isdB
      do k=1,g%ke
        do j=max(segment%HI%JsdB, g%HI%JsdB+1),min(segment%HI%JedB, g%HI%JedB-1)
          segment%grad_normal(j,1,k) = (uvel(i-1,j+1,k)-uvel(i-1,j,k)) * g%mask2dBu(i-1,j)
          segment%grad_normal(j,2,k) = (uvel(i,j+1,k)-uvel(i,j,k)) * g%mask2dBu(i,j)
        enddo
      enddo
      if (segment%oblique_tan) then
        do k=1,g%ke
          do j=max(segment%HI%jsd-1, g%HI%jsd),min(segment%HI%jed+1, g%HI%jed)
            segment%grad_tan(j,1,k) = (vvel(i-1,j,k)-vvel(i-1,j-1,k)) * g%mask2dT(i-1,j)
            segment%grad_tan(j,2,k) = (vvel(i,j,k)-vvel(i,j-1,k)) * g%mask2dT(i,j)
          enddo
        enddo
      endif
      if (segment%oblique_grad) then
        do k=1,g%ke
          do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
            segment%grad_gradient(j,1,k) = (((vvel(i-1,j,k) - vvel(i-2,j,k))*g%IdxBu(i-2,j)) - &
                 (vvel(i-1,j-1,k) - vvel(i-2,j-1,k))*g%IdxBu(i-2,j-1)) * g%mask2dCu(i-2,j)
            segment%grad_gradient(j,2,k) = (((vvel(i,j,k) - vvel(i-1,j,k))*g%IdxBu(i-1,j)) - &
                 (vvel(i,j-1,k) - vvel(i-1,j-1,k))*g%IdxBu(i-1,j-1)) * g%mask2dCu(i-1,j)
          enddo
        enddo
      endif
    else ! western segment
      i=segment%HI%isdB
      do k=1,g%ke
        do j=max(segment%HI%JsdB, g%HI%JsdB+1),min(segment%HI%JedB, g%HI%JedB-1)
          segment%grad_normal(j,1,k) = (uvel(i+1,j+1,k)-uvel(i+1,j,k)) * g%mask2dBu(i+1,j)
          segment%grad_normal(j,2,k) = (uvel(i,j+1,k)-uvel(i,j,k)) * g%mask2dBu(i,j)
        enddo
      enddo
      if (segment%oblique_tan) then
        do k=1,g%ke
          do j=max(segment%HI%jsd-1, g%HI%jsd),min(segment%HI%jed+1, g%HI%jed)
            segment%grad_tan(j,1,k) = (vvel(i+2,j,k)-vvel(i+2,j-1,k)) * g%mask2dT(i+2,j)
            segment%grad_tan(j,2,k) = (vvel(i+1,j,k)-vvel(i+1,j-1,k)) * g%mask2dT(i+1,j)
          enddo
        enddo
      endif
      if (segment%oblique_grad) then
        do k=1,g%ke
          do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
            segment%grad_gradient(j,1,k) = (((vvel(i+3,j,k) - vvel(i+2,j,k))*g%IdxBu(i+2,j)) - &
                 (vvel(i+3,j-1,k) - vvel(i+2,j-1,k))*g%IdxBu(i+2,j-1)) * g%mask2dCu(i+2,j)
            segment%grad_gradient(j,2,k) = (((vvel(i+2,j,k) - vvel(i+1,j,k))*g%IdxBu(i+1,j)) - &
                 (vvel(i+2,j-1,k) - vvel(i+1,j-1,k))*g%IdxBu(i+1,j-1)) * g%mask2dCu(i+1,j)
          enddo
        enddo
      endif
    endif
  elseif (segment%is_N_or_S) then
    if (segment%direction == obc_direction_n) then
      j=segment%HI%jsdB
      do k=1,g%ke
        do i=max(segment%HI%IsdB, g%HI%IsdB+1),min(segment%HI%IedB, g%HI%IedB-1)
          segment%grad_normal(i,1,k) = (vvel(i+1,j-1,k)-vvel(i,j-1,k)) * g%mask2dBu(i,j-1)
          segment%grad_normal(i,2,k) = (vvel(i+1,j,k)-vvel(i,j,k)) * g%mask2dBu(i,j)
        enddo
      enddo
      if (segment%oblique_tan) then
        do k=1,g%ke
          do i=max(segment%HI%isd-1, g%HI%isd),min(segment%HI%ied+1, g%HI%ied)
            segment%grad_tan(i,1,k) = (uvel(i,j-1,k)-uvel(i-1,j-1,k)) * g%mask2dT(i,j-1)
            segment%grad_tan(i,2,k) = (uvel(i,j,k)-uvel(i-1,j,k)) * g%mask2dT(i,j)
          enddo
        enddo
      endif
      if (segment%oblique_grad) then
        do k=1,g%ke
          do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
            segment%grad_gradient(i,1,k) = (((uvel(i,j-1,k) - uvel(i,j-2,k))*g%IdyBu(i,j-2)) - &
                 (uvel(i-1,j-1,k) - uvel(i-1,j-2,k))*g%IdyBu(i-1,j-2)) * g%mask2dCv(i,j-2)
            segment%grad_gradient(i,2,k) = (((uvel(i,j,k) - uvel(i,j-1,k))*g%IdyBu(i,j-1)) - &
                 (uvel(i-1,j,k) - uvel(i-1,j-1,k))*g%IdyBu(i-1,j-1)) * g%mask2dCv(i,j-1)
          enddo
        enddo
      endif
    else ! south segment
      j=segment%HI%jsdB
      do k=1,g%ke
        do i=max(segment%HI%IsdB, g%HI%IsdB+1),min(segment%HI%IedB, g%HI%IedB-1)
          segment%grad_normal(i,1,k) = (vvel(i+1,j+1,k)-vvel(i,j+1,k)) * g%mask2dBu(i,j+1)
          segment%grad_normal(i,2,k) = (vvel(i+1,j,k)-vvel(i,j,k)) * g%mask2dBu(i,j)
        enddo
      enddo
      if (segment%oblique_tan) then
        do k=1,g%ke
          do i=max(segment%HI%isd-1, g%HI%isd),min(segment%HI%ied+1, g%HI%ied)
            segment%grad_tan(i,1,k) = (uvel(i,j+2,k)-uvel(i-1,j+2,k)) * g%mask2dT(i,j+2)
            segment%grad_tan(i,2,k) = (uvel(i,j+1,k)-uvel(i-1,j+1,k)) * g%mask2dT(i,j+1)
          enddo
        enddo
      endif
      if (segment%oblique_grad) then
        do k=1,g%ke
          do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
            segment%grad_gradient(i,1,k) = (((uvel(i,j+3,k) - uvel(i,j+2,k))*g%IdyBu(i,j+2)) - &
                 (uvel(i-1,j+3,k) - uvel(i-1,j+2,k))*g%IdyBu(i-1,j+2)) * g%mask2dCv(i,j+2)
            segment%grad_gradient(i,2,k) = (((uvel(i,j+2,k) - uvel(i,j+1,k))*g%IdyBu(i,j+1)) - &
                 (uvel(i-1,j+2,k) - uvel(i-1,j+1,k))*g%IdyBu(i-1,j+1)) * g%mask2dCv(i,j+1)
          enddo
        enddo
      endif
    endif
  endif
 
end subroutine gradient_at_q_points
 
 
!> Sets the initial values of the tracer open boundary conditions.
!! Redoing this elsewhere.
subroutine set_tracer_data(OBC, tv, h, G, PF, tracer_Reg)
  type(ocean_grid_type),                     intent(inout) :: g !< Ocean grid structure
  type(ocean_obc_type),                      pointer       :: obc !< Open boundary structure
  type(thermo_var_ptrs),                     intent(inout) :: tv !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h !< Thickness
  type(param_file_type),                     intent(in)    :: pf !< Parameter file handle
  type(tracer_registry_type),                pointer       :: tracer_reg !< Tracer registry
  ! Local variables
  integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, n
  integer :: isd_off, jsd_off
  integer :: isdb, iedb, jsdb, jedb
  type(obc_segment_type), pointer :: segment  => null() ! pointer to segment type list
  character(len=40)  :: mdl = "set_tracer_data" ! This subroutine's name.
  character(len=200) :: filename, obc_file, inputdir ! Strings for file/path
 
  real :: temp_u(g%domain%niglobal+1,g%domain%njglobal)
  real :: temp_v(g%domain%niglobal,g%domain%njglobal+1)
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
 
  ! For now, there are no radiation conditions applied to the thicknesses, since
  ! the thicknesses might not be physically motivated.  Instead, sponges should be
  ! used to enforce the near-boundary layer structure.
 
  if (associated(tv%T)) then
 
    call pass_var(tv%T, g%Domain)
    call pass_var(tv%S, g%Domain)
 
    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 k=1,g%ke ;  do j=segment%HI%jsd,segment%HI%jed
          tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k)
        enddo ; enddo
      elseif (segment%direction == obc_direction_w) then
        i=segment%HI%IsdB
        do k=1,g%ke ;  do j=segment%HI%jsd,segment%HI%jed
          tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k)
        enddo ; enddo
      elseif (segment%direction == obc_direction_n) then
        j=segment%HI%JsdB
        do k=1,g%ke ;  do i=segment%HI%isd,segment%HI%ied
          tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k)
        enddo ; enddo
      elseif (segment%direction == obc_direction_s) then
        j=segment%HI%JsdB
        do k=1,g%ke ;  do i=segment%HI%isd,segment%HI%ied
          tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k)
        enddo ; enddo
      endif
    enddo
  endif
 
end subroutine set_tracer_data
 
!> Needs documentation
function lookup_seg_field(OBC_seg,field)
  type(obc_segment_type), pointer :: obc_seg !< OBC segment
  character(len=32), intent(in) :: field !< The field name
  integer :: lookup_seg_field
  ! Local variables
  integer :: n
 
  lookup_seg_field=-1
  do n=1,obc_seg%num_fields
   if (trim(field) == obc_seg%field(n)%name) then
     lookup_seg_field=n
     return
   endif
  enddo
 
end function lookup_seg_field
 
 
!> Allocate segment data fields
subroutine allocate_obc_segment_data(OBC, segment)
  type(ocean_obc_type),   pointer       :: OBC     !< Open boundary structure
  type(obc_segment_type), intent(inout) :: segment !< Open boundary segment
  ! Local variables
  integer :: isd, ied, jsd, jed
  integer :: IsdB, IedB, JsdB, JedB
  integer :: IscB, IecB, JscB, JecB
  character(len=40)  :: mdl = "allocate_OBC_segment_data" ! This subroutine's name.
 
  isd = segment%HI%isd ; ied = segment%HI%ied
  jsd = segment%HI%jsd ; jed = segment%HI%jed
  isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
  jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
  iscb = segment%HI%IscB ; iecb = segment%HI%IecB
  jscb = segment%HI%JscB ; jecb = segment%HI%JecB
 
  if (.not. segment%on_pe) return
 
  if (segment%is_E_or_W) then
    ! If these are just Flather, change update_OBC_segment_data accordingly
    allocate(segment%Cg(isdb:iedb,jsd:jed));                  segment%Cg(:,:)=0.
    allocate(segment%Htot(isdb:iedb,jsd:jed));                segment%Htot(:,:)=0.0
    allocate(segment%h(isdb:iedb,jsd:jed,obc%ke));            segment%h(:,:,:)=0.0
    allocate(segment%eta(isdb:iedb,jsd:jed));                 segment%eta(:,:)=0.0
    if (segment%radiation) then
      allocate(segment%rx_norm_rad(isdb:iedb,jsd:jed,obc%ke));  segment%rx_norm_rad(:,:,:)=0.0
    endif
    allocate(segment%normal_vel(isdb:iedb,jsd:jed,obc%ke));   segment%normal_vel(:,:,:)=0.0
    allocate(segment%normal_vel_bt(isdb:iedb,jsd:jed));       segment%normal_vel_bt(:,:)=0.0
    allocate(segment%normal_trans(isdb:iedb,jsd:jed,obc%ke)); segment%normal_trans(:,:,:)=0.0
    if (segment%nudged) then
      allocate(segment%nudged_normal_vel(isdb:iedb,jsd:jed,obc%ke)); segment%nudged_normal_vel(:,:,:)=0.0
    endif
    if (segment%radiation_tan .or. segment%nudged_tan .or. segment%specified_tan .or. &
        segment%oblique_tan .or. obc%computed_vorticity .or. obc%computed_strain) then
      allocate(segment%tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_vel(:,:,:)=0.0
    endif
    if (segment%nudged_tan) then
      allocate(segment%nudged_tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_vel(:,:,:)=0.0
    endif
    if (segment%nudged_grad) then
      allocate(segment%nudged_tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
    endif
    if (obc%specified_vorticity .or. obc%specified_strain .or. segment%radiation_grad .or. &
              segment%oblique_grad) then
      allocate(segment%tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_grad(:,:,:)=0.0
    endif
    if (segment%oblique) then
      allocate(segment%grad_normal(jsdb:jedb,2,obc%ke));      segment%grad_normal(:,:,:) = 0.0
      allocate(segment%rx_norm_obl(isdb:iedb,jsd:jed,obc%ke));  segment%rx_norm_obl(:,:,:)=0.0
      allocate(segment%ry_norm_obl(isdb:iedb,jsd:jed,obc%ke));  segment%ry_norm_obl(:,:,:)=0.0
      allocate(segment%cff_normal(isdb:iedb,jsd:jed,obc%ke)); segment%cff_normal(:,:,:)=0.0
    endif
    if (segment%oblique_tan) then
      allocate(segment%grad_tan(jsd-1:jed+1,2,obc%ke));           segment%grad_tan(:,:,:) = 0.0
    endif
    if (segment%oblique_grad) then
      allocate(segment%grad_gradient(jsd:jed,2,obc%ke));      segment%grad_gradient(:,:,:) = 0.0
    endif
  endif
 
  if (segment%is_N_or_S) then
    ! If these are just Flather, change update_OBC_segment_data accordingly
    allocate(segment%Cg(isd:ied,jsdb:jedb));                  segment%Cg(:,:)=0.
    allocate(segment%Htot(isd:ied,jsdb:jedb));                segment%Htot(:,:)=0.0
    allocate(segment%h(isd:ied,jsdb:jedb,obc%ke));            segment%h(:,:,:)=0.0
    allocate(segment%eta(isd:ied,jsdb:jedb));                 segment%eta(:,:)=0.0
    if (segment%radiation) then
      allocate(segment%ry_norm_rad(isd:ied,jsdb:jedb,obc%ke));  segment%ry_norm_rad(:,:,:)=0.0
    endif
    allocate(segment%normal_vel(isd:ied,jsdb:jedb,obc%ke));   segment%normal_vel(:,:,:)=0.0
    allocate(segment%normal_vel_bt(isd:ied,jsdb:jedb));       segment%normal_vel_bt(:,:)=0.0
    allocate(segment%normal_trans(isd:ied,jsdb:jedb,obc%ke)); segment%normal_trans(:,:,:)=0.0
    if (segment%nudged) then
      allocate(segment%nudged_normal_vel(isd:ied,jsdb:jedb,obc%ke)); segment%nudged_normal_vel(:,:,:)=0.0
    endif
    if (segment%radiation_tan .or. segment%nudged_tan .or. segment%specified_tan .or. &
        segment%oblique_tan .or. obc%computed_vorticity .or. obc%computed_strain) then
      allocate(segment%tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_vel(:,:,:)=0.0
    endif
    if (segment%nudged_tan) then
      allocate(segment%nudged_tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_vel(:,:,:)=0.0
    endif
    if (segment%nudged_grad) then
      allocate(segment%nudged_tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
    endif
    if (obc%specified_vorticity .or. obc%specified_strain .or. segment%radiation_grad .or. &
              segment%oblique_grad) then
      allocate(segment%tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_grad(:,:,:)=0.0
    endif
    if (segment%oblique) then
      allocate(segment%grad_normal(isdb:iedb,2,obc%ke));      segment%grad_normal(:,:,:) = 0.0
      allocate(segment%rx_norm_obl(isd:ied,jsdb:jedb,obc%ke));  segment%rx_norm_obl(:,:,:)=0.0
      allocate(segment%ry_norm_obl(isd:ied,jsdb:jedb,obc%ke));  segment%ry_norm_obl(:,:,:)=0.0
      allocate(segment%cff_normal(isd:ied,jsdb:jedb,obc%ke)); segment%cff_normal(:,:,:)=0.0
    endif
    if (segment%oblique_tan) then
      allocate(segment%grad_tan(isd-1:ied+1,2,obc%ke));           segment%grad_tan(:,:,:) = 0.0
    endif
    if (segment%oblique_grad) then
      allocate(segment%grad_gradient(isd:ied,2,obc%ke));      segment%grad_gradient(:,:,:) = 0.0
    endif
  endif
 
end subroutine allocate_obc_segment_data
 
!> Deallocate segment data fields
subroutine deallocate_obc_segment_data(OBC, segment)
  type(ocean_obc_type),   pointer       :: OBC     !< Open boundary structure
  type(obc_segment_type), intent(inout) :: segment !< Open boundary segment
  ! Local variables
  character(len=40)  :: mdl = "deallocate_OBC_segment_data" ! This subroutine's name.
 
  if (.not. segment%on_pe) return
 
  if (associated (segment%Cg)) deallocate(segment%Cg)
  if (associated (segment%Htot)) deallocate(segment%Htot)
  if (associated (segment%h)) deallocate(segment%h)
  if (associated (segment%eta)) deallocate(segment%eta)
  if (associated (segment%rx_norm_rad)) deallocate(segment%rx_norm_rad)
  if (associated (segment%ry_norm_rad)) deallocate(segment%ry_norm_rad)
  if (associated (segment%rx_norm_obl)) deallocate(segment%rx_norm_obl)
  if (associated (segment%ry_norm_obl)) deallocate(segment%ry_norm_obl)
  if (associated (segment%cff_normal)) deallocate(segment%cff_normal)
  if (associated (segment%grad_normal)) deallocate(segment%grad_normal)
  if (associated (segment%grad_tan)) deallocate(segment%grad_tan)
  if (associated (segment%grad_gradient)) deallocate(segment%grad_gradient)
  if (associated (segment%normal_vel)) deallocate(segment%normal_vel)
  if (associated (segment%normal_vel_bt)) deallocate(segment%normal_vel_bt)
  if (associated (segment%normal_trans)) deallocate(segment%normal_trans)
  if (associated (segment%nudged_normal_vel)) deallocate(segment%nudged_normal_vel)
  if (associated (segment%tangential_vel)) deallocate(segment%tangential_vel)
  if (associated (segment%nudged_tangential_vel)) deallocate(segment%nudged_tangential_vel)
  if (associated (segment%nudged_tangential_grad)) deallocate(segment%nudged_tangential_grad)
  if (associated (segment%tangential_grad)) deallocate(segment%tangential_grad)
  if (associated (segment%tr_Reg)) call segment_tracer_registry_end(segment%tr_Reg)
 
 
end subroutine deallocate_obc_segment_data
 
!> Set tangential velocities outside of open boundaries to silly values
!! (used for checking the interior state is independent of values outside
!! of the domain).
subroutine open_boundary_test_extern_uv(G, OBC, u, v)
  type(ocean_grid_type),                     intent(in)    :: g !< Ocean grid structure
  type(ocean_obc_type),                      pointer       :: obc !< Open boundary structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)),intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)),intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
  ! Local variables
  integer :: i, j, k, n
 
  if (.not. associated(obc)) return
 
  do n = 1, obc%number_of_segments
    do k = 1, g%ke
      if (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%IsdB, obc%segment(n)%HI%IedB
            u(i,j+1,k) = obc%silly_u
          enddo
        else
          do i = obc%segment(n)%HI%IsdB, obc%segment(n)%HI%IedB
            u(i,j,k) = obc%silly_u
          enddo
        endif
      elseif (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%JsdB, obc%segment(n)%HI%JedB
            v(i+1,j,k) = obc%silly_u
          enddo
        else
          do j = obc%segment(n)%HI%JsdB, obc%segment(n)%HI%JedB
            v(i,j,k) = obc%silly_u
          enddo
        endif
      endif
    enddo
  enddo
 
end subroutine open_boundary_test_extern_uv
 
!> Set thicknesses outside of open boundaries to silly values
!! (used for checking the interior state is independent of values outside
!! of the domain).
subroutine open_boundary_test_extern_h(G, GV, OBC, h)
  type(ocean_grid_type),                     intent(in)    :: g   !< Ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: gv  !<  Ocean vertical grid structure
  type(ocean_obc_type),                      pointer       :: obc !< Open boundary structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)),intent(inout) :: h   !< Layer thickness [H ~> m or kg m-2]
  ! Local variables
  real :: silly_h  ! A silly thickness for testing [H ~> m or kg m-2]
  integer :: i, j, k, n
 
  if (.not. associated(obc)) return
 
  silly_h = gv%Z_to_H*obc%silly_h
 
  do n = 1, obc%number_of_segments
    do k = 1, gv%ke
      if (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
            h(i,j+1,k) = silly_h
          enddo
        else
          do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
            h(i,j,k) = silly_h
          enddo
        endif
      elseif (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
            h(i+1,j,k) = silly_h
          enddo
        else
          do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
            h(i,j,k) = silly_h
          enddo
        endif
      endif
    enddo
  enddo
 
end subroutine open_boundary_test_extern_h
 
!> Update the OBC values on the segments.
subroutine update_obc_segment_data(G, GV, US, OBC, tv, h, Time)
  type(ocean_grid_type),                     intent(in)    :: g    !< Ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: gv   !<  Ocean vertical grid structure
  type(unit_scale_type),                     intent(in)    :: us   !< A dimensional unit scaling type
  type(ocean_obc_type),                      pointer       :: obc  !< Open boundary structure
  type(thermo_var_ptrs),                     intent(in)    :: tv   !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h    !< Thickness [m]
  type(time_type),                           intent(in)    :: time !< Model time
  ! Local variables
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed
  integer :: isdb, iedb, jsdb, jedb, n, m, nz
  character(len=40)  :: mdl = "set_OBC_segment_data" ! This subroutine's name.
  character(len=200) :: filename, obc_file, inputdir ! Strings for file/path
  type(obc_segment_type), pointer :: segment => null()
  integer, dimension(4) :: siz,siz2
  real :: sumh ! column sum of thicknesses [m]
  integer :: ni_seg, nj_seg  ! number of src gridpoints along the segments
  integer :: i2, j2          ! indices for referencing local domain array
  integer :: is_obc, ie_obc, js_obc, je_obc  ! segment indices within local domain
  integer :: ishift, jshift  ! offsets for staggered locations
  real, dimension(:,:), pointer :: seg_vel => null()  ! pointer to segment velocity array
  real, dimension(:,:), pointer :: seg_trans => null()  ! pointer to segment transport array
  real, dimension(:,:,:), allocatable :: tmp_buffer
  real, dimension(:), allocatable :: h_stack
  integer :: is_obc2, js_obc2
  real :: net_h_src, net_h_int, scl_fac
  real, pointer, dimension(:,:)   :: normal_trans_bt=>null() ! barotropic transport
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  nz=g%ke
 
  if (.not. associated(obc)) return
 
  do n = 1, obc%number_of_segments
    segment => obc%segment(n)
 
    if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain
 
    ni_seg = segment%ie_obc-segment%is_obc+1
    nj_seg = segment%je_obc-segment%js_obc+1
    is_obc = max(segment%is_obc,isd-1)
    ie_obc = min(segment%ie_obc,ied)
    js_obc = max(segment%js_obc,jsd-1)
    je_obc = min(segment%je_obc,jed)
 
! Calculate auxiliary fields at staggered locations.
! Segment indices are on q points:
!
!       |-----------|------------|-----------|-----------|  J_obc
!     Is_obc                                          Ie_obc
!
! i2 has to start at Is_obc+1 and end at Ie_obc.
! j2 is J_obc and jshift has to be +1 at both the north and south.
 
     ! calculate auxiliary fields at staggered locations
    ishift=0;jshift=0
    if (segment%is_E_or_W) then
      allocate(normal_trans_bt(segment%HI%IsdB:segment%HI%IedB,segment%HI%jsd:segment%HI%jed))
      normal_trans_bt(:,:)=0.0
      if (segment%direction == obc_direction_w) ishift=1
      i=segment%HI%IsdB
      do j=segment%HI%jsd,segment%HI%jed
        segment%Cg(i,j) = sqrt(gv%g_prime(1)*g%bathyT(i+ishift,j))
        segment%Htot(i,j)=0.0
        do k=1,g%ke
          segment%h(i,j,k) = h(i+ishift,j,k)
          segment%Htot(i,j)=segment%Htot(i,j)+segment%h(i,j,k)
        enddo
      enddo
    else! (segment%direction == OBC_DIRECTION_N .or. segment%direction == OBC_DIRECTION_S)
      allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB))
      normal_trans_bt(:,:)=0.0
      if (segment%direction == obc_direction_s) jshift=1
      j=segment%HI%JsdB
      do i=segment%HI%isd,segment%HI%ied
        segment%Cg(i,j) = sqrt(gv%g_prime(1)*g%bathyT(i,j+jshift))
        segment%Htot(i,j)=0.0
        do k=1,g%ke
          segment%h(i,j,k) = h(i,j+jshift,k)
          segment%Htot(i,j)=segment%Htot(i,j)+segment%h(i,j,k)
        enddo
      enddo
    endif
 
    allocate(h_stack(g%ke))
    h_stack(:) = 0.0
    do m = 1,segment%num_fields
      if (segment%field(m)%fid > 0) then
        siz(1)=size(segment%field(m)%buffer_src,1)
        siz(2)=size(segment%field(m)%buffer_src,2)
        siz(3)=size(segment%field(m)%buffer_src,3)
        if (.not.associated(segment%field(m)%buffer_dst)) then
          if (siz(3) /= segment%field(m)%nk_src) call mom_error(fatal,'nk_src inconsistency')
          if (segment%field(m)%nk_src > 1) then
            if (segment%is_E_or_W) then
              if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
                allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
              else
                allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
              endif
              if (segment%field(m)%name == 'U') then
                allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
                segment%field(m)%bt_vel(:,:)=0.0
              endif
            else
              if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
                allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
              else
                allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
              endif
              if (segment%field(m)%name == 'V') then
                allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
                segment%field(m)%bt_vel(:,:)=0.0
              endif
            endif
          else
            if (segment%is_E_or_W) then
              if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
                allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
              else
                allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,1))
              endif
              if (segment%field(m)%name == 'U') then
                allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
                segment%field(m)%bt_vel(:,:)=0.0
              endif
            else
              if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
                allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
              else
                allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,1))
              endif
              if (segment%field(m)%name == 'V') then
                allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
                segment%field(m)%bt_vel(:,:)=0.0
              endif
            endif
          endif
          segment%field(m)%buffer_dst(:,:,:)=0.0
        endif
        ! read source data interpolated to the current model time
        if (siz(1)==1) then
          if (obc%brushcutter_mode) then
            allocate(tmp_buffer(1,nj_seg*2-1,segment%field(m)%nk_src))  ! segment data is currrently on supergrid
          else
            allocate(tmp_buffer(1,nj_seg,segment%field(m)%nk_src))  ! segment data is currrently on supergrid
          endif
        else
          if (obc%brushcutter_mode) then
            allocate(tmp_buffer(ni_seg*2-1,1,segment%field(m)%nk_src))  ! segment data is currrently on supergrid
          else
            allocate(tmp_buffer(ni_seg,1,segment%field(m)%nk_src))  ! segment data is currrently on supergrid
          endif
        endif
 
        call time_interp_external(segment%field(m)%fid,time, tmp_buffer)
        if (obc%brushcutter_mode) then
          if (segment%is_E_or_W) then
            if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
              segment%field(m)%buffer_src(is_obc,:,:) = &
                  tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset)+1:2,:)
            else
              segment%field(m)%buffer_src(is_obc,:,:) = &
                  tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset):2,:)
            endif
          else
            if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
              segment%field(m)%buffer_src(:,js_obc,:) = &
                  tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset)+1:2,1,:)
            else
              segment%field(m)%buffer_src(:,js_obc,:) = &
                  tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset):2,1,:)
            endif
          endif
        else
          if (segment%is_E_or_W) then
            if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
              segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset+1,:)
            else
              segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset,:)
            endif
          else
            if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
              segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset+1,1,:)
            else
              segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset,1,:)
            endif
          endif
        endif
        if (segment%field(m)%nk_src > 1) then
          call time_interp_external(segment%field(m)%fid_dz,time, tmp_buffer)
          if (obc%brushcutter_mode) then
            if (segment%is_E_or_W) then
              if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
                segment%field(m)%dz_src(is_obc,:,:) = &
                    tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset)+1:2,:)
              else
                segment%field(m)%dz_src(is_obc,:,:) = &
                    tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset):2,:)
              endif
            else
              if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
                segment%field(m)%dz_src(:,js_obc,:) = &
                    tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset)+1:2,1,:)
              else
                segment%field(m)%dz_src(:,js_obc,:) = &
                    tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset):2,1,:)
              endif
            endif
          else
            if (segment%is_E_or_W) then
              if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
                segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset+1,:)
              else
                segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset,:)
              endif
            else
              if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
                segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset+1,1,:)
              else
                segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset,1,:)
              endif
            endif
          endif
 
          call adjustsegmentetatofitbathymetry(g,gv,us,segment,m)
 
          if (segment%is_E_or_W) then
            ishift=1
            if (segment%direction == obc_direction_e) ishift=0
            i=is_obc
            if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
              ! Do q points for the whole segment
              do j=max(js_obc,jsd),min(je_obc,jed-1)
                ! Using the h remapping approach
                ! Pretty sure we need to check for source/target grid consistency here
                segment%field(m)%buffer_dst(i,j,:)=0.0  ! initialize remap destination buffer
                if (g%mask2dCu(i,j)>0. .and. g%mask2dCu(i,j+1)>0.) then
                  h_stack(:) = 0.5*(h(i+ishift,j,:) + h(i+ishift,j+1,:))
                  call remapping_core_h(obc%remap_CS, &
                       segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
                       segment%field(m)%buffer_src(i,j,:), &
                       g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
                elseif (g%mask2dCu(i,j)>0.) then
                  h_stack(:) = h(i+ishift,j,:)
                  call remapping_core_h(obc%remap_CS, &
                       segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
                       segment%field(m)%buffer_src(i,j,:), &
                       g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
                elseif (g%mask2dCu(i,j+1)>0.) then
                  h_stack(:) = h(i+ishift,j+1,:)
                  call remapping_core_h(obc%remap_CS, &
                       segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
                       segment%field(m)%buffer_src(i,j,:), &
                       g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
                endif
              enddo
            else
              do j=js_obc+1,je_obc
                ! Using the h remapping approach
                ! Pretty sure we need to check for source/target grid consistency here
                segment%field(m)%buffer_dst(i,j,:)=0.0  ! initialize remap destination buffer
                if (g%mask2dCu(i,j)>0.) then
                  net_h_src = sum( segment%field(m)%dz_src(i,j,:) )
                  net_h_int = sum( h(i+ishift,j,:) )
                  scl_fac = net_h_int / net_h_src
                  call remapping_core_h(obc%remap_CS, &
                       segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,j,:), &
                       segment%field(m)%buffer_src(i,j,:), &
                       g%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(i,j,:))
                endif
              enddo
            endif
          else
            jshift=1
            if (segment%direction == obc_direction_n) jshift=0
            j=js_obc
            if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
              ! Do q points for the whole segment
              do i=max(is_obc,isd),min(ie_obc,ied-1)
                segment%field(m)%buffer_dst(i,j,:)=0.0  ! initialize remap destination buffer
                if (g%mask2dCv(i,j)>0. .and. g%mask2dCv(i+1,j)>0.) then
              ! Using the h remapping approach
              ! Pretty sure we need to check for source/target grid consistency here
                  h_stack(:) = 0.5*(h(i,j+jshift,:) + h(i+1,j+jshift,:))
                  call remapping_core_h(obc%remap_CS, &
                       segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
                       segment%field(m)%buffer_src(i,j,:), &
                       g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
                elseif (g%mask2dCv(i,j)>0.) then
                  h_stack(:) = h(i,j+jshift,:)
                  call remapping_core_h(obc%remap_CS, &
                       segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
                       segment%field(m)%buffer_src(i,j,:), &
                       g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
                elseif (g%mask2dCv(i+1,j)>0.) then
                  h_stack(:) = h(i+1,j+jshift,:)
                  call remapping_core_h(obc%remap_CS, &
                       segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
                       segment%field(m)%buffer_src(i,j,:), &
                       g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
                endif
              enddo
            else
              do i=is_obc+1,ie_obc
              ! Using the h remapping approach
              ! Pretty sure we need to check for source/target grid consistency here
                segment%field(m)%buffer_dst(i,j,:)=0.0  ! initialize remap destination buffer
                if (g%mask2dCv(i,j)>0.) then
                  net_h_src = sum( segment%field(m)%dz_src(i,j,:) )
                  net_h_int = sum( h(i,j+jshift,:) )
                  scl_fac = net_h_int / net_h_src
                  call remapping_core_h(obc%remap_CS, &
                       segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
                       segment%field(m)%buffer_src(i,j,:), &
                       g%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,j,:))
                endif
              enddo
            endif
          endif
        else  ! 2d data
          segment%field(m)%buffer_dst(:,:,1) = segment%field(m)%buffer_src(:,:,1)  ! initialize remap destination buffer
        endif
        deallocate(tmp_buffer)
      else ! fid <= 0 (Uniform value)
        if (.not. associated(segment%field(m)%buffer_dst)) then
          if (segment%is_E_or_W) then
            if (segment%field(m)%name == 'V') then
              allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
              allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc))
            elseif (segment%field(m)%name == 'U') then
              allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
              allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
            elseif (segment%field(m)%name == 'DVDX') then
              allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
            elseif (segment%field(m)%name == 'SSH') then
              allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
            else
              allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
            endif
          else
            if (segment%field(m)%name == 'U') then
              allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
              allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc))
            elseif (segment%field(m)%name == 'V') then
              allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
              allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
            elseif (segment%field(m)%name == 'DUDY') then
              allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
            elseif (segment%field(m)%name == 'SSH') then
              allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
            else
              allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
            endif
          endif
          segment%field(m)%buffer_dst(:,:,:) = segment%field(m)%value
          if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then
            segment%field(m)%bt_vel(:,:) = segment%field(m)%value
          endif
        endif
      endif
 
      if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then
        if (segment%field(m)%fid>0) then ! calculate external BT velocity and transport if needed
          if (trim(segment%field(m)%name) == 'U' .and. segment%is_E_or_W) then
            i=is_obc
            do j=js_obc+1,je_obc
              normal_trans_bt(i,j) = 0.0
              do k=1,g%ke
                segment%normal_vel(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)
                segment%normal_trans(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)*segment%h(i,j,k) * &
                          g%dyCu(i,j)
                normal_trans_bt(i,j) = normal_trans_bt(i,j) + segment%normal_trans(i,j,k)
              enddo
              segment%normal_vel_bt(i,j) = normal_trans_bt(i,j) / (max(segment%Htot(i,j),1.e-12) * g%dyCu(i,j))
              if (associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,j,:) = segment%normal_vel(i,j,:)
            enddo
          elseif (trim(segment%field(m)%name) == 'V' .and. segment%is_N_or_S) then
            j=js_obc
            do i=is_obc+1,ie_obc
              normal_trans_bt(i,j) = 0.0
              do k=1,g%ke
                segment%normal_vel(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)
                segment%normal_trans(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)*segment%h(i,j,k) * &
                          g%dxCv(i,j)
                normal_trans_bt(i,j) = normal_trans_bt(i,j) + segment%normal_trans(i,j,k)
              enddo
              segment%normal_vel_bt(i,j) = normal_trans_bt(i,j) / (max(segment%Htot(i,j),1.e-12) * g%dxCv(i,j))
              if (associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,j,:) = segment%normal_vel(i,j,:)
            enddo
          elseif (trim(segment%field(m)%name) == 'V' .and. segment%is_E_or_W .and. &
                  associated(segment%tangential_vel)) then
            i=is_obc
            do j=js_obc,je_obc
              do k=1,g%ke
                segment%tangential_vel(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)
              enddo
              if (associated(segment%nudged_tangential_vel)) &
                segment%nudged_tangential_vel(i,j,:) = segment%tangential_vel(i,j,:)
            enddo
          elseif (trim(segment%field(m)%name) == 'U' .and. segment%is_N_or_S .and. &
                  associated(segment%tangential_vel)) then
            j=js_obc
            do i=is_obc,ie_obc
              do k=1,g%ke
                segment%tangential_vel(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)
              enddo
              if (associated(segment%nudged_tangential_vel)) &
                segment%nudged_tangential_vel(i,j,:) = segment%tangential_vel(i,j,:)
            enddo
          elseif (trim(segment%field(m)%name) == 'DVDX' .and. segment%is_E_or_W .and. &
                  associated(segment%tangential_grad)) then
            i=is_obc
            do j=js_obc,je_obc
              do k=1,g%ke
                segment%tangential_grad(i,j,k) = us%T_to_s*segment%field(m)%buffer_dst(i,j,k)
              enddo
            enddo
          elseif (trim(segment%field(m)%name) == 'DUDY' .and. segment%is_N_or_S .and. &
                  associated(segment%tangential_grad)) then
            j=js_obc
            do i=is_obc,ie_obc
              do k=1,g%ke
                segment%tangential_grad(i,j,k) = us%T_to_s*segment%field(m)%buffer_dst(i,j,k)
              enddo
            enddo
          endif
        endif
      endif
 
      ! from this point on, data are entirely on segments - will
      ! write all segment loops as 2d loops.
      if (segment%is_E_or_W) then
        js_obc2 = js_obc+1
        is_obc2 = is_obc
      else
        js_obc2 = js_obc
        is_obc2 = is_obc+1
      endif
      if (segment%is_N_or_S) then
        is_obc2 = is_obc+1
        js_obc2 = js_obc
      else
        is_obc2 = is_obc
        js_obc2 = js_obc+1
      endif
 
      if (trim(segment%field(m)%name) == 'SSH') then
        do j=js_obc2,je_obc
          do i=is_obc2,ie_obc
            segment%eta(i,j) = segment%field(m)%buffer_dst(i,j,1)
          enddo
        enddo
      endif
 
      if (trim(segment%field(m)%name) == 'TEMP') then
        if (associated(segment%field(m)%buffer_dst)) then
          do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc
            segment%tr_Reg%Tr(1)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
          enddo ; enddo ; enddo
          if (.not. segment%tr_Reg%Tr(1)%is_initialized) then
            ! if the tracer reservoir has not yet been initialized, then set to external value.
            do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc
              segment%tr_Reg%Tr(1)%tres(i,j,k) = segment%tr_Reg%Tr(1)%t(i,j,k)
            enddo ; enddo ; enddo
            segment%tr_Reg%Tr(1)%is_initialized=.true.
          endif
        else
          segment%tr_Reg%Tr(1)%OBC_inflow_conc = segment%field(m)%value
        endif
      elseif (trim(segment%field(m)%name) == 'SALT') then
        if (associated(segment%field(m)%buffer_dst)) then
          do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc
            segment%tr_Reg%Tr(2)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
          enddo ; enddo ; enddo
          if (.not. segment%tr_Reg%Tr(2)%is_initialized) then
            !if the tracer reservoir has not yet been initialized, then set to external value.
            do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc
              segment%tr_Reg%Tr(2)%tres(i,j,k) = segment%tr_Reg%Tr(2)%t(i,j,k)
            enddo ; enddo ; enddo
            segment%tr_Reg%Tr(2)%is_initialized=.true.
          endif
        else
          segment%tr_Reg%Tr(2)%OBC_inflow_conc = segment%field(m)%value
        endif
      endif
 
    enddo ! end field loop
    deallocate(h_stack)
    deallocate(normal_trans_bt)
 
  enddo ! end segment loop
 
end subroutine update_obc_segment_data
 
!> register open boundary objects for boundary updates.
subroutine register_obc(name, param_file, Reg)
  character(len=32),     intent(in)  :: name        !< OBC name used for error messages
  type(param_file_type), intent(in)  :: param_file  !< file to parse for  model parameter values
  type(obc_registry_type), pointer   :: reg         !< pointer to the tracer registry
  integer :: nobc
  character(len=256) :: mesg    ! Message for error messages.
 
  if (.not. associated(reg)) call obc_registry_init(param_file, reg)
 
  if (reg%nobc>=max_fields_) then
    write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for &
        &all the open boundaries being registered via register_OBC.")') reg%nobc+1
    call mom_error(fatal,"MOM register_tracer: "//mesg)
  endif
  reg%nobc = reg%nobc + 1
  nobc     = reg%nobc
 
  reg%OB(nobc)%name = name
 
  if (reg%locked) call mom_error(fatal, &
      "MOM register_tracer was called for variable "//trim(reg%OB(nobc)%name)//&
      " with a locked tracer registry.")
 
end subroutine register_obc
 
!> This routine include declares and sets the variable "version".
subroutine obc_registry_init(param_file, Reg)
  type(param_file_type),   intent(in) :: param_file !< open file to parse for model parameters
  type(obc_registry_type), pointer    :: reg        !< pointer to OBC registry
 
  integer, save :: init_calls = 0
 
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_open_boundary" ! This module's name.
  character(len=256) :: mesg    ! Message for error messages.
 
  if (.not.associated(reg)) then ; allocate(reg)
  else ; return ; endif
 
  ! Read all relevant parameters and write them to the model log.
! call log_version(param_file, mdl,s version, "")
 
  init_calls = init_calls + 1
  if (init_calls > 1) then
    write(mesg,'("OBC_registry_init called ",I3, &
      &" times with different registry pointers.")') init_calls
    if (is_root_pe()) call mom_error(warning,"MOM_open_boundary"//mesg)
  endif
 
end subroutine obc_registry_init
 
!> Add file to OBC registry.
function register_file_obc(param_file, CS, OBC_Reg)
  type(param_file_type),    intent(in) :: param_file !< parameter file.
  type(file_obc_cs),        pointer    :: cs         !< file control structure.
  type(obc_registry_type),  pointer    :: obc_reg    !< OBC registry.
  logical                              :: register_file_obc
  character(len=32)  :: casename = "OBC file"        !< This case's name.
 
  if (associated(cs)) then
    call mom_error(warning, "register_file_OBC called with an "// &
                            "associated control structure.")
    return
  endif
  allocate(cs)
 
  ! Register the file for boundary updates.
  call register_obc(casename, param_file, obc_reg)
  register_file_obc = .true.
 
end function register_file_obc
 
!> Clean up the file OBC from registry.
subroutine file_obc_end(CS)
  type(file_obc_cs), pointer    :: cs   !< OBC file control structure.
 
  if (associated(cs)) then
    deallocate(cs)
  endif
end subroutine file_obc_end
 
!> Initialize the segment tracer registry.
subroutine segment_tracer_registry_init(param_file, segment)
  type(param_file_type),      intent(in)      :: param_file !< open file to parse for model parameters
  type(obc_segment_type), intent(inout)       :: segment    !<  the segment
 
  integer, save :: init_calls = 0
 
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "segment_tracer_registry_init" ! This routine's name.
  character(len=256) :: mesg    ! Message for error messages.
 
  if (.not.associated(segment%tr_Reg)) then
    allocate(segment%tr_Reg)
  else
    return
  endif
 
  init_calls = init_calls + 1
 
  ! Read all relevant parameters and write them to the model log.
  if (init_calls == 1) call log_version(param_file, mdl, version, "")
 
! Need to call once per segment with tracers...
! if (init_calls > 1) then
!   write(mesg,'("segment_tracer_registry_init called ",I3, &
!     &" times with different registry pointers.")') init_calls
!   if (is_root_pe()) call MOM_error(WARNING,"MOM_tracer"//mesg)
! endif
 
end subroutine segment_tracer_registry_init
 
subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, &
                                   OBC_scalar, OBC_array)
  type(verticalgrid_type), intent(in)   :: gv         !< ocean vertical grid structure
  type(tracer_type), target             :: tr_ptr     !< A target that can be used to set a pointer to the
                                                      !! stored value of tr. This target must be
                                                      !! an enduring part of the control structure,
                                                      !! because the tracer registry will use this memory,
                                                      !! but it also means that any updates to this
                                                      !! structure in the calling module will be
                                                      !! available subsequently to the tracer registry.
  type(param_file_type), intent(in)     :: param_file !< file to parse for  model parameter values
  type(obc_segment_type), intent(inout) :: segment    !< current segment data structure
  real, optional, intent(in)            :: obc_scalar !< If present, use scalar value for segment tracer
                                                      !! inflow concentration.
  logical, optional, intent(in)         :: obc_array  !< If true, use array values for segment tracer
                                                      !! inflow concentration.
 
 
! Local variables
  integer :: ntseg
  integer :: isd, ied, jsd, jed
  integer :: isdb, iedb, jsdb, jedb
  character(len=256) :: mesg    ! Message for error messages.
 
  call segment_tracer_registry_init(param_file, segment)
 
  if (segment%tr_Reg%ntseg>=max_fields_) then
    write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for &
        &all the tracers being registered via register_tracer.")') segment%tr_Reg%ntseg+1
    call mom_error(fatal,"MOM register_tracer: "//mesg)
  endif
  segment%tr_Reg%ntseg = segment%tr_Reg%ntseg + 1
  ntseg     = segment%tr_Reg%ntseg
 
  isd = segment%HI%isd ; ied = segment%HI%ied
  jsd = segment%HI%jsd ; jed = segment%HI%jed
  isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
  jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
 
  segment%tr_Reg%Tr(ntseg)%Tr => tr_ptr
  segment%tr_Reg%Tr(ntseg)%name = tr_ptr%name
 
  if (segment%tr_Reg%locked) call mom_error(fatal, &
      "MOM register_tracer was called for variable "//trim(segment%tr_Reg%Tr(ntseg)%name)//&
      " with a locked tracer registry.")
 
  if (present(obc_scalar)) segment%tr_Reg%Tr(ntseg)%OBC_inflow_conc = obc_scalar ! initialize tracer value later
  if (present(obc_array)) then
    if (segment%is_E_or_W) then
      allocate(segment%tr_Reg%Tr(ntseg)%t(isdb:iedb,jsd:jed,1:gv%ke));segment%tr_Reg%Tr(ntseg)%t(:,:,:)=0.0
      allocate(segment%tr_Reg%Tr(ntseg)%tres(isdb:iedb,jsd:jed,1:gv%ke));segment%tr_Reg%Tr(ntseg)%tres(:,:,:)=0.0
      segment%tr_Reg%Tr(ntseg)%is_initialized=.false.
    elseif (segment%is_N_or_S) then
      allocate(segment%tr_Reg%Tr(ntseg)%t(isd:ied,jsdb:jedb,1:gv%ke));segment%tr_Reg%Tr(ntseg)%t(:,:,:)=0.0
      allocate(segment%tr_Reg%Tr(ntseg)%tres(isd:ied,jsdb:jedb,1:gv%ke));segment%tr_Reg%Tr(ntseg)%tres(:,:,:)=0.0
      segment%tr_Reg%Tr(ntseg)%is_initialized=.false.
    endif
  endif
 
end subroutine register_segment_tracer
 
!> Clean up the segment tracer registry.
subroutine segment_tracer_registry_end(Reg)
  type(segment_tracer_registry_type), pointer :: reg        !< pointer to tracer registry
 
! Local variables
  integer n
 
  if (associated(reg)) then
    do n = 1, reg%ntseg
      if (associated(reg%Tr(n)%t)) deallocate(reg%Tr(n)%t)
    enddo
    deallocate(reg)
  endif
end subroutine segment_tracer_registry_end
 
subroutine register_temp_salt_segments(GV, OBC, tr_Reg, param_file)
  type(verticalgrid_type),    intent(in)    :: gv         !< ocean vertical grid structure
  type(ocean_obc_type),       pointer       :: obc        !< Open boundary structure
  type(tracer_registry_type), pointer       :: tr_reg     !< Tracer registry
  type(param_file_type),      intent(in)    :: param_file !< file to parse for  model parameter values
 
! Local variables
  integer :: isd, ied, isdb, iedb, jsd, jed, jsdb, jedb, nz, nf
  integer :: i, j, k, n
  character(len=32)  :: name
  type(obc_segment_type), pointer :: segment => null() ! pointer to segment type list
  type(tracer_type), pointer      :: tr_ptr => null()
 
  if (.not. associated(obc)) return
 
  do n=1, obc%number_of_segments
    segment=>obc%segment(n)
    if (.not. segment%on_pe) cycle
 
    if (associated(segment%tr_Reg)) &
         call mom_error(fatal,"register_temp_salt_segments: tracer array was previously allocated")
 
    name = 'temp'
    call tracer_name_lookup(tr_reg, tr_ptr, name)
    call register_segment_tracer(tr_ptr, param_file, gv, segment, &
                                 obc_array=segment%temp_segment_data_exists)
    name = 'salt'
    call tracer_name_lookup(tr_reg, tr_ptr, name)
    call register_segment_tracer(tr_ptr, param_file, gv, segment, &
                                 obc_array=segment%salt_segment_data_exists)
  enddo
 
end subroutine register_temp_salt_segments
 
subroutine fill_temp_salt_segments(G, OBC, tv)
  type(ocean_grid_type),      intent(inout) :: g          !< Ocean grid structure
  type(ocean_obc_type),       pointer       :: obc        !< Open boundary structure
  type(thermo_var_ptrs),      intent(inout) :: tv         !< Thermodynamics structure
 
! Local variables
  integer :: isd, ied, isdb, iedb, jsd, jed, jsdb, jedb, n, nz
  integer :: i, j, k
  type(obc_segment_type), pointer :: segment => null() ! pointer to segment type list
 
  if (.not. associated(obc)) return
  if (.not. associated(tv%T) .and. associated(tv%S)) return
  ! Both temperature and salinity fields
 
  call pass_var(tv%T, g%Domain)
  call pass_var(tv%S, g%Domain)
 
  nz = g%ke
 
  do n=1, obc%number_of_segments
    segment => obc%segment(n)
    if (.not. segment%on_pe) cycle
 
    isd = segment%HI%isd ; ied = segment%HI%ied
    jsd = segment%HI%jsd ; jed = segment%HI%jed
    isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
    jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
 
    ! Fill with T and S values
    if (segment%is_E_or_W) then
      i=segment%HI%IsdB
      do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
        if (segment%direction == obc_direction_w) then
          segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i+1,j,k)
          segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i+1,j,k)
        else
          segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j,k)
          segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j,k)
        endif
      enddo ; enddo
    else
      j=segment%HI%JsdB
      do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
        if (segment%direction == obc_direction_s) then
          segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j+1,k)
          segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j+1,k)
        else
          segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j,k)
          segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j,k)
        endif
      enddo ; enddo
    endif
    segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:)
    segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:)
  enddo
  call setup_obc_tracer_reservoirs(g, obc)
end subroutine fill_temp_salt_segments
 
!> Find the region outside of all open boundary segments and
!! make sure it is set to land mask. Gonna need to know global land
!! mask as well to get it right...
subroutine mask_outside_obcs(G, US, param_file, OBC)
  type(dyn_horgrid_type),       intent(inout) :: G          !< Ocean grid structure
  type(param_file_type),        intent(in)    :: param_file !< Parameter file handle
  type(ocean_obc_type),         pointer       :: OBC        !< Open boundary structure
  type(unit_scale_type),        intent(in)    :: US         !< A dimensional unit scaling type
 
  ! Local variables
  integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n
  integer :: i, j
  logical :: fatal_error = .false.
  real    :: min_depth
  integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2
  character(len=256) :: mesg    ! Message for error messages.
  type(obc_segment_type), pointer :: segment => null() ! pointer to segment type list
  real, allocatable, dimension(:,:) :: color, color2  ! For sorting inside from outside,
                                                      ! two different ways
 
  if (.not. associated(obc)) return
 
  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
                 units="m", default=0.0, scale=us%m_to_Z, do_not_log=.true.)
 
  allocate(color(g%isd:g%ied, g%jsd:g%jed)) ; color = 0
  allocate(color2(g%isd:g%ied, g%jsd:g%jed)) ; color2 = 0
 
 
  ! Paint a frame around the outside.
  do j=g%jsd,g%jed
    color(g%isd,j) = cedge
    color(g%ied,j) = cedge
    color2(g%isd,j) = cedge
    color2(g%ied,j) = cedge
  enddo
  do i=g%isd,g%ied
    color(i,g%jsd) = cedge
    color(i,g%jed) = cedge
    color2(i,g%jsd) = cedge
    color2(i,g%jed) = cedge
  enddo
 
  ! Set color to cland in the land. Note that this is before the land
  ! mask has been initialized, set mask values based on depth.
  do j=g%jsd,g%jed
    do i=g%isd,g%ied
      if (g%bathyT(i,j) <= min_depth) color(i,j) = cland
      if (g%bathyT(i,j) <= min_depth) color2(i,j) = cland
    enddo
  enddo
 
  do j=g%jsd,g%jed ; do i=g%IsdB+1,g%IedB-1
    if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
      if (color(i,j) == 0.0) color(i,j) = cout
      if (color(i+1,j) == 0.0) color(i+1,j) = cin
    elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
      if (color(i,j) == 0.0) color(i,j) = cin
      if (color(i+1,j) == 0.0) color(i+1,j) = cout
    endif
  enddo ; enddo
  do j=g%JsdB+1,g%JedB-1 ; do i=g%isd,g%ied
    if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
      if (color(i,j) == 0.0) color(i,j) = cout
      if (color(i,j+1) == 0.0) color(i,j+1) = cin
    elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
      if (color(i,j) == 0.0) color(i,j) = cin
      if (color(i,j+1) == 0.0) color(i,j+1) = cout
    endif
  enddo ; enddo
 
  do j=g%JsdB+1,g%JedB-1 ; do i=g%isd,g%ied
    if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
      if (color2(i,j) == 0.0) color2(i,j) = cout
      if (color2(i,j+1) == 0.0) color2(i,j+1) = cin
    elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
      if (color2(i,j) == 0.0) color2(i,j) = cin
      if (color2(i,j+1) == 0.0) color2(i,j+1) = cout
    endif
  enddo ; enddo
  do j=g%jsd,g%jed ; do i=g%IsdB+1,g%IedB-1
    if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
      if (color2(i,j) == 0.0) color2(i,j) = cout
      if (color2(i+1,j) == 0.0) color2(i+1,j) = cin
    elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
      if (color2(i,j) == 0.0) color2(i,j) = cin
      if (color2(i+1,j) == 0.0) color2(i+1,j) = cout
    endif
  enddo ; enddo
 
  ! Do the flood fill until there are no more uncolored cells.
  call flood_fill(g, color, cin, cout, cland)
  call flood_fill2(g, color2, cin, cout, cland)
 
  ! Use the color to set outside to min_depth on this process.
  do j=g%jsd,g%jed ; do i=g%isd,g%ied
    if (color(i,j) /= color2(i,j)) then
      fatal_error = .true.
      write(mesg,'("MOM_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", &
          "the masking of the outside grid points.")') i, j
      call mom_error(warning,"MOM register_tracer: "//mesg, all_print=.true.)
    endif
    if (color(i,j) == cout) g%bathyT(i,j) = min_depth
  enddo ; enddo
  if (fatal_error) call mom_error(fatal, &
      "MOM_open_boundary: inconsistent OBC segments.")
 
  deallocate(color)
  deallocate(color2)
end subroutine mask_outside_obcs
 
!> flood the cin, cout values
subroutine flood_fill(G, color, cin, cout, cland)
  type(dyn_horgrid_type), intent(inout) :: G      !< Ocean grid structure
  real, dimension(:,:),   intent(inout) :: color  !< For sorting inside from outside
  integer, intent(in) :: cin    !< color for inside the domain
  integer, intent(in) :: cout   !< color for outside the domain
  integer, intent(in) :: cland  !< color for inside the land mask
 
! Local variables
  integer :: i, j, ncount
 
  ncount = 1
  do while (ncount > 0)
    ncount = 0
    do j=g%jsd+1,g%jed-1
      do i=g%isd+1,g%ied-1
        if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0) then
          color(i,j) = color(i-1,j)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0) then
          color(i,j) = color(i+1,j)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0) then
          color(i,j) = color(i,j-1)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0) then
          color(i,j) = color(i,j+1)
          ncount = ncount + 1
        endif
      enddo
    enddo
    do j=g%jed-1,g%jsd+1,-1
      do i=g%ied-1,g%isd+1,-1
        if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0) then
          color(i,j) = color(i-1,j)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0) then
          color(i,j) = color(i+1,j)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0) then
          color(i,j) = color(i,j-1)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0) then
          color(i,j) = color(i,j+1)
          ncount = ncount + 1
        endif
      enddo
    enddo
    call pass_var(color, g%Domain)
    call sum_across_pes(ncount)
  enddo
 
end subroutine flood_fill
 
!> flood the cin, cout values
subroutine flood_fill2(G, color, cin, cout, cland)
  type(dyn_horgrid_type), intent(inout) :: G       !< Ocean grid structure
  real, dimension(:,:),   intent(inout) :: color   !< For sorting inside from outside
  integer, intent(in) :: cin    !< color for inside the domain
  integer, intent(in) :: cout   !< color for outside the domain
  integer, intent(in) :: cland  !< color for inside the land mask
 
! Local variables
  integer :: i, j, ncount
 
  ncount = 1
  do while (ncount > 0)
    ncount = 0
    do i=g%isd+1,g%ied-1
      do j=g%jsd+1,g%jed-1
        if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0) then
          color(i,j) = color(i-1,j)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0) then
          color(i,j) = color(i+1,j)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0) then
          color(i,j) = color(i,j-1)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0) then
          color(i,j) = color(i,j+1)
          ncount = ncount + 1
        endif
      enddo
    enddo
    do i=g%ied-1,g%isd+1,-1
      do j=g%jed-1,g%jsd+1,-1
        if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0) then
          color(i,j) = color(i-1,j)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0) then
          color(i,j) = color(i+1,j)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0) then
          color(i,j) = color(i,j-1)
          ncount = ncount + 1
        endif
        if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0) then
          color(i,j) = color(i,j+1)
          ncount = ncount + 1
        endif
      enddo
    enddo
    call pass_var(color, g%Domain)
    call sum_across_pes(ncount)
  enddo
 
end subroutine flood_fill2
 
!> Register OBC segment data for restarts
subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart_CSp, &
                                           use_temperature)
  type(hor_index_type),    intent(in) :: hi !< Horizontal indices
  type(verticalgrid_type), pointer    :: gv !< Container for vertical grid information
  type(ocean_obc_type),    pointer    :: obc !< OBC data structure, data intent(inout)
  type(tracer_registry_type), pointer :: reg !< pointer to tracer registry
  type(param_file_type),   intent(in) :: param_file !< Parameter file handle
  type(mom_restart_cs),    pointer    :: restart_csp !< Restart structure, data intent(inout)
  logical,                 intent(in) :: use_temperature !< If true, T and S are used
  ! Local variables
  type(vardesc) :: vd
  integer       :: m, n
  character(len=100) :: mesg
  type(obc_segment_type), pointer :: segment=>null()
 
  if (.not. associated(obc)) &
       call mom_error(fatal, "open_boundary_register_restarts: Called with "//&
                      "uninitialized OBC control structure")
 
  if (associated(obc%rx_normal) .or. associated(obc%ry_normal) .or. &
      associated(obc%rx_oblique) .or. associated(obc%ry_oblique) .or. associated(obc%cff_normal)) &
       call mom_error(fatal, "open_boundary_register_restarts: Restart "//&
                      "arrays were previously allocated")
 
  if (associated(obc%tres_x) .or. associated(obc%tres_y)) &
       call mom_error(fatal, "open_boundary_register_restarts: Restart "//&
                      "arrays were previously allocated")
 
  ! *** This is a temporary work around for restarts with OBC segments.
  ! This implementation uses 3D arrays solely for restarts. We need
  ! to be able to add 2D ( x,z or y,z ) data to restarts to avoid using
  ! so much memory and disk space. ***
  if (obc%radiation_BCs_exist_globally) then
    allocate(obc%rx_normal(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke))
    obc%rx_normal(:,:,:) = 0.0
    vd = var_desc("rx_normal", "m s-1", "Normal Phase Speed for EW radiation OBCs", 'u', 'L')
    call register_restart_field(obc%rx_normal, vd, .false., restart_csp)
    allocate(obc%ry_normal(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke))
    obc%ry_normal(:,:,:) = 0.0
    vd = var_desc("ry_normal", "m s-1", "Normal Phase Speed for NS radiation OBCs", 'v', 'L')
    call register_restart_field(obc%ry_normal, vd, .false., restart_csp)
  endif
  if (obc%oblique_BCs_exist_globally) then
    allocate(obc%rx_oblique(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke))
    obc%rx_oblique(:,:,:) = 0.0
    vd = var_desc("rx_oblique", "m2 s-2", "Radiation Speed Squared for EW oblique OBCs", 'u', 'L')
    call register_restart_field(obc%rx_oblique, vd, .false., restart_csp)
    allocate(obc%ry_oblique(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke))
    obc%ry_oblique(:,:,:) = 0.0
    vd = var_desc("ry_oblique", "m2 s-2", "Radiation Speed Squared for NS oblique OBCs", 'v', 'L')
    call register_restart_field(obc%ry_oblique, vd, .false., restart_csp)
    allocate(obc%cff_normal(hi%IsdB:hi%IedB,hi%jsdB:hi%jedB,gv%ke))
    obc%cff_normal(:,:,:) = 0.0
    vd = var_desc("cff_normal", "m2 s-2", "denominator for oblique OBCs", 'q', 'L')
    call register_restart_field(obc%cff_normal, vd, .false., restart_csp)
  endif
 
  if (reg%ntr == 0) return
  if (.not. associated(obc%tracer_x_reservoirs_used)) then
    obc%ntr = reg%ntr
    allocate(obc%tracer_x_reservoirs_used(reg%ntr))
    allocate(obc%tracer_y_reservoirs_used(reg%ntr))
    obc%tracer_x_reservoirs_used(:) = .false.
    obc%tracer_y_reservoirs_used(:) = .false.
    call parse_for_tracer_reservoirs(obc, param_file, use_temperature)
  else
    ! This would be coming from user code such as DOME.
    if (obc%ntr /= reg%ntr) then
!        call MOM_error(FATAL, "open_boundary_regiser_restarts: Inconsistent value for ntr")
      write(mesg,'("Inconsistent values for ntr ", I8," and ",I8,".")') obc%ntr, reg%ntr
      call mom_error(warning, 'open_boundary_register_restarts: '//mesg)
    endif
  endif
 
  ! Still painfully inefficient, now in four dimensions.
  if (any(obc%tracer_x_reservoirs_used)) then
    allocate(obc%tres_x(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke,obc%ntr))
    obc%tres_x(:,:,:,:) = 0.0
    do m=1,obc%ntr
      if (obc%tracer_x_reservoirs_used(m)) then
        write(mesg,'("tres_x_",I3.3)') m
        vd = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L')
        call register_restart_field(obc%tres_x(:,:,:,m), vd, .false., restart_csp)
      endif
    enddo
  endif
  if (any(obc%tracer_y_reservoirs_used)) then
    allocate(obc%tres_y(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke,obc%ntr))
    obc%tres_y(:,:,:,:) = 0.0
    do m=1,obc%ntr
      if (obc%tracer_y_reservoirs_used(m)) then
        write(mesg,'("tres_y_",I3.3)') m
        vd = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L')
        call register_restart_field(obc%tres_y(:,:,:,m), vd, .false., restart_csp)
      endif
    enddo
  endif
 
end subroutine open_boundary_register_restarts
 
!> Update the OBC tracer reservoirs after the tracers have been updated.
subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
  type(ocean_grid_type),                     intent(in) :: g   !< The ocean's grid structure
  type(verticalgrid_type),                   intent(in) :: gv  !<  Ocean vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uhr !< accumulated volume/mass flux through
                                                               !! the zonal face [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vhr !< accumulated volume/mass flux through
                                                               !! the meridional face [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h   !< layer thickness after advection
                                                               !! [H ~> m or kg m-2]
  type(ocean_obc_type),                      pointer    :: obc !< Open boundary structure
  real,                                      intent(in) :: dt  !< time increment [T ~> s]
  type(tracer_registry_type),                pointer    :: reg !< pointer to tracer registry
  ! Local variables
  type(obc_segment_type), pointer :: segment=>null()
  real :: u_l_in, u_l_out ! The zonal distance moved in or out of a cell [L ~> m]
  real :: v_l_in, v_l_out ! The meridional distance moved in or out of a cell [L ~> m]
  real :: fac1            ! The denominator of the expression for tracer updates [nondim]
  integer :: i, j, k, m, n, ntr, nz
  integer :: ishift, idir, jshift, jdir
 
  nz = gv%ke
  ntr = reg%ntr
  if (associated(obc)) then ; if (obc%OBC_pe) then ; do n=1,obc%number_of_segments
    segment=>obc%segment(n)
    if (.not. associated(segment%tr_Reg)) cycle
    if (segment%is_E_or_W) then
      do j=segment%HI%jsd,segment%HI%jed
        i = segment%HI%IsdB
        ! ishift+I corresponds to the nearest interior tracer cell index
        ! idir switches the sign of the flow so that positive is into the reservoir
        if (segment%direction == obc_direction_w) then
          ishift = 1 ; idir = -1
        else
          ishift = 0 ; idir = 1
        endif
        ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
        do m=1,ntr ; if (associated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
          u_l_out = max(0.0, (idir*uhr(i,j,k))*segment%Tr_InvLscale_out / (h(i+ishift,j,k)*g%dyCu(i,j)))
          u_l_in  = min(0.0, (idir*uhr(i,j,k))*segment%Tr_InvLscale_in  / (h(i+ishift,j,k)*g%dyCu(i,j)))
          fac1 = 1.0 + (u_l_out-u_l_in)
          segment%tr_Reg%Tr(m)%tres(i,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,j,k) + &
                            (u_l_out*reg%Tr(m)%t(i+ishift,j,k) - &
                             u_l_in*segment%tr_Reg%Tr(m)%t(i,j,k)))
          if (associated(obc%tres_x)) obc%tres_x(i,j,k,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
        enddo ; endif ; enddo
      enddo
    else
      do i=segment%HI%isd,segment%HI%ied
        j = segment%HI%JsdB
        ! jshift+J corresponds to the nearest interior tracer cell index
        ! jdir switches the sign of the flow so that positive is into the reservoir
        if (segment%direction == obc_direction_s) then
          jshift = 1 ; jdir = -1
        else
          jshift = 0 ; jdir = 1
        endif
        ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
        do m=1,ntr ; if (associated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
          v_l_out = max(0.0, (jdir*vhr(i,j,k))*segment%Tr_InvLscale_out / (h(i,j+jshift,k)*g%dxCv(i,j)))
          v_l_in  = min(0.0, (jdir*vhr(i,j,k))*segment%Tr_InvLscale_in  / (h(i,j+jshift,k)*g%dxCv(i,j)))
          fac1 = 1.0 + (v_l_out-v_l_in)
          segment%tr_Reg%Tr(m)%tres(i,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,j,k) + &
                            (v_l_out*reg%Tr(m)%t(i,j+jshift,k) - &
                             v_l_in*segment%tr_Reg%Tr(m)%t(i,j,k)))
          if (associated(obc%tres_y)) obc%tres_y(i,j,k,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
        enddo ; endif ; enddo
      enddo
    endif
  enddo ; endif ; endif
 
end subroutine update_segment_tracer_reservoirs
 
!> Adjust interface heights to fit the bathymetry and diagnose layer thickness.
!!
!! If the bottom most interface is below the topography then the bottom-most
!! layers are contracted to GV%Angstrom_m.
!! If the bottom most interface is above the topography then the entire column
!! is dilated (expanded) to fill the void.
!!   @remark{There is a (hard-wired) "tolerance" parameter such that the
!! criteria for adjustment must equal or exceed 10cm.}
subroutine adjustsegmentetatofitbathymetry(G, GV, US, segment,fld)
  type(ocean_grid_type),                      intent(in)    :: G   !< The ocean's grid structure
  type(verticalgrid_type),                    intent(in)    :: GV  !< The ocean's vertical grid structure
  type(unit_scale_type),                      intent(in)    :: US  !< A dimensional unit scaling type
  type(obc_segment_type),                     intent(inout) :: segment !< pointer to segment type
  integer,                                    intent(in)    :: fld  !< field index to adjust thickness
  ! Local variables
  integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
  integer :: n
  real, allocatable, dimension(:,:,:) :: eta ! Segment source data interface heights, [Z -> m]
  real :: hTolerance = 0.1 !<  Tolerance to exceed adjustment criteria [Z ~> m]
  real :: hTmp, eTmp, dilate
  character(len=100) :: mesg
 
  htolerance = 0.1*us%m_to_Z
 
  nz = size(segment%field(fld)%dz_src,3)
 
  if (segment%is_E_or_W) then
     ! segment thicknesses are defined at cell face centers.
     is = segment%HI%isdB ; ie = segment%HI%iedB
     js = segment%HI%jsd ; je = segment%HI%jed
  else
     is = segment%HI%isd ; ie = segment%HI%ied
     js = segment%HI%jsdB ; je = segment%HI%jedB
  endif
  allocate(eta(is:ie,js:je,nz+1))
  contractions=0; dilations=0
  do j=js,je ; do i=is,ie
    eta(i,j,1)=0.0   ! segment data are assumed to be located on a static grid
    ! For remapping calls, the entire column will be dilated
    ! by a factor equal to the ratio of the sum of the geopotential referenced
    ! source data thicknesses, and the current model thicknesses. This could be
    ! an issue to be addressed, for instance if we are placing open boundaries
    ! under ice shelf cavities.
    do k=2,nz+1
      eta(i,j,k)=eta(i,j,k-1)-segment%field(fld)%dz_src(i,j,k-1)
    enddo
    ! The normal slope at the boundary is zero by a
    ! previous call to open_boundary_impose_normal_slope
    do k=nz+1,1,-1
      if (-eta(i,j,k) > segment%Htot(i,j) + htolerance) then
         eta(i,j,k) = -segment%Htot(i,j)
         contractions = contractions + 1
      endif
    enddo
 
    do k=1,nz
      ! Collapse layers to thinnest possible if the thickness less than
      ! the thinnest possible (or negative).
      if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) then
         eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
         segment%field(fld)%dz_src(i,j,k) = gv%Angstrom_Z
      else
         segment%field(fld)%dz_src(i,j,k) = (eta(i,j,k) - eta(i,j,k+1))
      endif
    enddo
 
    !   The whole column is dilated to accommodate deeper topography than
    ! the bathymetry would indicate.
    if (-eta(i,j,nz+1) < segment%Htot(i,j) - htolerance) then
       dilations = dilations + 1
       ! expand bottom-most cell only
       eta(i,j,nz+1) = -segment%Htot(i,j)
       segment%field(fld)%dz_src(i,j,nz)= eta(i,j,nz)-eta(i,j,nz+1)
       ! if (eta(i,j,1) <= eta(i,j,nz+1)) then
       !   do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = (eta(i,j,1) + G%bathyT(i,j)) / real(nz) ; enddo
       ! else
       !   dilate = (eta(i,j,1) + G%bathyT(i,j)) / (eta(i,j,1) - eta(i,j,nz+1))
       !   do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k) * dilate ; enddo
       ! endif
       !do k=nz,2,-1 ; eta(i,j,K) = eta(i,j,K+1) + segment%field(fld)%dz_src(i,j,k) ; enddo
    endif
    ! Now convert thicknesses to units of H.
    do k=1,nz
      segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k)*gv%Z_to_H
    enddo
  enddo; enddo
 
  ! can not do communication call here since only PEs on the current segment are here
 
  ! call sum_across_PEs(contractions)
  ! if ((contractions > 0) .and. (is_root_pe())) then
  !    write(mesg,'("Thickness OBCs were contracted ",'// &
  !         '"to fit topography in ",I8," places.")') contractions
  !    call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
  ! endif
  ! call sum_across_PEs(dilations)
  ! if ((dilations > 0) .and. (is_root_pe())) then
  !    write(mesg,'("Thickness OBCs were dilated ",'// &
  !         '"to fit topography in ",I8," places.")') dilations
  !    call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
  ! endif
  deallocate(eta)
 
 
 
end subroutine adjustsegmentetatofitbathymetry
 
!> \namespace mom_open_boundary
!! This module implements some aspects of internal open boundary
!! conditions in MOM.
!!
!! A small fragment of the grid is shown below:
!!
!!    j+1  x ^ x ^ x   At x:  q, CoriolisBu
!!    j+1  > o > o >   At ^:  v, tauy
!!    j    x ^ x ^ x   At >:  u, taux
!!    j    > o > o >   At o:  h, bathyT, buoy, tr, T, S, Rml, ustar
!!    j-1  x ^ x ^ x
!!        i-1  i  i+1  At x & ^:
!!           i  i+1    At > & o:
!!
!! The boundaries always run through q grid points (x).
 
end module mom_open_boundary