MOM6
mom_open_boundary Module Reference

Detailed Description

Controls where open boundary conditions are applied.

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

Data Types

type  file_obc_cs
 Control structure for open boundaries that read from files. Probably lots to update here. More...
 
type  obc_registry_type
 Type to carry basic OBC information needed for updating values. More...
 
type  obc_segment_data_type
 Open boundary segment data from files (mostly). More...
 
type  obc_segment_tracer_type
 Tracer on OBC segment data structure, for putting into a segment tracer registry. More...
 
type  obc_segment_type
 Open boundary segment data structure. More...
 
type  obc_struct_type
 Type to carry something (what] for the OBC registry. More...
 
type  ocean_obc_type
 Open-boundary data. More...
 
type  segment_tracer_registry_type
 Registry type for tracers on segments. More...
 

Functions/Subroutines

subroutine, public open_boundary_config (G, US, param_file, OBC)
 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. More...
 
subroutine initialize_segment_data (G, OBC, PF)
 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. More...
 
subroutine setup_segment_indices (G, seg, Is_obc, Ie_obc, Js_obc, Je_obc)
 Define indices for segment and store in hor_index_type using global segment bounds corresponding to q-points. More...
 
subroutine setup_u_point_obc (OBC, G, US, segment_str, l_seg, PF, reentrant_y)
 Parse an OBC_SEGMENT_%%% string starting with "I=" and configure placement and type of OBC accordingly. More...
 
subroutine setup_v_point_obc (OBC, G, US, segment_str, l_seg, PF, reentrant_x)
 Parse an OBC_SEGMENT_%%% string starting with "J=" and configure placement and type of OBC accordingly. More...
 
subroutine parse_segment_str (ni_global, nj_global, segment_str, l, m, n, action_str, reentrant)
 Parse an OBC_SEGMENT_%%% string. More...
 
subroutine parse_segment_data_str (segment_str, var, value, filenam, fieldnam, fields, num_fields, debug)
 Parse an OBC_SEGMENT_%%_DATA string. More...
 
subroutine parse_for_tracer_reservoirs (OBC, PF, use_temperature)
 Parse all the OBC_SEGMENT_%%_DATA strings again to see which need tracer reservoirs (all pes need to know). More...
 
subroutine parse_segment_param_real (segment_str, var, param_value, debug)
 Parse an OBC_SEGMENT_%%_PARAMS string. More...
 
subroutine, public open_boundary_init (G, GV, US, param_file, OBC, restart_CSp)
 Initialize open boundary control structure and do any necessary rescaling of OBC fields that have been read from a restart file. More...
 
logical function, public open_boundary_query (OBC, apply_open_OBC, apply_specified_OBC, apply_Flather_OBC, apply_nudged_OBC, needs_ext_seg_data)
 
subroutine open_boundary_dealloc (OBC)
 Deallocate open boundary data. More...
 
subroutine, public open_boundary_end (OBC)
 Close open boundary data. More...
 
subroutine, public open_boundary_impose_normal_slope (OBC, G, depth)
 Sets the slope of bathymetry normal to an open bounndary to zero. More...
 
subroutine, public open_boundary_impose_land_mask (OBC, G, areaCu, areaCv, US)
 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. More...
 
subroutine setup_obc_tracer_reservoirs (G, OBC)
 Make sure the OBC tracer reservoirs are initialized. More...
 
subroutine, public radiation_open_bdry_conds (OBC, u_new, u_old, v_new, v_old, G, US, dt)
 Apply radiation conditions to 3D u,v at open boundaries. More...
 
subroutine, public open_boundary_apply_normal_flow (OBC, G, u, v)
 Applies OBC values stored in segments to 3d u,v fields. More...
 
subroutine, public open_boundary_zero_normal_flow (OBC, G, u, v)
 Applies zero values to 3d u,v fields on OBC segments. More...
 
subroutine gradient_at_q_points (G, segment, uvel, vvel)
 Calculate the tangential gradient of the normal flow at the boundary q-points. More...
 
subroutine, public set_tracer_data (OBC, tv, h, G, PF, tracer_Reg)
 Sets the initial values of the tracer open boundary conditions. Redoing this elsewhere. More...
 
integer function lookup_seg_field (OBC_seg, field)
 Needs documentation. More...
 
subroutine allocate_obc_segment_data (OBC, segment)
 Allocate segment data fields. More...
 
subroutine deallocate_obc_segment_data (OBC, segment)
 Deallocate segment data fields. More...
 
subroutine, public open_boundary_test_extern_uv (G, OBC, u, v)
 Set tangential velocities outside of open boundaries to silly values (used for checking the interior state is independent of values outside of the domain). More...
 
subroutine, public open_boundary_test_extern_h (G, GV, OBC, h)
 Set thicknesses outside of open boundaries to silly values (used for checking the interior state is independent of values outside of the domain). More...
 
subroutine, public update_obc_segment_data (G, GV, US, OBC, tv, h, Time)
 Update the OBC values on the segments. More...
 
subroutine, public register_obc (name, param_file, Reg)
 register open boundary objects for boundary updates. More...
 
subroutine, public obc_registry_init (param_file, Reg)
 This routine include declares and sets the variable "version". More...
 
logical function, public register_file_obc (param_file, CS, OBC_Reg)
 Add file to OBC registry. More...
 
subroutine, public file_obc_end (CS)
 Clean up the file OBC from registry. More...
 
subroutine, public segment_tracer_registry_init (param_file, segment)
 Initialize the segment tracer registry. More...
 
subroutine, public register_segment_tracer (tr_ptr, param_file, GV, segment, OBC_scalar, OBC_array)
 
subroutine, public segment_tracer_registry_end (Reg)
 Clean up the segment tracer registry. More...
 
subroutine, public register_temp_salt_segments (GV, OBC, tr_Reg, param_file)
 
subroutine, public fill_temp_salt_segments (G, OBC, tv)
 
subroutine mask_outside_obcs (G, US, param_file, OBC)
 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... More...
 
subroutine flood_fill (G, color, cin, cout, cland)
 flood the cin, cout values More...
 
subroutine flood_fill2 (G, color, cin, cout, cland)
 flood the cin, cout values More...
 
subroutine, public open_boundary_register_restarts (HI, GV, OBC, Reg, param_file, restart_CSp, use_temperature)
 Register OBC segment data for restarts. More...
 
subroutine, public update_segment_tracer_reservoirs (G, GV, uhr, vhr, h, OBC, dt, Reg)
 Update the OBC tracer reservoirs after the tracers have been updated. More...
 
subroutine adjustsegmentetatofitbathymetry (G, GV, US, segment, fld)
 Adjust interface heights to fit the bathymetry and diagnose layer thickness. More...
 

Variables

integer, parameter, public obc_none = 0
 Indicates the use of no open boundary. More...
 
integer, parameter, public obc_simple = 1
 Indicates the use of a simple inflow open boundary. More...
 
integer, parameter, public obc_wall = 2
 Indicates the use of a closed wall. More...
 
integer, parameter, public obc_flather = 3
 Indicates the use of a Flather open boundary. More...
 
integer, parameter, public obc_radiation = 4
 Indicates the use of a radiation open boundary. More...
 
integer, parameter, public obc_direction_n = 100
 Indicates the boundary is an effective northern boundary. More...
 
integer, parameter, public obc_direction_s = 200
 Indicates the boundary is an effective southern boundary. More...
 
integer, parameter, public obc_direction_e = 300
 Indicates the boundary is an effective eastern boundary. More...
 
integer, parameter, public obc_direction_w = 400
 Indicates the boundary is an effective western boundary. More...
 
integer, parameter max_obc_fields = 100
 Maximum number of data fields needed for OBC segments. More...
 
integer id_clock_pass
 A CPU time clock. More...
 
character(len=40) mdl = "MOM_open_boundary"
 This module's name. More...
 

Function/Subroutine Documentation

◆ adjustsegmentetatofitbathymetry()

subroutine mom_open_boundary::adjustsegmentetatofitbathymetry ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(obc_segment_type), intent(inout)  segment,
integer, intent(in)  fld 
)
private

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 GVAngstrom_m. If the bottom most interface is above the topography then the entire column is dilated (expanded) to fill the void.

Remarks
{There is a (hard-wired) "tolerance" parameter such that the criteria for adjustment must equal or exceed 10cm.}
Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in,out]segmentpointer to segment type
[in]fldfield index to adjust thickness

Definition at line 4593 of file MOM_open_boundary.F90.

4593  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
4594  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
4595  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
4596  type(OBC_segment_type), intent(inout) :: segment !< pointer to segment type
4597  integer, intent(in) :: fld !< field index to adjust thickness
4598  ! Local variables
4599  integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
4600  integer :: n
4601  real, allocatable, dimension(:,:,:) :: eta ! Segment source data interface heights, [Z -> m]
4602  real :: hTolerance = 0.1 !< Tolerance to exceed adjustment criteria [Z ~> m]
4603  real :: hTmp, eTmp, dilate
4604  character(len=100) :: mesg
4605 
4606  htolerance = 0.1*us%m_to_Z
4607 
4608  nz = size(segment%field(fld)%dz_src,3)
4609 
4610  if (segment%is_E_or_W) then
4611  ! segment thicknesses are defined at cell face centers.
4612  is = segment%HI%isdB ; ie = segment%HI%iedB
4613  js = segment%HI%jsd ; je = segment%HI%jed
4614  else
4615  is = segment%HI%isd ; ie = segment%HI%ied
4616  js = segment%HI%jsdB ; je = segment%HI%jedB
4617  endif
4618  allocate(eta(is:ie,js:je,nz+1))
4619  contractions=0; dilations=0
4620  do j=js,je ; do i=is,ie
4621  eta(i,j,1)=0.0 ! segment data are assumed to be located on a static grid
4622  ! For remapping calls, the entire column will be dilated
4623  ! by a factor equal to the ratio of the sum of the geopotential referenced
4624  ! source data thicknesses, and the current model thicknesses. This could be
4625  ! an issue to be addressed, for instance if we are placing open boundaries
4626  ! under ice shelf cavities.
4627  do k=2,nz+1
4628  eta(i,j,k)=eta(i,j,k-1)-segment%field(fld)%dz_src(i,j,k-1)
4629  enddo
4630  ! The normal slope at the boundary is zero by a
4631  ! previous call to open_boundary_impose_normal_slope
4632  do k=nz+1,1,-1
4633  if (-eta(i,j,k) > segment%Htot(i,j) + htolerance) then
4634  eta(i,j,k) = -segment%Htot(i,j)
4635  contractions = contractions + 1
4636  endif
4637  enddo
4638 
4639  do k=1,nz
4640  ! Collapse layers to thinnest possible if the thickness less than
4641  ! the thinnest possible (or negative).
4642  if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) then
4643  eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
4644  segment%field(fld)%dz_src(i,j,k) = gv%Angstrom_Z
4645  else
4646  segment%field(fld)%dz_src(i,j,k) = (eta(i,j,k) - eta(i,j,k+1))
4647  endif
4648  enddo
4649 
4650  ! The whole column is dilated to accommodate deeper topography than
4651  ! the bathymetry would indicate.
4652  if (-eta(i,j,nz+1) < segment%Htot(i,j) - htolerance) then
4653  dilations = dilations + 1
4654  ! expand bottom-most cell only
4655  eta(i,j,nz+1) = -segment%Htot(i,j)
4656  segment%field(fld)%dz_src(i,j,nz)= eta(i,j,nz)-eta(i,j,nz+1)
4657  ! if (eta(i,j,1) <= eta(i,j,nz+1)) then
4658  ! do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = (eta(i,j,1) + G%bathyT(i,j)) / real(nz) ; enddo
4659  ! else
4660  ! dilate = (eta(i,j,1) + G%bathyT(i,j)) / (eta(i,j,1) - eta(i,j,nz+1))
4661  ! do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k) * dilate ; enddo
4662  ! endif
4663  !do k=nz,2,-1 ; eta(i,j,K) = eta(i,j,K+1) + segment%field(fld)%dz_src(i,j,k) ; enddo
4664  endif
4665  ! Now convert thicknesses to units of H.
4666  do k=1,nz
4667  segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k)*gv%Z_to_H
4668  enddo
4669  enddo; enddo
4670 
4671  ! can not do communication call here since only PEs on the current segment are here
4672 
4673  ! call sum_across_PEs(contractions)
4674  ! if ((contractions > 0) .and. (is_root_pe())) then
4675  ! write(mesg,'("Thickness OBCs were contracted ",'// &
4676  ! '"to fit topography in ",I8," places.")') contractions
4677  ! call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
4678  ! endif
4679  ! call sum_across_PEs(dilations)
4680  ! if ((dilations > 0) .and. (is_root_pe())) then
4681  ! write(mesg,'("Thickness OBCs were dilated ",'// &
4682  ! '"to fit topography in ",I8," places.")') dilations
4683  ! call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
4684  ! endif
4685  deallocate(eta)
4686 
4687 
4688 

Referenced by update_obc_segment_data().

Here is the caller graph for this function:

◆ allocate_obc_segment_data()

subroutine mom_open_boundary::allocate_obc_segment_data ( type(ocean_obc_type), pointer  OBC,
type(obc_segment_type), intent(inout)  segment 
)
private

Allocate segment data fields.

Parameters
obcOpen boundary structure
[in,out]segmentOpen boundary segment

Definition at line 3171 of file MOM_open_boundary.F90.

3171  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
3172  type(OBC_segment_type), intent(inout) :: segment !< Open boundary segment
3173  ! Local variables
3174  integer :: isd, ied, jsd, jed
3175  integer :: IsdB, IedB, JsdB, JedB
3176  integer :: IscB, IecB, JscB, JecB
3177  character(len=40) :: mdl = "allocate_OBC_segment_data" ! This subroutine's name.
3178 
3179  isd = segment%HI%isd ; ied = segment%HI%ied
3180  jsd = segment%HI%jsd ; jed = segment%HI%jed
3181  isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
3182  jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
3183  iscb = segment%HI%IscB ; iecb = segment%HI%IecB
3184  jscb = segment%HI%JscB ; jecb = segment%HI%JecB
3185 
3186  if (.not. segment%on_pe) return
3187 
3188  if (segment%is_E_or_W) then
3189  ! If these are just Flather, change update_OBC_segment_data accordingly
3190  allocate(segment%Cg(isdb:iedb,jsd:jed)); segment%Cg(:,:)=0.
3191  allocate(segment%Htot(isdb:iedb,jsd:jed)); segment%Htot(:,:)=0.0
3192  allocate(segment%h(isdb:iedb,jsd:jed,obc%ke)); segment%h(:,:,:)=0.0
3193  allocate(segment%eta(isdb:iedb,jsd:jed)); segment%eta(:,:)=0.0
3194  if (segment%radiation) then
3195  allocate(segment%rx_norm_rad(isdb:iedb,jsd:jed,obc%ke)); segment%rx_norm_rad(:,:,:)=0.0
3196  endif
3197  allocate(segment%normal_vel(isdb:iedb,jsd:jed,obc%ke)); segment%normal_vel(:,:,:)=0.0
3198  allocate(segment%normal_vel_bt(isdb:iedb,jsd:jed)); segment%normal_vel_bt(:,:)=0.0
3199  allocate(segment%normal_trans(isdb:iedb,jsd:jed,obc%ke)); segment%normal_trans(:,:,:)=0.0
3200  if (segment%nudged) then
3201  allocate(segment%nudged_normal_vel(isdb:iedb,jsd:jed,obc%ke)); segment%nudged_normal_vel(:,:,:)=0.0
3202  endif
3203  if (segment%radiation_tan .or. segment%nudged_tan .or. segment%specified_tan .or. &
3204  segment%oblique_tan .or. obc%computed_vorticity .or. obc%computed_strain) then
3205  allocate(segment%tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_vel(:,:,:)=0.0
3206  endif
3207  if (segment%nudged_tan) then
3208  allocate(segment%nudged_tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_vel(:,:,:)=0.0
3209  endif
3210  if (segment%nudged_grad) then
3211  allocate(segment%nudged_tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
3212  endif
3213  if (obc%specified_vorticity .or. obc%specified_strain .or. segment%radiation_grad .or. &
3214  segment%oblique_grad) then
3215  allocate(segment%tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_grad(:,:,:)=0.0
3216  endif
3217  if (segment%oblique) then
3218  allocate(segment%grad_normal(jsdb:jedb,2,obc%ke)); segment%grad_normal(:,:,:) = 0.0
3219  allocate(segment%rx_norm_obl(isdb:iedb,jsd:jed,obc%ke)); segment%rx_norm_obl(:,:,:)=0.0
3220  allocate(segment%ry_norm_obl(isdb:iedb,jsd:jed,obc%ke)); segment%ry_norm_obl(:,:,:)=0.0
3221  allocate(segment%cff_normal(isdb:iedb,jsd:jed,obc%ke)); segment%cff_normal(:,:,:)=0.0
3222  endif
3223  if (segment%oblique_tan) then
3224  allocate(segment%grad_tan(jsd-1:jed+1,2,obc%ke)); segment%grad_tan(:,:,:) = 0.0
3225  endif
3226  if (segment%oblique_grad) then
3227  allocate(segment%grad_gradient(jsd:jed,2,obc%ke)); segment%grad_gradient(:,:,:) = 0.0
3228  endif
3229  endif
3230 
3231  if (segment%is_N_or_S) then
3232  ! If these are just Flather, change update_OBC_segment_data accordingly
3233  allocate(segment%Cg(isd:ied,jsdb:jedb)); segment%Cg(:,:)=0.
3234  allocate(segment%Htot(isd:ied,jsdb:jedb)); segment%Htot(:,:)=0.0
3235  allocate(segment%h(isd:ied,jsdb:jedb,obc%ke)); segment%h(:,:,:)=0.0
3236  allocate(segment%eta(isd:ied,jsdb:jedb)); segment%eta(:,:)=0.0
3237  if (segment%radiation) then
3238  allocate(segment%ry_norm_rad(isd:ied,jsdb:jedb,obc%ke)); segment%ry_norm_rad(:,:,:)=0.0
3239  endif
3240  allocate(segment%normal_vel(isd:ied,jsdb:jedb,obc%ke)); segment%normal_vel(:,:,:)=0.0
3241  allocate(segment%normal_vel_bt(isd:ied,jsdb:jedb)); segment%normal_vel_bt(:,:)=0.0
3242  allocate(segment%normal_trans(isd:ied,jsdb:jedb,obc%ke)); segment%normal_trans(:,:,:)=0.0
3243  if (segment%nudged) then
3244  allocate(segment%nudged_normal_vel(isd:ied,jsdb:jedb,obc%ke)); segment%nudged_normal_vel(:,:,:)=0.0
3245  endif
3246  if (segment%radiation_tan .or. segment%nudged_tan .or. segment%specified_tan .or. &
3247  segment%oblique_tan .or. obc%computed_vorticity .or. obc%computed_strain) then
3248  allocate(segment%tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_vel(:,:,:)=0.0
3249  endif
3250  if (segment%nudged_tan) then
3251  allocate(segment%nudged_tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_vel(:,:,:)=0.0
3252  endif
3253  if (segment%nudged_grad) then
3254  allocate(segment%nudged_tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
3255  endif
3256  if (obc%specified_vorticity .or. obc%specified_strain .or. segment%radiation_grad .or. &
3257  segment%oblique_grad) then
3258  allocate(segment%tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_grad(:,:,:)=0.0
3259  endif
3260  if (segment%oblique) then
3261  allocate(segment%grad_normal(isdb:iedb,2,obc%ke)); segment%grad_normal(:,:,:) = 0.0
3262  allocate(segment%rx_norm_obl(isd:ied,jsdb:jedb,obc%ke)); segment%rx_norm_obl(:,:,:)=0.0
3263  allocate(segment%ry_norm_obl(isd:ied,jsdb:jedb,obc%ke)); segment%ry_norm_obl(:,:,:)=0.0
3264  allocate(segment%cff_normal(isd:ied,jsdb:jedb,obc%ke)); segment%cff_normal(:,:,:)=0.0
3265  endif
3266  if (segment%oblique_tan) then
3267  allocate(segment%grad_tan(isd-1:ied+1,2,obc%ke)); segment%grad_tan(:,:,:) = 0.0
3268  endif
3269  if (segment%oblique_grad) then
3270  allocate(segment%grad_gradient(isd:ied,2,obc%ke)); segment%grad_gradient(:,:,:) = 0.0
3271  endif
3272  endif
3273 

Referenced by setup_u_point_obc(), and setup_v_point_obc().

Here is the caller graph for this function:

◆ deallocate_obc_segment_data()

subroutine mom_open_boundary::deallocate_obc_segment_data ( type(ocean_obc_type), pointer  OBC,
type(obc_segment_type), intent(inout)  segment 
)
private

Deallocate segment data fields.

Parameters
obcOpen boundary structure
[in,out]segmentOpen boundary segment

Definition at line 3278 of file MOM_open_boundary.F90.

3278  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
3279  type(OBC_segment_type), intent(inout) :: segment !< Open boundary segment
3280  ! Local variables
3281  character(len=40) :: mdl = "deallocate_OBC_segment_data" ! This subroutine's name.
3282 
3283  if (.not. segment%on_pe) return
3284 
3285  if (associated (segment%Cg)) deallocate(segment%Cg)
3286  if (associated (segment%Htot)) deallocate(segment%Htot)
3287  if (associated (segment%h)) deallocate(segment%h)
3288  if (associated (segment%eta)) deallocate(segment%eta)
3289  if (associated (segment%rx_norm_rad)) deallocate(segment%rx_norm_rad)
3290  if (associated (segment%ry_norm_rad)) deallocate(segment%ry_norm_rad)
3291  if (associated (segment%rx_norm_obl)) deallocate(segment%rx_norm_obl)
3292  if (associated (segment%ry_norm_obl)) deallocate(segment%ry_norm_obl)
3293  if (associated (segment%cff_normal)) deallocate(segment%cff_normal)
3294  if (associated (segment%grad_normal)) deallocate(segment%grad_normal)
3295  if (associated (segment%grad_tan)) deallocate(segment%grad_tan)
3296  if (associated (segment%grad_gradient)) deallocate(segment%grad_gradient)
3297  if (associated (segment%normal_vel)) deallocate(segment%normal_vel)
3298  if (associated (segment%normal_vel_bt)) deallocate(segment%normal_vel_bt)
3299  if (associated (segment%normal_trans)) deallocate(segment%normal_trans)
3300  if (associated (segment%nudged_normal_vel)) deallocate(segment%nudged_normal_vel)
3301  if (associated (segment%tangential_vel)) deallocate(segment%tangential_vel)
3302  if (associated (segment%nudged_tangential_vel)) deallocate(segment%nudged_tangential_vel)
3303  if (associated (segment%nudged_tangential_grad)) deallocate(segment%nudged_tangential_grad)
3304  if (associated (segment%tangential_grad)) deallocate(segment%tangential_grad)
3305  if (associated (segment%tr_Reg)) call segment_tracer_registry_end(segment%tr_Reg)
3306 
3307 

References segment_tracer_registry_end().

Referenced by open_boundary_dealloc().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ file_obc_end()

subroutine, public mom_open_boundary::file_obc_end ( type(file_obc_cs), pointer  CS)

Clean up the file OBC from registry.

Parameters
csOBC file control structure.

Definition at line 3972 of file MOM_open_boundary.F90.

3972  type(file_OBC_CS), pointer :: CS !< OBC file control structure.
3973 
3974  if (associated(cs)) then
3975  deallocate(cs)
3976  endif

Referenced by mom_boundary_update::obc_register_end().

Here is the caller graph for this function:

◆ fill_temp_salt_segments()

subroutine, public mom_open_boundary::fill_temp_salt_segments ( type(ocean_grid_type), intent(inout)  G,
type(ocean_obc_type), pointer  OBC,
type(thermo_var_ptrs), intent(inout)  tv 
)
Parameters
[in,out]gOcean grid structure
obcOpen boundary structure
[in,out]tvThermodynamics structure

Definition at line 4122 of file MOM_open_boundary.F90.

4122  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
4123  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
4124  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
4125 
4126 ! Local variables
4127  integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n, nz
4128  integer :: i, j, k
4129  type(OBC_segment_type), pointer :: segment => null() ! pointer to segment type list
4130 
4131  if (.not. associated(obc)) return
4132  if (.not. associated(tv%T) .and. associated(tv%S)) return
4133  ! Both temperature and salinity fields
4134 
4135  call pass_var(tv%T, g%Domain)
4136  call pass_var(tv%S, g%Domain)
4137 
4138  nz = g%ke
4139 
4140  do n=1, obc%number_of_segments
4141  segment => obc%segment(n)
4142  if (.not. segment%on_pe) cycle
4143 
4144  isd = segment%HI%isd ; ied = segment%HI%ied
4145  jsd = segment%HI%jsd ; jed = segment%HI%jed
4146  isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
4147  jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
4148 
4149  ! Fill with T and S values
4150  if (segment%is_E_or_W) then
4151  i=segment%HI%IsdB
4152  do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
4153  if (segment%direction == obc_direction_w) then
4154  segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i+1,j,k)
4155  segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i+1,j,k)
4156  else
4157  segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j,k)
4158  segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j,k)
4159  endif
4160  enddo ; enddo
4161  else
4162  j=segment%HI%JsdB
4163  do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
4164  if (segment%direction == obc_direction_s) then
4165  segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j+1,k)
4166  segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j+1,k)
4167  else
4168  segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j,k)
4169  segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j,k)
4170  endif
4171  enddo ; enddo
4172  endif
4173  segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:)
4174  segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:)
4175  enddo
4176  call setup_obc_tracer_reservoirs(g, obc)

References obc_direction_s, obc_direction_w, and setup_obc_tracer_reservoirs().

Here is the call graph for this function:

◆ flood_fill()

subroutine mom_open_boundary::flood_fill ( type(dyn_horgrid_type), intent(inout)  G,
real, dimension(:,:), intent(inout)  color,
integer, intent(in)  cin,
integer, intent(in)  cout,
integer, intent(in)  cland 
)
private

flood the cin, cout values

Parameters
[in,out]gOcean grid structure
[in,out]colorFor sorting inside from outside
[in]cincolor for inside the domain
[in]coutcolor for outside the domain
[in]clandcolor for inside the land mask

Definition at line 4292 of file MOM_open_boundary.F90.

4292  type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure
4293  real, dimension(:,:), intent(inout) :: color !< For sorting inside from outside
4294  integer, intent(in) :: cin !< color for inside the domain
4295  integer, intent(in) :: cout !< color for outside the domain
4296  integer, intent(in) :: cland !< color for inside the land mask
4297 
4298 ! Local variables
4299  integer :: i, j, ncount
4300 
4301  ncount = 1
4302  do while (ncount > 0)
4303  ncount = 0
4304  do j=g%jsd+1,g%jed-1
4305  do i=g%isd+1,g%ied-1
4306  if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0) then
4307  color(i,j) = color(i-1,j)
4308  ncount = ncount + 1
4309  endif
4310  if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0) then
4311  color(i,j) = color(i+1,j)
4312  ncount = ncount + 1
4313  endif
4314  if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0) then
4315  color(i,j) = color(i,j-1)
4316  ncount = ncount + 1
4317  endif
4318  if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0) then
4319  color(i,j) = color(i,j+1)
4320  ncount = ncount + 1
4321  endif
4322  enddo
4323  enddo
4324  do j=g%jed-1,g%jsd+1,-1
4325  do i=g%ied-1,g%isd+1,-1
4326  if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0) then
4327  color(i,j) = color(i-1,j)
4328  ncount = ncount + 1
4329  endif
4330  if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0) then
4331  color(i,j) = color(i+1,j)
4332  ncount = ncount + 1
4333  endif
4334  if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0) then
4335  color(i,j) = color(i,j-1)
4336  ncount = ncount + 1
4337  endif
4338  if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0) then
4339  color(i,j) = color(i,j+1)
4340  ncount = ncount + 1
4341  endif
4342  enddo
4343  enddo
4344  call pass_var(color, g%Domain)
4345  call sum_across_pes(ncount)
4346  enddo
4347 

Referenced by mask_outside_obcs().

Here is the caller graph for this function:

◆ flood_fill2()

subroutine mom_open_boundary::flood_fill2 ( type(dyn_horgrid_type), intent(inout)  G,
real, dimension(:,:), intent(inout)  color,
integer, intent(in)  cin,
integer, intent(in)  cout,
integer, intent(in)  cland 
)
private

flood the cin, cout values

Parameters
[in,out]gOcean grid structure
[in,out]colorFor sorting inside from outside
[in]cincolor for inside the domain
[in]coutcolor for outside the domain
[in]clandcolor for inside the land mask

Definition at line 4352 of file MOM_open_boundary.F90.

4352  type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure
4353  real, dimension(:,:), intent(inout) :: color !< For sorting inside from outside
4354  integer, intent(in) :: cin !< color for inside the domain
4355  integer, intent(in) :: cout !< color for outside the domain
4356  integer, intent(in) :: cland !< color for inside the land mask
4357 
4358 ! Local variables
4359  integer :: i, j, ncount
4360 
4361  ncount = 1
4362  do while (ncount > 0)
4363  ncount = 0
4364  do i=g%isd+1,g%ied-1
4365  do j=g%jsd+1,g%jed-1
4366  if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0) then
4367  color(i,j) = color(i-1,j)
4368  ncount = ncount + 1
4369  endif
4370  if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0) then
4371  color(i,j) = color(i+1,j)
4372  ncount = ncount + 1
4373  endif
4374  if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0) then
4375  color(i,j) = color(i,j-1)
4376  ncount = ncount + 1
4377  endif
4378  if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0) then
4379  color(i,j) = color(i,j+1)
4380  ncount = ncount + 1
4381  endif
4382  enddo
4383  enddo
4384  do i=g%ied-1,g%isd+1,-1
4385  do j=g%jed-1,g%jsd+1,-1
4386  if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0) then
4387  color(i,j) = color(i-1,j)
4388  ncount = ncount + 1
4389  endif
4390  if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0) then
4391  color(i,j) = color(i+1,j)
4392  ncount = ncount + 1
4393  endif
4394  if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0) then
4395  color(i,j) = color(i,j-1)
4396  ncount = ncount + 1
4397  endif
4398  if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0) then
4399  color(i,j) = color(i,j+1)
4400  ncount = ncount + 1
4401  endif
4402  enddo
4403  enddo
4404  call pass_var(color, g%Domain)
4405  call sum_across_pes(ncount)
4406  enddo
4407 

Referenced by mask_outside_obcs().

Here is the caller graph for this function:

◆ gradient_at_q_points()

subroutine mom_open_boundary::gradient_at_q_points ( type(ocean_grid_type), intent(in)  G,
type(obc_segment_type), pointer  segment,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  uvel,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  vvel 
)
private

Calculate the tangential gradient of the normal flow at the boundary q-points.

Parameters
[in]gOcean grid structure
segmentOBC segment structure
[in]uvelzonal velocity [L T-1 ~> m s-1]
[in]vvelmeridional velocity [L T-1 ~> m s-1]

Definition at line 2966 of file MOM_open_boundary.F90.

2966  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
2967  type(OBC_segment_type), pointer :: segment !< OBC segment structure
2968  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uvel !< zonal velocity [L T-1 ~> m s-1]
2969  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vvel !< meridional velocity [L T-1 ~> m s-1]
2970  integer :: i,j,k
2971 
2972  if (.not. segment%on_pe) return
2973 
2974  if (segment%is_E_or_W) then
2975  if (segment%direction == obc_direction_e) then
2976  i=segment%HI%isdB
2977  do k=1,g%ke
2978  do j=max(segment%HI%JsdB, g%HI%JsdB+1),min(segment%HI%JedB, g%HI%JedB-1)
2979  segment%grad_normal(j,1,k) = (uvel(i-1,j+1,k)-uvel(i-1,j,k)) * g%mask2dBu(i-1,j)
2980  segment%grad_normal(j,2,k) = (uvel(i,j+1,k)-uvel(i,j,k)) * g%mask2dBu(i,j)
2981  enddo
2982  enddo
2983  if (segment%oblique_tan) then
2984  do k=1,g%ke
2985  do j=max(segment%HI%jsd-1, g%HI%jsd),min(segment%HI%jed+1, g%HI%jed)
2986  segment%grad_tan(j,1,k) = (vvel(i-1,j,k)-vvel(i-1,j-1,k)) * g%mask2dT(i-1,j)
2987  segment%grad_tan(j,2,k) = (vvel(i,j,k)-vvel(i,j-1,k)) * g%mask2dT(i,j)
2988  enddo
2989  enddo
2990  endif
2991  if (segment%oblique_grad) then
2992  do k=1,g%ke
2993  do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
2994  segment%grad_gradient(j,1,k) = (((vvel(i-1,j,k) - vvel(i-2,j,k))*g%IdxBu(i-2,j)) - &
2995  (vvel(i-1,j-1,k) - vvel(i-2,j-1,k))*g%IdxBu(i-2,j-1)) * g%mask2dCu(i-2,j)
2996  segment%grad_gradient(j,2,k) = (((vvel(i,j,k) - vvel(i-1,j,k))*g%IdxBu(i-1,j)) - &
2997  (vvel(i,j-1,k) - vvel(i-1,j-1,k))*g%IdxBu(i-1,j-1)) * g%mask2dCu(i-1,j)
2998  enddo
2999  enddo
3000  endif
3001  else ! western segment
3002  i=segment%HI%isdB
3003  do k=1,g%ke
3004  do j=max(segment%HI%JsdB, g%HI%JsdB+1),min(segment%HI%JedB, g%HI%JedB-1)
3005  segment%grad_normal(j,1,k) = (uvel(i+1,j+1,k)-uvel(i+1,j,k)) * g%mask2dBu(i+1,j)
3006  segment%grad_normal(j,2,k) = (uvel(i,j+1,k)-uvel(i,j,k)) * g%mask2dBu(i,j)
3007  enddo
3008  enddo
3009  if (segment%oblique_tan) then
3010  do k=1,g%ke
3011  do j=max(segment%HI%jsd-1, g%HI%jsd),min(segment%HI%jed+1, g%HI%jed)
3012  segment%grad_tan(j,1,k) = (vvel(i+2,j,k)-vvel(i+2,j-1,k)) * g%mask2dT(i+2,j)
3013  segment%grad_tan(j,2,k) = (vvel(i+1,j,k)-vvel(i+1,j-1,k)) * g%mask2dT(i+1,j)
3014  enddo
3015  enddo
3016  endif
3017  if (segment%oblique_grad) then
3018  do k=1,g%ke
3019  do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
3020  segment%grad_gradient(j,1,k) = (((vvel(i+3,j,k) - vvel(i+2,j,k))*g%IdxBu(i+2,j)) - &
3021  (vvel(i+3,j-1,k) - vvel(i+2,j-1,k))*g%IdxBu(i+2,j-1)) * g%mask2dCu(i+2,j)
3022  segment%grad_gradient(j,2,k) = (((vvel(i+2,j,k) - vvel(i+1,j,k))*g%IdxBu(i+1,j)) - &
3023  (vvel(i+2,j-1,k) - vvel(i+1,j-1,k))*g%IdxBu(i+1,j-1)) * g%mask2dCu(i+1,j)
3024  enddo
3025  enddo
3026  endif
3027  endif
3028  elseif (segment%is_N_or_S) then
3029  if (segment%direction == obc_direction_n) then
3030  j=segment%HI%jsdB
3031  do k=1,g%ke
3032  do i=max(segment%HI%IsdB, g%HI%IsdB+1),min(segment%HI%IedB, g%HI%IedB-1)
3033  segment%grad_normal(i,1,k) = (vvel(i+1,j-1,k)-vvel(i,j-1,k)) * g%mask2dBu(i,j-1)
3034  segment%grad_normal(i,2,k) = (vvel(i+1,j,k)-vvel(i,j,k)) * g%mask2dBu(i,j)
3035  enddo
3036  enddo
3037  if (segment%oblique_tan) then
3038  do k=1,g%ke
3039  do i=max(segment%HI%isd-1, g%HI%isd),min(segment%HI%ied+1, g%HI%ied)
3040  segment%grad_tan(i,1,k) = (uvel(i,j-1,k)-uvel(i-1,j-1,k)) * g%mask2dT(i,j-1)
3041  segment%grad_tan(i,2,k) = (uvel(i,j,k)-uvel(i-1,j,k)) * g%mask2dT(i,j)
3042  enddo
3043  enddo
3044  endif
3045  if (segment%oblique_grad) then
3046  do k=1,g%ke
3047  do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
3048  segment%grad_gradient(i,1,k) = (((uvel(i,j-1,k) - uvel(i,j-2,k))*g%IdyBu(i,j-2)) - &
3049  (uvel(i-1,j-1,k) - uvel(i-1,j-2,k))*g%IdyBu(i-1,j-2)) * g%mask2dCv(i,j-2)
3050  segment%grad_gradient(i,2,k) = (((uvel(i,j,k) - uvel(i,j-1,k))*g%IdyBu(i,j-1)) - &
3051  (uvel(i-1,j,k) - uvel(i-1,j-1,k))*g%IdyBu(i-1,j-1)) * g%mask2dCv(i,j-1)
3052  enddo
3053  enddo
3054  endif
3055  else ! south segment
3056  j=segment%HI%jsdB
3057  do k=1,g%ke
3058  do i=max(segment%HI%IsdB, g%HI%IsdB+1),min(segment%HI%IedB, g%HI%IedB-1)
3059  segment%grad_normal(i,1,k) = (vvel(i+1,j+1,k)-vvel(i,j+1,k)) * g%mask2dBu(i,j+1)
3060  segment%grad_normal(i,2,k) = (vvel(i+1,j,k)-vvel(i,j,k)) * g%mask2dBu(i,j)
3061  enddo
3062  enddo
3063  if (segment%oblique_tan) then
3064  do k=1,g%ke
3065  do i=max(segment%HI%isd-1, g%HI%isd),min(segment%HI%ied+1, g%HI%ied)
3066  segment%grad_tan(i,1,k) = (uvel(i,j+2,k)-uvel(i-1,j+2,k)) * g%mask2dT(i,j+2)
3067  segment%grad_tan(i,2,k) = (uvel(i,j+1,k)-uvel(i-1,j+1,k)) * g%mask2dT(i,j+1)
3068  enddo
3069  enddo
3070  endif
3071  if (segment%oblique_grad) then
3072  do k=1,g%ke
3073  do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
3074  segment%grad_gradient(i,1,k) = (((uvel(i,j+3,k) - uvel(i,j+2,k))*g%IdyBu(i,j+2)) - &
3075  (uvel(i-1,j+3,k) - uvel(i-1,j+2,k))*g%IdyBu(i-1,j+2)) * g%mask2dCv(i,j+2)
3076  segment%grad_gradient(i,2,k) = (((uvel(i,j+2,k) - uvel(i,j+1,k))*g%IdyBu(i,j+1)) - &
3077  (uvel(i-1,j+2,k) - uvel(i-1,j+1,k))*g%IdyBu(i-1,j+1)) * g%mask2dCv(i,j+1)
3078  enddo
3079  enddo
3080  endif
3081  endif
3082  endif
3083 

References obc_direction_e, and obc_direction_n.

Referenced by radiation_open_bdry_conds().

Here is the caller graph for this function:

◆ initialize_segment_data()

subroutine mom_open_boundary::initialize_segment_data ( type(dyn_horgrid_type), intent(in)  G,
type(ocean_obc_type), intent(inout)  OBC,
type(param_file_type), intent(in)  PF 
)
private

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.

Parameters
[in]gOcean grid structure
[in,out]obcOpen boundary control structure
[in]pfParameter file handle

Definition at line 541 of file MOM_open_boundary.F90.

541  use mpp_mod, only : mpp_pe, mpp_set_current_pelist, mpp_get_current_pelist,mpp_npes
542 
543  type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure
544  type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure
545  type(param_file_type), intent(in) :: PF !< Parameter file handle
546 
547  integer :: n,m,num_fields
548  character(len=256) :: segstr, filename
549  character(len=20) :: segnam, suffix
550  character(len=32) :: varnam, fieldname
551  real :: value
552  character(len=32), dimension(MAX_OBC_FIELDS) :: fields ! segment field names
553  character(len=128) :: inputdir
554  type(OBC_segment_type), pointer :: segment => null() ! pointer to segment type list
555  character(len=32) :: remappingScheme
556  character(len=256) :: mesg ! Message for error messages.
557  logical :: check_reconstruction, check_remapping, force_bounds_in_subcell
558  integer, dimension(4) :: siz,siz2
559  integer :: is, ie, js, je
560  integer :: isd, ied, jsd, jed
561  integer :: IsdB, IedB, JsdB, JedB
562  integer, dimension(:), allocatable :: saved_pelist
563  integer :: current_pe
564  integer, dimension(1) :: single_pelist
565  !will be able to dynamically switch between sub-sampling refined grid data or model grid
566 
567  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
568 
569  ! There is a problem with the order of the OBC initialization
570  ! with respect to ALE_init. Currently handling this by copying the
571  ! param file so that I can use it later in step_MOM in order to finish
572  ! initializing segments on the first step.
573 
574  call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
575  inputdir = slasher(inputdir)
576 
577  call get_param(pf, mdl, "REMAPPING_SCHEME", remappingscheme, &
578  "This sets the reconstruction scheme used "//&
579  "for vertical remapping for all variables. "//&
580  "It can be one of the following schemes: \n"//&
581  trim(remappingschemesdoc), default=remappingdefaultscheme,do_not_log=.true.)
582  call get_param(pf, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
583  "If true, cell-by-cell reconstructions are checked for "//&
584  "consistency and if non-monotonicity or an inconsistency is "//&
585  "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
586  call get_param(pf, mdl, "FATAL_CHECK_REMAPPING", check_remapping, &
587  "If true, the results of remapping are checked for "//&
588  "conservation and new extrema and if an inconsistency is "//&
589  "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
590  call get_param(pf, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
591  "If true, the values on the intermediate grid used for remapping "//&
592  "are forced to be bounded, which might not be the case due to "//&
593  "round off.", default=.false.,do_not_log=.true.)
594  call get_param(pf, mdl, "BRUSHCUTTER_MODE", obc%brushcutter_mode, &
595  "If true, read external OBC data on the supergrid.", &
596  default=.false.)
597 
598  allocate(obc%remap_CS)
599  call initialize_remapping(obc%remap_CS, remappingscheme, boundary_extrapolation = .false., &
600  check_reconstruction=check_reconstruction, &
601  check_remapping=check_remapping, force_bounds_in_subcell=force_bounds_in_subcell)
602 
603  if (obc%user_BCs_set_globally) return
604 
605  ! Try this here just for the documentation. It is repeated below.
606  do n=1, obc%number_of_segments
607  segment => obc%segment(n)
608  write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n
609  call get_param(pf, mdl, segnam, segstr, 'OBC segment docs')
610  enddo
611 
612  !< temporarily disable communication in order to read segment data independently
613 
614  allocate(saved_pelist(0:mpp_npes()-1))
615  call mpp_get_current_pelist(saved_pelist)
616  current_pe = mpp_pe()
617  single_pelist(1) = current_pe
618  call mpp_set_current_pelist(single_pelist)
619 
620  do n=1, obc%number_of_segments
621  segment => obc%segment(n)
622  if (.not. segment%values_needed) cycle
623 
624  write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n
625  write(suffix,"('_segment_',i3.3)") n
626  ! needs documentation !! Yet, unsafe for now, causes grief for
627  ! MOM_parameter_docs in circle_obcs on two processes.
628 ! call get_param(PF, mdl, segnam, segstr, 'xyz')
629  ! Clear out any old values
630  segstr = ''
631  call get_param(pf, mdl, segnam, segstr)
632  if (segstr == '') then
633  write(mesg,'("No OBC_SEGMENT_XXX_DATA string for OBC segment ",I3)') n
634  call mom_error(fatal, mesg)
635  endif
636 
637  call parse_segment_data_str(trim(segstr), fields=fields, num_fields=num_fields)
638  if (num_fields == 0) then
639  call mom_mesg('initialize_segment_data: num_fields = 0')
640  cycle ! cycle to next segment
641  endif
642 
643  allocate(segment%field(num_fields))
644  segment%num_fields = num_fields
645 
646  segment%temp_segment_data_exists=.false.
647  segment%salt_segment_data_exists=.false.
648 !!
649 ! CODE HERE FOR OTHER OPTIONS (CLAMPED, NUDGED,..)
650 !!
651 
652  isd = segment%HI%isd ; ied = segment%HI%ied
653  jsd = segment%HI%jsd ; jed = segment%HI%jed
654  isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
655  jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
656 
657  do m=1,num_fields
658  call parse_segment_data_str(trim(segstr), var=trim(fields(m)), value=value, filenam=filename, fieldnam=fieldname)
659  if (trim(filename) /= 'none') then
660  obc%update_OBC = .true. ! Data is assumed to be time-dependent if we are reading from file
661  obc%needs_IO_for_data = .true. ! At least one segment is using I/O for OBC data
662 ! segment%values_needed = .true. ! Indicates that i/o will be needed for this segment
663  segment%field(m)%name = trim(fields(m))
664  if (segment%field(m)%name == 'TEMP') then
665  segment%temp_segment_data_exists=.true.
666  segment%t_values_needed = .false.
667  endif
668  if (segment%field(m)%name == 'SALT') then
669  segment%salt_segment_data_exists=.true.
670  segment%s_values_needed = .false.
671  endif
672  filename = trim(inputdir)//trim(filename)
673  fieldname = trim(fieldname)//trim(suffix)
674  call field_size(filename,fieldname,siz,no_domain=.true.)
675 ! if (siz(4) == 1) segment%values_needed = .false.
676  if (segment%on_pe) then
677  if (obc%brushcutter_mode .and. (modulo(siz(1),2) == 0 .or. modulo(siz(2),2) == 0)) then
678  call mom_error(fatal,'segment data are not on the supergrid')
679  endif
680  siz2(1)=1
681 
682  if (siz(1)>1) then
683  if (obc%brushcutter_mode) then
684  siz2(1)=(siz(1)-1)/2
685  else
686  siz2(1)=siz(1)
687  endif
688  endif
689  siz2(2)=1
690  if (siz(2)>1) then
691  if (obc%brushcutter_mode) then
692  siz2(2)=(siz(2)-1)/2
693  else
694  siz2(2)=siz(2)
695  endif
696  endif
697  siz2(3)=siz(3)
698 
699  if (segment%is_E_or_W) then
700  if (segment%field(m)%name == 'V') then
701  allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
702  segment%v_values_needed = .false.
703  else if (segment%field(m)%name == 'DVDX') then
704  allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
705  segment%g_values_needed = .false.
706  else
707  allocate(segment%field(m)%buffer_src(isdb:iedb,jsd:jed,siz2(3)))
708  if (segment%field(m)%name == 'U') then
709  segment%u_values_needed = .false.
710  else if (segment%field(m)%name == 'SSH') then
711  segment%z_values_needed = .false.
712  else if (segment%field(m)%name == 'TEMP') then
713  segment%t_values_needed = .false.
714  else if (segment%field(m)%name == 'SALT') then
715  segment%s_values_needed = .false.
716  endif
717  endif
718  else
719  if (segment%field(m)%name == 'U') then
720  allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
721  segment%u_values_needed = .false.
722  else if (segment%field(m)%name == 'DUDY') then
723  allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
724  segment%g_values_needed = .false.
725  else
726  allocate(segment%field(m)%buffer_src(isd:ied,jsdb:jedb,siz2(3)))
727  if (segment%field(m)%name == 'V') then
728  segment%v_values_needed = .false.
729  else if (segment%field(m)%name == 'SSH') then
730  segment%z_values_needed = .false.
731  else if (segment%field(m)%name == 'TEMP') then
732  segment%t_values_needed = .false.
733  else if (segment%field(m)%name == 'SALT') then
734  segment%s_values_needed = .false.
735  endif
736  endif
737  endif
738  segment%field(m)%buffer_src(:,:,:)=0.0
739  segment%field(m)%fid = init_external_field(trim(filename),&
740  trim(fieldname),ignore_axis_atts=.true.,threading=single_file)
741  if (siz(3) > 1) then
742  fieldname = 'dz_'//trim(fieldname)
743  call field_size(filename,fieldname,siz,no_domain=.true.)
744  if (segment%is_E_or_W) then
745  if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
746  allocate(segment%field(m)%dz_src(isdb:iedb,jsdb:jedb,siz(3)))
747  else
748  allocate(segment%field(m)%dz_src(isdb:iedb,jsd:jed,siz(3)))
749  endif
750  else
751  if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
752  allocate(segment%field(m)%dz_src(isdb:iedb,jsdb:jedb,siz(3)))
753  else
754  allocate(segment%field(m)%dz_src(isd:ied,jsdb:jedb,siz(3)))
755  endif
756  endif
757  segment%field(m)%dz_src(:,:,:)=0.0
758  segment%field(m)%nk_src=siz(3)
759  segment%field(m)%fid_dz = init_external_field(trim(filename),trim(fieldname),&
760  ignore_axis_atts=.true.,threading=single_file)
761  else
762  segment%field(m)%nk_src=1
763  endif
764  endif
765  else
766  segment%field(m)%fid = -1
767  segment%field(m)%value = value
768  segment%field(m)%name = trim(fields(m))
769  if (segment%field(m)%name == 'U') then
770  segment%u_values_needed = .false.
771  elseif (segment%field(m)%name == 'V') then
772  segment%v_values_needed = .false.
773  elseif (segment%field(m)%name == 'SSH') then
774  segment%z_values_needed = .false.
775  elseif (segment%field(m)%name == 'TEMP') then
776  segment%t_values_needed = .false.
777  elseif (segment%field(m)%name == 'SALT') then
778  segment%s_values_needed = .false.
779  elseif (segment%field(m)%name == 'DVDX' .or. segment%field(m)%name == 'DUDY') then
780  segment%g_values_needed = .false.
781  endif
782  endif
783  enddo
784  if (segment%u_values_needed .or. segment%v_values_needed .or. &
785  segment%t_values_needed .or. segment%s_values_needed .or. &
786  segment%z_values_needed .or. segment%g_values_needed) then
787  write(mesg,'("Values needed for OBC segment ",I3)') n
788  call mom_error(fatal, mesg)
789  endif
790  enddo
791 
792  call mpp_set_current_pelist(saved_pelist)
793 

References mom_remapping::initialize_remapping(), mdl, mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), and parse_segment_data_str().

Referenced by open_boundary_config().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ lookup_seg_field()

integer function mom_open_boundary::lookup_seg_field ( type(obc_segment_type), pointer  OBC_seg,
character(len=32), intent(in)  field 
)
private

Needs documentation.

Parameters
obc_segOBC segment
[in]fieldThe field name

Definition at line 3152 of file MOM_open_boundary.F90.

3152  type(OBC_segment_type), pointer :: OBC_seg !< OBC segment
3153  character(len=32), intent(in) :: field !< The field name
3154  integer :: lookup_seg_field
3155  ! Local variables
3156  integer :: n
3157 
3158  lookup_seg_field=-1
3159  do n=1,obc_seg%num_fields
3160  if (trim(field) == obc_seg%field(n)%name) then
3161  lookup_seg_field=n
3162  return
3163  endif
3164  enddo
3165 

◆ mask_outside_obcs()

subroutine mom_open_boundary::mask_outside_obcs ( type(dyn_horgrid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(ocean_obc_type), pointer  OBC 
)
private

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

Parameters
[in,out]gOcean grid structure
[in]param_fileParameter file handle
obcOpen boundary structure
[in]usA dimensional unit scaling type

Definition at line 4183 of file MOM_open_boundary.F90.

4183  type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure
4184  type(param_file_type), intent(in) :: param_file !< Parameter file handle
4185  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
4186  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
4187 
4188  ! Local variables
4189  integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n
4190  integer :: i, j
4191  logical :: fatal_error = .false.
4192  real :: min_depth
4193  integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2
4194  character(len=256) :: mesg ! Message for error messages.
4195  type(OBC_segment_type), pointer :: segment => null() ! pointer to segment type list
4196  real, allocatable, dimension(:,:) :: color, color2 ! For sorting inside from outside,
4197  ! two different ways
4198 
4199  if (.not. associated(obc)) return
4200 
4201  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
4202  units="m", default=0.0, scale=us%m_to_Z, do_not_log=.true.)
4203 
4204  allocate(color(g%isd:g%ied, g%jsd:g%jed)) ; color = 0
4205  allocate(color2(g%isd:g%ied, g%jsd:g%jed)) ; color2 = 0
4206 
4207 
4208  ! Paint a frame around the outside.
4209  do j=g%jsd,g%jed
4210  color(g%isd,j) = cedge
4211  color(g%ied,j) = cedge
4212  color2(g%isd,j) = cedge
4213  color2(g%ied,j) = cedge
4214  enddo
4215  do i=g%isd,g%ied
4216  color(i,g%jsd) = cedge
4217  color(i,g%jed) = cedge
4218  color2(i,g%jsd) = cedge
4219  color2(i,g%jed) = cedge
4220  enddo
4221 
4222  ! Set color to cland in the land. Note that this is before the land
4223  ! mask has been initialized, set mask values based on depth.
4224  do j=g%jsd,g%jed
4225  do i=g%isd,g%ied
4226  if (g%bathyT(i,j) <= min_depth) color(i,j) = cland
4227  if (g%bathyT(i,j) <= min_depth) color2(i,j) = cland
4228  enddo
4229  enddo
4230 
4231  do j=g%jsd,g%jed ; do i=g%IsdB+1,g%IedB-1
4232  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
4233  if (color(i,j) == 0.0) color(i,j) = cout
4234  if (color(i+1,j) == 0.0) color(i+1,j) = cin
4235  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
4236  if (color(i,j) == 0.0) color(i,j) = cin
4237  if (color(i+1,j) == 0.0) color(i+1,j) = cout
4238  endif
4239  enddo ; enddo
4240  do j=g%JsdB+1,g%JedB-1 ; do i=g%isd,g%ied
4241  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
4242  if (color(i,j) == 0.0) color(i,j) = cout
4243  if (color(i,j+1) == 0.0) color(i,j+1) = cin
4244  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
4245  if (color(i,j) == 0.0) color(i,j) = cin
4246  if (color(i,j+1) == 0.0) color(i,j+1) = cout
4247  endif
4248  enddo ; enddo
4249 
4250  do j=g%JsdB+1,g%JedB-1 ; do i=g%isd,g%ied
4251  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
4252  if (color2(i,j) == 0.0) color2(i,j) = cout
4253  if (color2(i,j+1) == 0.0) color2(i,j+1) = cin
4254  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
4255  if (color2(i,j) == 0.0) color2(i,j) = cin
4256  if (color2(i,j+1) == 0.0) color2(i,j+1) = cout
4257  endif
4258  enddo ; enddo
4259  do j=g%jsd,g%jed ; do i=g%IsdB+1,g%IedB-1
4260  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
4261  if (color2(i,j) == 0.0) color2(i,j) = cout
4262  if (color2(i+1,j) == 0.0) color2(i+1,j) = cin
4263  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
4264  if (color2(i,j) == 0.0) color2(i,j) = cin
4265  if (color2(i+1,j) == 0.0) color2(i+1,j) = cout
4266  endif
4267  enddo ; enddo
4268 
4269  ! Do the flood fill until there are no more uncolored cells.
4270  call flood_fill(g, color, cin, cout, cland)
4271  call flood_fill2(g, color2, cin, cout, cland)
4272 
4273  ! Use the color to set outside to min_depth on this process.
4274  do j=g%jsd,g%jed ; do i=g%isd,g%ied
4275  if (color(i,j) /= color2(i,j)) then
4276  fatal_error = .true.
4277  write(mesg,'("MOM_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", &
4278  "the masking of the outside grid points.")') i, j
4279  call mom_error(warning,"MOM register_tracer: "//mesg, all_print=.true.)
4280  endif
4281  if (color(i,j) == cout) g%bathyT(i,j) = min_depth
4282  enddo ; enddo
4283  if (fatal_error) call mom_error(fatal, &
4284  "MOM_open_boundary: inconsistent OBC segments.")
4285 
4286  deallocate(color)
4287  deallocate(color2)

References flood_fill(), flood_fill2(), mdl, obc_direction_e, obc_direction_n, obc_direction_s, and obc_direction_w.

Referenced by open_boundary_config().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ obc_registry_init()

subroutine, public mom_open_boundary::obc_registry_init ( type(param_file_type), intent(in)  param_file,
type(obc_registry_type), pointer  Reg 
)

This routine include declares and sets the variable "version".

Parameters
[in]param_fileopen file to parse for model parameters
regpointer to OBC registry

Definition at line 3925 of file MOM_open_boundary.F90.

3925  type(param_file_type), intent(in) :: param_file !< open file to parse for model parameters
3926  type(OBC_registry_type), pointer :: Reg !< pointer to OBC registry
3927 
3928  integer, save :: init_calls = 0
3929 
3930 #include "version_variable.h"
3931  character(len=40) :: mdl = "MOM_open_boundary" ! This module's name.
3932  character(len=256) :: mesg ! Message for error messages.
3933 
3934  if (.not.associated(reg)) then ; allocate(reg)
3935  else ; return ; endif
3936 
3937  ! Read all relevant parameters and write them to the model log.
3938 ! call log_version(param_file, mdl,s version, "")
3939 
3940  init_calls = init_calls + 1
3941  if (init_calls > 1) then
3942  write(mesg,'("OBC_registry_init called ",I3, &
3943  &" times with different registry pointers.")') init_calls
3944  if (is_root_pe()) call mom_error(warning,"MOM_open_boundary"//mesg)
3945  endif
3946 

Referenced by register_obc().

Here is the caller graph for this function:

◆ open_boundary_apply_normal_flow()

subroutine, public mom_open_boundary::open_boundary_apply_normal_flow ( type(ocean_obc_type), pointer  OBC,
type(ocean_grid_type), intent(inout)  G,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)  v 
)

Applies OBC values stored in segments to 3d u,v fields.

Parameters
obcOpen boundary control structure
[in,out]gOcean grid structure
[in,out]uu field to update on open boundaries [L T-1 ~> m s-1]
[in,out]vv field to update on open boundaries [L T-1 ~> m s-1]

Definition at line 2898 of file MOM_open_boundary.F90.

2898  ! Arguments
2899  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
2900  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
2901  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< u field to update on open
2902  !! boundaries [L T-1 ~> m s-1]
2903  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< v field to update on open
2904  !! boundaries [L T-1 ~> m s-1]
2905  ! Local variables
2906  integer :: i, j, k, n
2907  type(OBC_segment_type), pointer :: segment => null()
2908 
2909  if (.not.associated(obc)) return ! Bail out if OBC is not available
2910 
2911  do n=1,obc%number_of_segments
2912  segment => obc%segment(n)
2913  if (.not. segment%on_pe) then
2914  cycle
2915  elseif (segment%radiation .or. segment%oblique .or. segment%gradient) then
2916  if (segment%is_E_or_W) then
2917  i=segment%HI%IsdB
2918  do k=1,g%ke ; do j=segment%HI%jsd,segment%HI%jed
2919  u(i,j,k) = segment%normal_vel(i,j,k)
2920  enddo ; enddo
2921  elseif (segment%is_N_or_S) then
2922  j=segment%HI%JsdB
2923  do k=1,g%ke ; do i=segment%HI%isd,segment%HI%ied
2924  v(i,j,k) = segment%normal_vel(i,j,k)
2925  enddo ; enddo
2926  endif
2927  endif
2928  enddo
2929 

Referenced by radiation_open_bdry_conds().

Here is the caller graph for this function:

◆ open_boundary_config()

subroutine, public mom_open_boundary::open_boundary_config ( type(dyn_horgrid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(ocean_obc_type), pointer  OBC 
)

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.

Parameters
[in,out]gOcean grid structure
[in]usA dimensional unit scaling type
[in]param_fileParameter file handle
obcOpen boundary control structure

Definition at line 317 of file MOM_open_boundary.F90.

317  type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure
318  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
319  type(param_file_type), intent(in) :: param_file !< Parameter file handle
320  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
321  ! Local variables
322  integer :: l ! For looping over segments
323  logical :: debug_OBC, debug, mask_outside, reentrant_x, reentrant_y
324  character(len=15) :: segment_param_str ! The run-time parameter name for each segment
325  character(len=100) :: segment_str ! The contents (rhs) for parameter "segment_param_str"
326  character(len=200) :: config1 ! String for OBC_USER_CONFIG
327  real :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m]
328  allocate(obc)
329 
330  call log_version(param_file, mdl, version, &
331  "Controls where open boundaries are located, what kind of boundary condition "//&
332  "to impose, and what data to apply, if any.")
333  call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", obc%number_of_segments, &
334  "The number of open boundary segments.", &
335  default=0)
336  call get_param(param_file, mdl, "G_EARTH", obc%g_Earth, &
337  "The gravitational acceleration of the Earth.", &
338  units="m s-2", default = 9.80)
339  call get_param(param_file, mdl, "OBC_USER_CONFIG", config1, &
340  "A string that sets how the open boundary conditions are "//&
341  " configured: \n", default="none", do_not_log=.true.)
342  call get_param(param_file, mdl, "NK", obc%ke, &
343  "The number of model layers", default=0, do_not_log=.true.)
344 
345  if (config1 /= "none" .and. config1 /= "dyed_obcs") obc%user_BCs_set_globally = .true.
346 
347  if (obc%number_of_segments > 0) then
348  call get_param(param_file, mdl, "OBC_ZERO_VORTICITY", obc%zero_vorticity, &
349  "If true, sets relative vorticity to zero on open boundaries.", &
350  default=.false.)
351  call get_param(param_file, mdl, "OBC_FREESLIP_VORTICITY", obc%freeslip_vorticity, &
352  "If true, sets the normal gradient of tangential velocity to "//&
353  "zero in the relative vorticity on open boundaries. This cannot "//&
354  "be true if another OBC_XXX_VORTICITY option is True.", default=.true.)
355  call get_param(param_file, mdl, "OBC_COMPUTED_VORTICITY", obc%computed_vorticity, &
356  "If true, uses the external values of tangential velocity "//&
357  "in the relative vorticity on open boundaries. This cannot "//&
358  "be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
359  call get_param(param_file, mdl, "OBC_SPECIFIED_VORTICITY", obc%specified_vorticity, &
360  "If true, uses the external values of tangential velocity "//&
361  "in the relative vorticity on open boundaries. This cannot "//&
362  "be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
363  if ((obc%zero_vorticity .and. obc%freeslip_vorticity) .or. &
364  (obc%zero_vorticity .and. obc%computed_vorticity) .or. &
365  (obc%zero_vorticity .and. obc%specified_vorticity) .or. &
366  (obc%freeslip_vorticity .and. obc%computed_vorticity) .or. &
367  (obc%freeslip_vorticity .and. obc%specified_vorticity) .or. &
368  (obc%computed_vorticity .and. obc%specified_vorticity)) &
369  call mom_error(fatal, "MOM_open_boundary.F90, open_boundary_config:\n"//&
370  "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"//&
371  "and OBC_IMPORTED_VORTICITY can be True at once.")
372  call get_param(param_file, mdl, "OBC_ZERO_STRAIN", obc%zero_strain, &
373  "If true, sets the strain used in the stress tensor to zero on open boundaries.", &
374  default=.false.)
375  call get_param(param_file, mdl, "OBC_FREESLIP_STRAIN", obc%freeslip_strain, &
376  "If true, sets the normal gradient of tangential velocity to "//&
377  "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
378  "be true if another OBC_XXX_STRAIN option is True.", default=.true.)
379  call get_param(param_file, mdl, "OBC_COMPUTED_STRAIN", obc%computed_strain, &
380  "If true, sets the normal gradient of tangential velocity to "//&
381  "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
382  "be true if another OBC_XXX_STRAIN option is True.", default=.false.)
383  call get_param(param_file, mdl, "OBC_SPECIFIED_STRAIN", obc%specified_strain, &
384  "If true, sets the normal gradient of tangential velocity to "//&
385  "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
386  "be true if another OBC_XXX_STRAIN option is True.", default=.false.)
387  if ((obc%zero_strain .and. obc%freeslip_strain) .or. &
388  (obc%zero_strain .and. obc%computed_strain) .or. &
389  (obc%zero_strain .and. obc%specified_strain) .or. &
390  (obc%freeslip_strain .and. obc%computed_strain) .or. &
391  (obc%freeslip_strain .and. obc%specified_strain) .or. &
392  (obc%computed_strain .and. obc%specified_strain)) &
393  call mom_error(fatal, "MOM_open_boundary.F90, open_boundary_config: \n"//&
394  "Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN, OBC_COMPUTED_STRAIN \n"//&
395  "and OBC_IMPORTED_STRAIN can be True at once.")
396  call get_param(param_file, mdl, "OBC_ZERO_BIHARMONIC", obc%zero_biharmonic, &
397  "If true, zeros the Laplacian of flow on open boundaries in the biharmonic "//&
398  "viscosity term.", default=.false.)
399  call get_param(param_file, mdl, "MASK_OUTSIDE_OBCS", mask_outside, &
400  "If true, set the areas outside open boundaries to be land.", &
401  default=.false.)
402  call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
403  call get_param(param_file, mdl, "DEBUG_OBC", debug_obc, default=.false.)
404  if (debug_obc .or. debug) &
405  call log_param(param_file, mdl, "DEBUG_OBC", debug_obc, &
406  "If true, do additional calls to help debug the performance "//&
407  "of the open boundary condition code.", default=.false., &
408  debuggingparam=.true.)
409 
410  call get_param(param_file, mdl, "OBC_SILLY_THICK", obc%silly_h, &
411  "A silly value of thicknesses used outside of open boundary "//&
412  "conditions for debugging.", units="m", default=0.0, scale=us%m_to_Z, &
413  do_not_log=.not.debug_obc, debuggingparam=.true.)
414  call get_param(param_file, mdl, "OBC_SILLY_VEL", obc%silly_u, &
415  "A silly value of velocities used outside of open boundary "//&
416  "conditions for debugging.", units="m/s", default=0.0, scale=us%m_s_to_L_T, &
417  do_not_log=.not.debug_obc, debuggingparam=.true.)
418  reentrant_x = .false.
419  call get_param(param_file, mdl, "REENTRANT_X", reentrant_x, default=.true.)
420  reentrant_y = .false.
421  call get_param(param_file, mdl, "REENTRANT_Y", reentrant_y, default=.false.)
422 
423  ! Allocate everything
424  ! Note the 0-segment is needed when %segnum_u/v(:,:) = 0
425  allocate(obc%segment(0:obc%number_of_segments))
426  do l=0,obc%number_of_segments
427  obc%segment(l)%Flather = .false.
428  obc%segment(l)%radiation = .false.
429  obc%segment(l)%radiation_tan = .false.
430  obc%segment(l)%radiation_grad = .false.
431  obc%segment(l)%oblique = .false.
432  obc%segment(l)%oblique_tan = .false.
433  obc%segment(l)%oblique_grad = .false.
434  obc%segment(l)%nudged = .false.
435  obc%segment(l)%nudged_tan = .false.
436  obc%segment(l)%nudged_grad = .false.
437  obc%segment(l)%specified = .false.
438  obc%segment(l)%specified_tan = .false.
439  obc%segment(l)%open = .false.
440  obc%segment(l)%gradient = .false.
441  obc%segment(l)%values_needed = .false.
442  obc%segment(l)%u_values_needed = .false.
443  obc%segment(l)%v_values_needed = .false.
444  obc%segment(l)%t_values_needed = .false.
445  obc%segment(l)%s_values_needed = .false.
446  obc%segment(l)%z_values_needed = .false.
447  obc%segment(l)%g_values_needed = .false.
448  obc%segment(l)%direction = obc_none
449  obc%segment(l)%is_N_or_S = .false.
450  obc%segment(l)%is_E_or_W = .false.
451  obc%segment(l)%is_E_or_W_2 = .false.
452  obc%segment(l)%Velocity_nudging_timescale_in = 0.0
453  obc%segment(l)%Velocity_nudging_timescale_out = 0.0
454  obc%segment(l)%num_fields = 0
455  enddo
456  allocate(obc%segnum_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; obc%segnum_u(:,:) = obc_none
457  allocate(obc%segnum_v(g%isd:g%ied,g%JsdB:g%JedB)) ; obc%segnum_v(:,:) = obc_none
458 
459  do l = 1, obc%number_of_segments
460  write(segment_param_str(1:15),"('OBC_SEGMENT_',i3.3)") l
461  call get_param(param_file, mdl, segment_param_str, segment_str, &
462  "Documentation needs to be dynamic?????", &
463  fail_if_missing=.true.)
464  segment_str = remove_spaces(segment_str)
465  if (segment_str(1:2) == 'I=') then
466  call setup_u_point_obc(obc, g, us, segment_str, l, param_file, reentrant_y)
467  elseif (segment_str(1:2) == 'J=') then
468  call setup_v_point_obc(obc, g, us, segment_str, l, param_file, reentrant_x)
469  else
470  call mom_error(fatal, "MOM_open_boundary.F90, open_boundary_config: "//&
471  "Unable to interpret "//segment_param_str//" = "//trim(segment_str))
472  endif
473  enddo
474 
475  ! if (open_boundary_query(OBC, needs_ext_seg_data=.true.)) &
476  call initialize_segment_data(g, obc, param_file)
477 
478  if (open_boundary_query(obc, apply_open_obc=.true.)) then
479  call get_param(param_file, mdl, "OBC_RADIATION_MAX", obc%rx_max, &
480  "The maximum magnitude of the baroclinic radiation "//&
481  "velocity (or speed of characteristics). This is only "//&
482  "used if one of the open boundary segments is using Orlanski.", &
483  units="m s-1", default=10.0)
484  call get_param(param_file, mdl, "OBC_RAD_VEL_WT", obc%gamma_uv, &
485  "The relative weighting for the baroclinic radiation "//&
486  "velocities (or speed of characteristics) at the new "//&
487  "time level (1) or the running mean (0) for velocities. "//&
488  "Valid values range from 0 to 1. This is only used if "//&
489  "one of the open boundary segments is using Orlanski.", &
490  units="nondim", default=0.3)
491  endif
492 
493  lscale_in = 0.
494  lscale_out = 0.
495  if (open_boundary_query(obc, apply_open_obc=.true.)) then
496  call get_param(param_file, mdl, "OBC_TRACER_RESERVOIR_LENGTH_SCALE_OUT ", lscale_out, &
497  "An effective length scale for restoring the tracer concentration "//&
498  "at the boundaries to externally imposed values when the flow "//&
499  "is exiting the domain.", units="m", default=0.0, scale=us%m_to_L)
500 
501  call get_param(param_file, mdl, "OBC_TRACER_RESERVOIR_LENGTH_SCALE_IN ", lscale_in, &
502  "An effective length scale for restoring the tracer concentration "//&
503  "at the boundaries to values from the interior when the flow "//&
504  "is entering the domain.", units="m", default=0.0, scale=us%m_to_L)
505  endif
506 
507  if (mask_outside) call mask_outside_obcs(g, us, param_file, obc)
508 
509  ! All tracers are using the same restoring length scale for now, but we may want to make this
510  ! tracer-specific in the future for example, in cases where certain tracers are poorly constrained
511  ! by data while others are well constrained - MJH.
512  do l = 1, obc%number_of_segments
513  obc%segment(l)%Tr_InvLscale_in = 0.0
514  if (lscale_in>0.) obc%segment(l)%Tr_InvLscale_in = 1.0/lscale_in
515  obc%segment(l)%Tr_InvLscale_out = 0.0
516  if (lscale_out>0.) obc%segment(l)%Tr_InvLscale_out = 1.0/lscale_out
517  enddo
518 
519  endif ! OBC%number_of_segments > 0
520 
521  ! Safety check
522  if ((obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally) .and. &
523  .not.g%symmetric ) call mom_error(fatal, &
524  "MOM_open_boundary, open_boundary_config: "//&
525  "Symmetric memory must be used when using Flather OBCs.")
526 
527  if (.not.(obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
528  obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) then
529  ! No open boundaries have been requested
530  call open_boundary_dealloc(obc)
531  else
532  ! Need this for ocean_only mode boundary interpolation.
533  call time_interp_external_init()
534  endif
535 

References initialize_segment_data(), mask_outside_obcs(), mdl, mom_error_handler::mom_error(), obc_none, open_boundary_dealloc(), open_boundary_query(), mom_string_functions::remove_spaces(), setup_u_point_obc(), and setup_v_point_obc().

Here is the call graph for this function:

◆ open_boundary_dealloc()

subroutine mom_open_boundary::open_boundary_dealloc ( type(ocean_obc_type), pointer  OBC)
private

Deallocate open boundary data.

Parameters
obcOpen boundary control structure

Definition at line 1563 of file MOM_open_boundary.F90.

1563  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
1564  type(OBC_segment_type), pointer :: segment => null()
1565  integer :: n
1566 
1567  if (.not. associated(obc)) return
1568 
1569  do n=1, obc%number_of_segments
1570  segment => obc%segment(n)
1571  call deallocate_obc_segment_data(obc, segment)
1572  enddo
1573  if (associated(obc%segment)) deallocate(obc%segment)
1574  if (associated(obc%segnum_u)) deallocate(obc%segnum_u)
1575  if (associated(obc%segnum_v)) deallocate(obc%segnum_v)
1576  if (associated(obc%rx_normal)) deallocate(obc%rx_normal)
1577  if (associated(obc%ry_normal)) deallocate(obc%ry_normal)
1578  if (associated(obc%rx_oblique)) deallocate(obc%rx_oblique)
1579  if (associated(obc%ry_oblique)) deallocate(obc%ry_oblique)
1580  if (associated(obc%cff_normal)) deallocate(obc%cff_normal)
1581  if (associated(obc%tres_x)) deallocate(obc%tres_x)
1582  if (associated(obc%tres_y)) deallocate(obc%tres_y)
1583  deallocate(obc)

References deallocate_obc_segment_data().

Referenced by open_boundary_config(), and open_boundary_end().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ open_boundary_end()

subroutine, public mom_open_boundary::open_boundary_end ( type(ocean_obc_type), pointer  OBC)

Close open boundary data.

Parameters
obcOpen boundary control structure

Definition at line 1588 of file MOM_open_boundary.F90.

1588  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
1589  call open_boundary_dealloc(obc)

References open_boundary_dealloc().

Here is the call graph for this function:

◆ open_boundary_impose_land_mask()

subroutine, public mom_open_boundary::open_boundary_impose_land_mask ( type(ocean_obc_type), pointer  OBC,
type(dyn_horgrid_type), intent(inout)  G,
real, dimension(szib_(g),szj_(g)), intent(inout)  areaCu,
real, dimension(szi_(g),szjb_(g)), intent(inout)  areaCv,
type(unit_scale_type), intent(in)  US 
)

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.

Parameters
obcOpen boundary control structure
[in,out]gOcean grid structure
[in]usA dimensional unit scaling type
[in,out]areacuArea of a u-cell [L2 ~> m2]
[in,out]areacvArea of a u-cell [L2 ~> m2]

Definition at line 1639 of file MOM_open_boundary.F90.

1639  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
1640  type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure
1641  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1642  real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: areaCu !< Area of a u-cell [L2 ~> m2]
1643  real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: areaCv !< Area of a u-cell [L2 ~> m2]
1644  ! Local variables
1645  integer :: i, j, n
1646  type(OBC_segment_type), pointer :: segment => null()
1647  logical :: any_U, any_V
1648 
1649  if (.not.associated(obc)) return
1650 
1651  do n=1,obc%number_of_segments
1652  segment=>obc%segment(n)
1653  if (.not. segment%on_pe) cycle
1654  if (segment%is_E_or_W) then
1655  ! Sweep along u-segments and delete the OBC for blocked points.
1656  ! Also, mask all points outside.
1657  i=segment%HI%IsdB
1658  do j=segment%HI%jsd,segment%HI%jed
1659  if (g%mask2dCu(i,j) == 0) obc%segnum_u(i,j) = obc_none
1660  if (segment%direction == obc_direction_w) then
1661  g%mask2dT(i,j) = 0
1662  else
1663  g%mask2dT(i+1,j) = 0
1664  endif
1665  enddo
1666  do j=segment%HI%JsdB+1,segment%HI%JedB-1
1667  if (segment%direction == obc_direction_w) then
1668  g%mask2dCv(i,j) = 0
1669  else
1670  g%mask2dCv(i+1,j) = 0
1671  endif
1672  enddo
1673  else
1674  ! Sweep along v-segments and delete the OBC for blocked points.
1675  j=segment%HI%JsdB
1676  do i=segment%HI%isd,segment%HI%ied
1677  if (g%mask2dCv(i,j) == 0) obc%segnum_v(i,j) = obc_none
1678  if (segment%direction == obc_direction_s) then
1679  g%mask2dT(i,j) = 0
1680  else
1681  g%mask2dT(i,j+1) = 0
1682  endif
1683  enddo
1684  do i=segment%HI%IsdB+1,segment%HI%IedB-1
1685  if (segment%direction == obc_direction_s) then
1686  g%mask2dCu(i,j) = 0
1687  else
1688  g%mask2dCu(i,j+1) = 0
1689  endif
1690  enddo
1691  endif
1692  enddo
1693 
1694  do n=1,obc%number_of_segments
1695  segment=>obc%segment(n)
1696  if (.not. segment%on_pe .or. .not. segment%specified) cycle
1697  if (segment%is_E_or_W) then
1698  ! Sweep along u-segments and for %specified BC points reset the u-point area which was masked out
1699  i=segment%HI%IsdB
1700  do j=segment%HI%jsd,segment%HI%jed
1701  if (segment%direction == obc_direction_e) then
1702  areacu(i,j) = g%areaT(i,j) ! Both of these are in [L2]
1703  else ! West
1704  areacu(i,j) = g%areaT(i+1,j) ! Both of these are in [L2]
1705  endif
1706  enddo
1707  else
1708  ! Sweep along v-segments and for %specified BC points reset the v-point area which was masked out
1709  j=segment%HI%JsdB
1710  do i=segment%HI%isd,segment%HI%ied
1711  if (segment%direction == obc_direction_s) then
1712  areacv(i,j) = g%areaT(i,j+1) ! Both of these are in [L2]
1713  else ! North
1714  areacu(i,j) = g%areaT(i,j) ! Both of these are in [L2]
1715  endif
1716  enddo
1717  endif
1718  enddo
1719 
1720  ! G%mask2du will be open wherever bathymetry allows it.
1721  ! Bathymetry outside of the open boundary was adjusted to match
1722  ! the bathymetry inside so these points will be open unless the
1723  ! bathymetry inside the boundary was too shallow and flagged as land.
1724  any_u = .false.
1725  any_v = .false.
1726  do n=1,obc%number_of_segments
1727  segment=>obc%segment(n)
1728  if (.not. segment%on_pe) cycle
1729  if (segment%is_E_or_W) then
1730  i=segment%HI%IsdB
1731  do j=segment%HI%jsd,segment%HI%jed
1732  if (obc%segnum_u(i,j) /= obc_none) any_u = .true.
1733  enddo
1734  else
1735  j=segment%HI%JsdB
1736  do i=segment%HI%isd,segment%HI%ied
1737  if (obc%segnum_v(i,j) /= obc_none) any_v = .true.
1738  enddo
1739  endif
1740  enddo
1741 
1742  obc%OBC_pe = .true.
1743  if (.not.(any_u .or. any_v)) obc%OBC_pe = .false.
1744 

References obc_direction_e, obc_direction_s, obc_direction_w, and obc_none.

Referenced by mom_fixed_initialization::mom_initialize_fixed().

Here is the caller graph for this function:

◆ open_boundary_impose_normal_slope()

subroutine, public mom_open_boundary::open_boundary_impose_normal_slope ( type(ocean_obc_type), pointer  OBC,
type(dyn_horgrid_type), intent(in)  G,
real, dimension(szi_(g),szj_(g)), intent(inout)  depth 
)

Sets the slope of bathymetry normal to an open bounndary to zero.

Parameters
obcOpen boundary control structure
[in]gOcean grid structure
[in,out]depthBathymetry at h-points

Definition at line 1594 of file MOM_open_boundary.F90.

1594  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
1595  type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure
1596  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points
1597  ! Local variables
1598  integer :: i, j, n
1599  type(OBC_segment_type), pointer :: segment => null()
1600 
1601  if (.not.associated(obc)) return
1602 
1603  if (.not.(obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
1604  obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
1605  return
1606 
1607  do n=1,obc%number_of_segments
1608  segment=>obc%segment(n)
1609  if (.not. segment%on_pe) cycle
1610  if (segment%direction == obc_direction_e) then
1611  i=segment%HI%IsdB
1612  do j=segment%HI%jsd,segment%HI%jed
1613  depth(i+1,j) = depth(i,j)
1614  enddo
1615  elseif (segment%direction == obc_direction_w) then
1616  i=segment%HI%IsdB
1617  do j=segment%HI%jsd,segment%HI%jed
1618  depth(i,j) = depth(i+1,j)
1619  enddo
1620  elseif (segment%direction == obc_direction_n) then
1621  j=segment%HI%JsdB
1622  do i=segment%HI%isd,segment%HI%ied
1623  depth(i,j+1) = depth(i,j)
1624  enddo
1625  elseif (segment%direction == obc_direction_s) then
1626  j=segment%HI%JsdB
1627  do i=segment%HI%isd,segment%HI%ied
1628  depth(i,j) = depth(i,j+1)
1629  enddo
1630  endif
1631  enddo
1632 

References obc_direction_e, obc_direction_n, obc_direction_s, and obc_direction_w.

◆ open_boundary_init()

subroutine, public mom_open_boundary::open_boundary_init ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(ocean_obc_type), pointer  OBC,
type(mom_restart_cs), pointer  restart_CSp 
)

Initialize open boundary control structure and do any necessary rescaling of OBC fields that have been read from a restart file.

Parameters
[in]gOcean grid structure
[in]gvContainer for vertical grid information
[in]usA dimensional unit scaling type
[in]param_fileParameter file handle
obcOpen boundary control structure
restart_cspRestart structure, data intent(inout)

Definition at line 1479 of file MOM_open_boundary.F90.

1479  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1480  type(verticalGrid_type), intent(in) :: GV !< Container for vertical grid information
1481  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1482  type(param_file_type), intent(in) :: param_file !< Parameter file handle
1483  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
1484  type(MOM_restart_CS), pointer :: restart_CSp !< Restart structure, data intent(inout)
1485 
1486  ! Local variables
1487  real :: vel2_rescale ! A rescaling factor for squared velocities from the representation in
1488  ! a restart file to the internal representation in this run.
1489  integer :: i, j, k, isd, ied, jsd, jed, nz
1490  integer :: IsdB, IedB, JsdB, JedB
1491  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
1492  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1493 
1494  if (.not.associated(obc)) return
1495 
1496  id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=clock_routine)
1497 
1498  ! The rx_normal and ry_normal arrays used with radiation OBCs are currently in units of grid
1499  ! points per timestep, but if this were to be corrected to [L T-1 ~> m s-1] or [T-1 ~> s-1] to
1500  ! permit timesteps to change between calls to the OBC code, the following would be needed:
1501 ! if ( OBC%radiation_BCs_exist_globally .and. (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. &
1502 ! ((US%m_to_L * US%s_to_T_restart) /= (US%m_to_L_restart * US%s_to_T)) ) then
1503 ! vel_rescale = (US%m_to_L * US%s_to_T_restart) / (US%m_to_L_restart * US%s_to_T)
1504 ! if (query_initialized(OBC%rx_normal, "rx_normal", restart_CSp)) then
1505 ! do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB
1506 ! OBC%rx_normal(I,j,k) = vel_rescale * OBC%rx_normal(I,j,k)
1507 ! enddo ; enddo ; enddo
1508 ! endif
1509 ! if (query_initialized(OBC%ry_normal, "ry_normal", restart_CSp)) then
1510 ! do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied
1511 ! OBC%ry_normal(i,J,k) = vel_rescale * OBC%ry_normal(i,J,k)
1512 ! enddo ; enddo ; enddo
1513 ! endif
1514 ! endif
1515 
1516  ! The oblique boundary condition terms have units of [L2 T-2 ~> m2 s-2] and may need to be rescaled.
1517  if ( obc%oblique_BCs_exist_globally .and. (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1518  ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) ) then
1519  vel2_rescale = (us%m_to_L * us%s_to_T_restart)**2 / (us%m_to_L_restart * us%s_to_T)**2
1520  if (query_initialized(obc%rx_oblique, "rx_oblique", restart_csp)) then
1521  do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
1522  obc%rx_oblique(i,j,k) = vel2_rescale * obc%rx_oblique(i,j,k)
1523  enddo ; enddo ; enddo
1524  endif
1525  if (query_initialized(obc%ry_oblique, "ry_oblique", restart_csp)) then
1526  do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
1527  obc%ry_oblique(i,j,k) = vel2_rescale * obc%ry_oblique(i,j,k)
1528  enddo ; enddo ; enddo
1529  endif
1530  if (query_initialized(obc%cff_normal, "cff_normal", restart_csp)) then
1531  do k=1,nz ; do j=jsdb,jedb ; do i=isdb,iedb
1532  obc%cff_normal(i,j,k) = vel2_rescale * obc%cff_normal(i,j,k)
1533  enddo ; enddo ; enddo
1534  endif
1535  endif
1536 

References id_clock_pass.

◆ open_boundary_query()

logical function, public mom_open_boundary::open_boundary_query ( type(ocean_obc_type), pointer  OBC,
logical, intent(in), optional  apply_open_OBC,
logical, intent(in), optional  apply_specified_OBC,
logical, intent(in), optional  apply_Flather_OBC,
logical, intent(in), optional  apply_nudged_OBC,
logical, intent(in), optional  needs_ext_seg_data 
)
Parameters
obcOpen boundary control structure
[in]apply_open_obcReturns True if open_*_BCs_exist_globally is true
[in]apply_specified_obcReturns True if specified_*_BCs_exist_globally is true
[in]apply_flather_obcReturns True if Flather_*_BCs_exist_globally is true
[in]apply_nudged_obcReturns True if nudged_*_BCs_exist_globally is true
[in]needs_ext_seg_dataReturns True if external segment data needed

Definition at line 1541 of file MOM_open_boundary.F90.

1541  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
1542  logical, optional, intent(in) :: apply_open_OBC !< Returns True if open_*_BCs_exist_globally is true
1543  logical, optional, intent(in) :: apply_specified_OBC !< Returns True if specified_*_BCs_exist_globally is true
1544  logical, optional, intent(in) :: apply_Flather_OBC !< Returns True if Flather_*_BCs_exist_globally is true
1545  logical, optional, intent(in) :: apply_nudged_OBC !< Returns True if nudged_*_BCs_exist_globally is true
1546  logical, optional, intent(in) :: needs_ext_seg_data !< Returns True if external segment data needed
1547  open_boundary_query = .false.
1548  if (.not. associated(obc)) return
1549  if (present(apply_open_obc)) open_boundary_query = obc%open_u_BCs_exist_globally .or. &
1550  obc%open_v_BCs_exist_globally
1551  if (present(apply_specified_obc)) open_boundary_query = obc%specified_u_BCs_exist_globally .or. &
1552  obc%specified_v_BCs_exist_globally
1553  if (present(apply_flather_obc)) open_boundary_query = obc%Flather_u_BCs_exist_globally .or. &
1554  obc%Flather_v_BCs_exist_globally
1555  if (present(apply_nudged_obc)) open_boundary_query = obc%nudged_u_BCs_exist_globally .or. &
1556  obc%nudged_v_BCs_exist_globally
1557  if (present(needs_ext_seg_data)) open_boundary_query = obc%needs_IO_for_data
1558 

Referenced by open_boundary_config().

Here is the caller graph for this function:

◆ open_boundary_register_restarts()

subroutine, public mom_open_boundary::open_boundary_register_restarts ( type(hor_index_type), intent(in)  HI,
type(verticalgrid_type), pointer  GV,
type(ocean_obc_type), pointer  OBC,
type(tracer_registry_type), pointer  Reg,
type(param_file_type), intent(in)  param_file,
type(mom_restart_cs), pointer  restart_CSp,
logical, intent(in)  use_temperature 
)

Register OBC segment data for restarts.

Parameters
[in]hiHorizontal indices
gvContainer for vertical grid information
obcOBC data structure, data intent(inout)
regpointer to tracer registry
[in]param_fileParameter file handle
restart_cspRestart structure, data intent(inout)
[in]use_temperatureIf true, T and S are used

Definition at line 4413 of file MOM_open_boundary.F90.

4413  type(hor_index_type), intent(in) :: HI !< Horizontal indices
4414  type(verticalGrid_type), pointer :: GV !< Container for vertical grid information
4415  type(ocean_OBC_type), pointer :: OBC !< OBC data structure, data intent(inout)
4416  type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry
4417  type(param_file_type), intent(in) :: param_file !< Parameter file handle
4418  type(MOM_restart_CS), pointer :: restart_CSp !< Restart structure, data intent(inout)
4419  logical, intent(in) :: use_temperature !< If true, T and S are used
4420  ! Local variables
4421  type(vardesc) :: vd
4422  integer :: m, n
4423  character(len=100) :: mesg
4424  type(OBC_segment_type), pointer :: segment=>null()
4425 
4426  if (.not. associated(obc)) &
4427  call mom_error(fatal, "open_boundary_register_restarts: Called with "//&
4428  "uninitialized OBC control structure")
4429 
4430  if (associated(obc%rx_normal) .or. associated(obc%ry_normal) .or. &
4431  associated(obc%rx_oblique) .or. associated(obc%ry_oblique) .or. associated(obc%cff_normal)) &
4432  call mom_error(fatal, "open_boundary_register_restarts: Restart "//&
4433  "arrays were previously allocated")
4434 
4435  if (associated(obc%tres_x) .or. associated(obc%tres_y)) &
4436  call mom_error(fatal, "open_boundary_register_restarts: Restart "//&
4437  "arrays were previously allocated")
4438 
4439  ! *** This is a temporary work around for restarts with OBC segments.
4440  ! This implementation uses 3D arrays solely for restarts. We need
4441  ! to be able to add 2D ( x,z or y,z ) data to restarts to avoid using
4442  ! so much memory and disk space. ***
4443  if (obc%radiation_BCs_exist_globally) then
4444  allocate(obc%rx_normal(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke))
4445  obc%rx_normal(:,:,:) = 0.0
4446  vd = var_desc("rx_normal", "m s-1", "Normal Phase Speed for EW radiation OBCs", 'u', 'L')
4447  call register_restart_field(obc%rx_normal, vd, .false., restart_csp)
4448  allocate(obc%ry_normal(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke))
4449  obc%ry_normal(:,:,:) = 0.0
4450  vd = var_desc("ry_normal", "m s-1", "Normal Phase Speed for NS radiation OBCs", 'v', 'L')
4451  call register_restart_field(obc%ry_normal, vd, .false., restart_csp)
4452  endif
4453  if (obc%oblique_BCs_exist_globally) then
4454  allocate(obc%rx_oblique(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke))
4455  obc%rx_oblique(:,:,:) = 0.0
4456  vd = var_desc("rx_oblique", "m2 s-2", "Radiation Speed Squared for EW oblique OBCs", 'u', 'L')
4457  call register_restart_field(obc%rx_oblique, vd, .false., restart_csp)
4458  allocate(obc%ry_oblique(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke))
4459  obc%ry_oblique(:,:,:) = 0.0
4460  vd = var_desc("ry_oblique", "m2 s-2", "Radiation Speed Squared for NS oblique OBCs", 'v', 'L')
4461  call register_restart_field(obc%ry_oblique, vd, .false., restart_csp)
4462  allocate(obc%cff_normal(hi%IsdB:hi%IedB,hi%jsdB:hi%jedB,gv%ke))
4463  obc%cff_normal(:,:,:) = 0.0
4464  vd = var_desc("cff_normal", "m2 s-2", "denominator for oblique OBCs", 'q', 'L')
4465  call register_restart_field(obc%cff_normal, vd, .false., restart_csp)
4466  endif
4467 
4468  if (reg%ntr == 0) return
4469  if (.not. associated(obc%tracer_x_reservoirs_used)) then
4470  obc%ntr = reg%ntr
4471  allocate(obc%tracer_x_reservoirs_used(reg%ntr))
4472  allocate(obc%tracer_y_reservoirs_used(reg%ntr))
4473  obc%tracer_x_reservoirs_used(:) = .false.
4474  obc%tracer_y_reservoirs_used(:) = .false.
4475  call parse_for_tracer_reservoirs(obc, param_file, use_temperature)
4476  else
4477  ! This would be coming from user code such as DOME.
4478  if (obc%ntr /= reg%ntr) then
4479 ! call MOM_error(FATAL, "open_boundary_regiser_restarts: Inconsistent value for ntr")
4480  write(mesg,'("Inconsistent values for ntr ", I8," and ",I8,".")') obc%ntr, reg%ntr
4481  call mom_error(warning, 'open_boundary_register_restarts: '//mesg)
4482  endif
4483  endif
4484 
4485  ! Still painfully inefficient, now in four dimensions.
4486  if (any(obc%tracer_x_reservoirs_used)) then
4487  allocate(obc%tres_x(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke,obc%ntr))
4488  obc%tres_x(:,:,:,:) = 0.0
4489  do m=1,obc%ntr
4490  if (obc%tracer_x_reservoirs_used(m)) then
4491  write(mesg,'("tres_x_",I3.3)') m
4492  vd = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L')
4493  call register_restart_field(obc%tres_x(:,:,:,m), vd, .false., restart_csp)
4494  endif
4495  enddo
4496  endif
4497  if (any(obc%tracer_y_reservoirs_used)) then
4498  allocate(obc%tres_y(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke,obc%ntr))
4499  obc%tres_y(:,:,:,:) = 0.0
4500  do m=1,obc%ntr
4501  if (obc%tracer_y_reservoirs_used(m)) then
4502  write(mesg,'("tres_y_",I3.3)') m
4503  vd = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L')
4504  call register_restart_field(obc%tres_y(:,:,:,m), vd, .false., restart_csp)
4505  endif
4506  enddo
4507  endif
4508 

References parse_for_tracer_reservoirs().

Here is the call graph for this function:

◆ open_boundary_test_extern_h()

subroutine, public mom_open_boundary::open_boundary_test_extern_h ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(ocean_obc_type), pointer  OBC,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(inout)  h 
)

Set thicknesses outside of open boundaries to silly values (used for checking the interior state is independent of values outside of the domain).

Parameters
[in]gOcean grid structure
[in]gvOcean vertical grid structure
obcOpen boundary structure
[in,out]hLayer thickness [H ~> m or kg m-2]

Definition at line 3357 of file MOM_open_boundary.F90.

3357  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
3358  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
3359  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
3360  real, dimension(SZI_(G),SZJ_(G), SZK_(GV)),intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
3361  ! Local variables
3362  real :: silly_h ! A silly thickness for testing [H ~> m or kg m-2]
3363  integer :: i, j, k, n
3364 
3365  if (.not. associated(obc)) return
3366 
3367  silly_h = gv%Z_to_H*obc%silly_h
3368 
3369  do n = 1, obc%number_of_segments
3370  do k = 1, gv%ke
3371  if (obc%segment(n)%is_N_or_S) then
3372  j = obc%segment(n)%HI%JsdB
3373  if (obc%segment(n)%direction == obc_direction_n) then
3374  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
3375  h(i,j+1,k) = silly_h
3376  enddo
3377  else
3378  do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
3379  h(i,j,k) = silly_h
3380  enddo
3381  endif
3382  elseif (obc%segment(n)%is_E_or_W) then
3383  i = obc%segment(n)%HI%IsdB
3384  if (obc%segment(n)%direction == obc_direction_e) then
3385  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
3386  h(i+1,j,k) = silly_h
3387  enddo
3388  else
3389  do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
3390  h(i,j,k) = silly_h
3391  enddo
3392  endif
3393  endif
3394  enddo
3395  enddo
3396 

References obc_direction_e, and obc_direction_n.

Referenced by mom_dynamics_split_rk2::step_mom_dyn_split_rk2().

Here is the caller graph for this function:

◆ open_boundary_test_extern_uv()

subroutine, public mom_open_boundary::open_boundary_test_extern_uv ( type(ocean_grid_type), intent(in)  G,
type(ocean_obc_type), pointer  OBC,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout)  v 
)

Set tangential velocities outside of open boundaries to silly values (used for checking the interior state is independent of values outside of the domain).

Parameters
[in]gOcean grid structure
obcOpen boundary structure
[in,out]uZonal velocity [L T-1 ~> m s-1]
[in,out]vMeridional velocity [L T-1 ~> m s-1]

Definition at line 3314 of file MOM_open_boundary.F90.

3314  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
3315  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
3316  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)),intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
3317  real, dimension(SZI_(G),SZJB_(G), SZK_(G)),intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
3318  ! Local variables
3319  integer :: i, j, k, n
3320 
3321  if (.not. associated(obc)) return
3322 
3323  do n = 1, obc%number_of_segments
3324  do k = 1, g%ke
3325  if (obc%segment(n)%is_N_or_S) then
3326  j = obc%segment(n)%HI%JsdB
3327  if (obc%segment(n)%direction == obc_direction_n) then
3328  do i = obc%segment(n)%HI%IsdB, obc%segment(n)%HI%IedB
3329  u(i,j+1,k) = obc%silly_u
3330  enddo
3331  else
3332  do i = obc%segment(n)%HI%IsdB, obc%segment(n)%HI%IedB
3333  u(i,j,k) = obc%silly_u
3334  enddo
3335  endif
3336  elseif (obc%segment(n)%is_E_or_W) then
3337  i = obc%segment(n)%HI%IsdB
3338  if (obc%segment(n)%direction == obc_direction_e) then
3339  do j = obc%segment(n)%HI%JsdB, obc%segment(n)%HI%JedB
3340  v(i+1,j,k) = obc%silly_u
3341  enddo
3342  else
3343  do j = obc%segment(n)%HI%JsdB, obc%segment(n)%HI%JedB
3344  v(i,j,k) = obc%silly_u
3345  enddo
3346  endif
3347  endif
3348  enddo
3349  enddo
3350 

References obc_direction_e, and obc_direction_n.

◆ open_boundary_zero_normal_flow()

subroutine, public mom_open_boundary::open_boundary_zero_normal_flow ( type(ocean_obc_type), pointer  OBC,
type(ocean_grid_type), intent(inout)  G,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)  v 
)

Applies zero values to 3d u,v fields on OBC segments.

Parameters
obcOpen boundary control structure
[in,out]gOcean grid structure
[in,out]uu field to update on open boundaries
[in,out]vv field to update on open boundaries

Definition at line 2934 of file MOM_open_boundary.F90.

2934  ! Arguments
2935  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
2936  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
2937  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< u field to update on open boundaries
2938  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< v field to update on open boundaries
2939  ! Local variables
2940  integer :: i, j, k, n
2941  type(OBC_segment_type), pointer :: segment => null()
2942 
2943  if (.not.associated(obc)) return ! Bail out if OBC is not available
2944 
2945  do n=1,obc%number_of_segments
2946  segment => obc%segment(n)
2947  if (.not. segment%on_pe) then
2948  cycle
2949  elseif (segment%is_E_or_W) then
2950  i=segment%HI%IsdB
2951  do k=1,g%ke ; do j=segment%HI%jsd,segment%HI%jed
2952  u(i,j,k) = 0.
2953  enddo ; enddo
2954  elseif (segment%is_N_or_S) then
2955  j=segment%HI%JsdB
2956  do k=1,g%ke ; do i=segment%HI%isd,segment%HI%ied
2957  v(i,j,k) = 0.
2958  enddo ; enddo
2959  endif
2960  enddo
2961 

Referenced by mom_dynamics_unsplit::step_mom_dyn_unsplit(), and mom_dynamics_unsplit_rk2::step_mom_dyn_unsplit_rk2().

Here is the caller graph for this function:

◆ parse_for_tracer_reservoirs()

subroutine mom_open_boundary::parse_for_tracer_reservoirs ( type(ocean_obc_type), intent(inout)  OBC,
type(param_file_type), intent(in)  PF,
logical, intent(in)  use_temperature 
)
private

Parse all the OBC_SEGMENT_%%_DATA strings again to see which need tracer reservoirs (all pes need to know).

Parameters
[in,out]obcOpen boundary control structure
[in]pfParameter file handle
[in]use_temperatureIf true, T and S are used

Definition at line 1332 of file MOM_open_boundary.F90.

1332  type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure
1333  type(param_file_type), intent(in) :: PF !< Parameter file handle
1334  logical, intent(in) :: use_temperature !< If true, T and S are used
1335 
1336  ! Local variables
1337  integer :: n,m,num_fields
1338  character(len=256) :: segstr, filename
1339  character(len=20) :: segnam, suffix
1340  character(len=32) :: varnam, fieldname
1341  real :: value
1342  character(len=32), dimension(MAX_OBC_FIELDS) :: fields ! segment field names
1343  type(OBC_segment_type), pointer :: segment => null() ! pointer to segment type list
1344  character(len=256) :: mesg ! Message for error messages.
1345 
1346  do n=1, obc%number_of_segments
1347  segment => obc%segment(n)
1348  write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n
1349  write(suffix,"('_segment_',i3.3)") n
1350  ! Clear out any old values
1351  segstr = ''
1352  call get_param(pf, mdl, segnam, segstr)
1353  if (segstr == '') cycle
1354 
1355  call parse_segment_data_str(trim(segstr), fields=fields, num_fields=num_fields)
1356  if (num_fields == 0) cycle
1357 
1358  ! At this point, just search for TEMP and SALT as tracers 1 and 2.
1359  do m=1,num_fields
1360  call parse_segment_data_str(trim(segstr), var=trim(fields(m)), value=value, filenam=filename, fieldnam=fieldname)
1361  if (trim(filename) /= 'none') then
1362  if (fields(m) == 'TEMP') then
1363  if (segment%is_E_or_W_2) then
1364  obc%tracer_x_reservoirs_used(1) = .true.
1365  else
1366  obc%tracer_y_reservoirs_used(1) = .true.
1367  endif
1368  endif
1369  if (fields(m) == 'SALT') then
1370  if (segment%is_E_or_W_2) then
1371  obc%tracer_x_reservoirs_used(2) = .true.
1372  else
1373  obc%tracer_y_reservoirs_used(2) = .true.
1374  endif
1375  endif
1376  endif
1377  enddo
1378  ! Alternately, set first two to true if use_temperature is true
1379  if (use_temperature) then
1380  if (segment%is_E_or_W_2) then
1381  obc%tracer_x_reservoirs_used(1) = .true.
1382  obc%tracer_x_reservoirs_used(2) = .true.
1383  else
1384  obc%tracer_y_reservoirs_used(1) = .true.
1385  obc%tracer_y_reservoirs_used(2) = .true.
1386  endif
1387  endif
1388  enddo
1389 
1390  return
1391 

References mdl, and parse_segment_data_str().

Referenced by open_boundary_register_restarts().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ parse_segment_data_str()

subroutine mom_open_boundary::parse_segment_data_str ( character(len=*), intent(in)  segment_str,
character(len=*), intent(in), optional  var,
real, intent(out), optional  value,
character(len=*), intent(out), optional  filenam,
character(len=*), intent(out), optional  fieldnam,
character(len=*), dimension(max_obc_fields), intent(out), optional  fields,
integer, intent(out), optional  num_fields,
logical, intent(in), optional  debug 
)
private

Parse an OBC_SEGMENT_%%_DATA string.

Parameters
[in]segment_strA string in form of "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..."
[in]varThe name of the variable for which parameters are needed
[out]filenamThe name of the input file if using "file" method
[out]fieldnamThe name of the variable in the input file if using "file" method
[out]valueA constant value if using the "value" method
[out]fieldsList of fieldnames for each segment
[out]num_fieldsThe number of fields in the segment data
[in]debugIf present and true, write verbose debugging messages

Definition at line 1243 of file MOM_open_boundary.F90.

1243  character(len=*), intent(in) :: segment_str !< A string in form of
1244  !! "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..."
1245  character(len=*), optional, intent(in) :: var !< The name of the variable for which parameters are needed
1246  character(len=*), optional, intent(out) :: filenam !< The name of the input file if using "file" method
1247  character(len=*), optional, intent(out) :: fieldnam !< The name of the variable in the input file if using
1248  !! "file" method
1249  real, optional, intent(out) :: value !< A constant value if using the "value" method
1250  character(len=*), dimension(MAX_OBC_FIELDS), &
1251  optional, intent(out) :: fields !< List of fieldnames for each segment
1252  integer, optional, intent(out) :: num_fields !< The number of fields in the segment data
1253  logical, optional, intent(in) :: debug !< If present and true, write verbose debugging messages
1254  ! Local variables
1255  character(len=128) :: word1, word2, word3, method
1256  integer :: lword, nfields, n, m
1257  logical :: continue,dbg
1258  character(len=32), dimension(MAX_OBC_FIELDS) :: flds
1259 
1260  nfields=0
1261  continue=.true.
1262  dbg=.false.
1263  if (PRESENT(debug)) dbg=debug
1264 
1265  do while (continue)
1266  word1 = extract_word(segment_str,',',nfields+1)
1267  if (trim(word1) == '') exit
1268  nfields=nfields+1
1269  word2 = extract_word(word1,'=',1)
1270  flds(nfields) = trim(word2)
1271  enddo
1272 
1273  if (PRESENT(fields)) then
1274  do n=1,nfields
1275  fields(n) = flds(n)
1276  enddo
1277  endif
1278 
1279  if (PRESENT(num_fields)) then
1280  num_fields=nfields
1281  return
1282  endif
1283 
1284  m=0
1285  if (PRESENT(var)) then
1286  do n=1,nfields
1287  if (trim(var)==trim(flds(n))) then
1288  m=n
1289  exit
1290  endif
1291  enddo
1292  if (m==0) then
1293  call abort()
1294  endif
1295 
1296  ! Process first word which will start with the fieldname
1297  word3 = extract_word(segment_str,',',m)
1298  word1 = extract_word(word3,':',1)
1299 ! if (trim(word1) == '') exit
1300  word2 = extract_word(word1,'=',1)
1301  if (trim(word2) == trim(var)) then
1302  method=trim(extract_word(word1,'=',2))
1303  lword=len_trim(method)
1304  if (method(lword-3:lword) == 'file') then
1305  ! raise an error id filename/fieldname not in argument list
1306  word1 = extract_word(word3,':',2)
1307  filenam = extract_word(word1,'(',1)
1308  fieldnam = extract_word(word1,'(',2)
1309  lword=len_trim(fieldnam)
1310  fieldnam = fieldnam(1:lword-1) ! remove trailing parenth
1311  value=-999.
1312  elseif (method(lword-4:lword) == 'value') then
1313  filenam = 'none'
1314  fieldnam = 'none'
1315  word1 = extract_word(word3,':',2)
1316  lword=len_trim(word1)
1317  read(word1(1:lword),*,end=986,err=987) value
1318  endif
1319  endif
1320  endif
1321 
1322  return
1323  986 call mom_error(fatal,'End of record while parsing segment data specification! '//trim(segment_str))
1324  987 call mom_error(fatal,'Error while parsing segment data specification! '//trim(segment_str))
1325 

References mom_string_functions::extract_word(), and mom_error_handler::mom_error().

Referenced by initialize_segment_data(), and parse_for_tracer_reservoirs().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ parse_segment_param_real()

subroutine mom_open_boundary::parse_segment_param_real ( character(len=*), intent(in)  segment_str,
character(len=*), intent(in)  var,
real, intent(out)  param_value,
logical, intent(in), optional  debug 
)
private

Parse an OBC_SEGMENT_%%_PARAMS string.

Parameters
[in]segment_strA string in form of "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..."
[in]varThe name of the variable for which parameters are needed
[out]param_valueThe value of the parameter
[in]debugIf present and true, write verbose debugging messages

Definition at line 1396 of file MOM_open_boundary.F90.

1396  character(len=*), intent(in) :: segment_str !< A string in form of
1397  !! "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..."
1398  character(len=*), intent(in) :: var !< The name of the variable for which parameters are needed
1399  real, intent(out) :: param_value !< The value of the parameter
1400  logical, optional, intent(in) :: debug !< If present and true, write verbose debugging messages
1401  ! Local variables
1402  character(len=128) :: word1, word2, word3, method
1403  integer :: lword, nfields, n, m
1404  logical :: continue,dbg
1405  character(len=32), dimension(MAX_OBC_FIELDS) :: flds
1406 
1407  nfields=0
1408  continue=.true.
1409  dbg=.false.
1410  if (PRESENT(debug)) dbg=debug
1411 
1412  do while (continue)
1413  word1 = extract_word(segment_str,',',nfields+1)
1414  if (trim(word1) == '') exit
1415  nfields=nfields+1
1416  word2 = extract_word(word1,'=',1)
1417  flds(nfields) = trim(word2)
1418  enddo
1419 
1420  ! if (PRESENT(fields)) then
1421  ! do n=1,nfields
1422  ! fields(n) = flds(n)
1423  ! enddo
1424  ! endif
1425 
1426  ! if (PRESENT(num_fields)) then
1427  ! num_fields=nfields
1428  ! return
1429  ! endif
1430 
1431  m=0
1432 ! if (PRESENT(var)) then
1433  do n=1,nfields
1434  if (trim(var)==trim(flds(n))) then
1435  m=n
1436  exit
1437  endif
1438  enddo
1439  if (m==0) then
1440  call abort()
1441  endif
1442 
1443  ! Process first word which will start with the fieldname
1444  word3 = extract_word(segment_str,',',m)
1445 ! word1 = extract_word(word3,':',1)
1446 ! if (trim(word1) == '') exit
1447  word2 = extract_word(word1,'=',1)
1448  if (trim(word2) == trim(var)) then
1449  method=trim(extract_word(word1,'=',2))
1450  lword=len_trim(method)
1451  read(method(1:lword),*,err=987) param_value
1452  ! if (method(lword-3:lword) == 'file') then
1453  ! ! raise an error id filename/fieldname not in argument list
1454  ! word1 = extract_word(word3,':',2)
1455  ! filenam = extract_word(word1,'(',1)
1456  ! fieldnam = extract_word(word1,'(',2)
1457  ! lword=len_trim(fieldnam)
1458  ! fieldnam = fieldnam(1:lword-1) ! remove trailing parenth
1459  ! value=-999.
1460  ! elseif (method(lword-4:lword) == 'value') then
1461  ! filenam = 'none'
1462  ! fieldnam = 'none'
1463  ! word1 = extract_word(word3,':',2)
1464  ! lword=len_trim(word1)
1465  ! read(word1(1:lword),*,end=986,err=987) value
1466  ! endif
1467  endif
1468 ! endif
1469 
1470  return
1471  986 call mom_error(fatal,'End of record while parsing segment data specification! '//trim(segment_str))
1472  987 call mom_error(fatal,'Error while parsing segment parameter specification! '//trim(segment_str))
1473 

◆ parse_segment_str()

subroutine mom_open_boundary::parse_segment_str ( integer, intent(in)  ni_global,
integer, intent(in)  nj_global,
character(len=*), intent(in)  segment_str,
integer, intent(out)  l,
integer, intent(out)  m,
integer, intent(out)  n,
character(len=*), dimension(:), intent(out)  action_str,
logical, intent(in)  reentrant 
)
private

Parse an OBC_SEGMENT_%%% string.

Parameters
[in]ni_globalNumber of h-points in zonal direction
[in]nj_globalNumber of h-points in meridional direction
[in]segment_strA string in form of "I=l,J=m:n,string" or "J=l,I=m,n,string"
[out]lThe value of I=l, if segment_str begins with I=l, or the value of J=l
[out]mThe value of J=m, if segment_str begins with I=, or the value of I=m
[out]nThe value of J=n, if segment_str begins with I=, or the value of I=n
[out]action_strThe "string" part of segment_str
[in]reentrantis domain reentrant in relevant direction?

Definition at line 1129 of file MOM_open_boundary.F90.

1129  integer, intent(in) :: ni_global !< Number of h-points in zonal direction
1130  integer, intent(in) :: nj_global !< Number of h-points in meridional direction
1131  character(len=*), intent(in) :: segment_str !< A string in form of "I=l,J=m:n,string" or "J=l,I=m,n,string"
1132  integer, intent(out) :: l !< The value of I=l, if segment_str begins with I=l, or the value of J=l
1133  integer, intent(out) :: m !< The value of J=m, if segment_str begins with I=, or the value of I=m
1134  integer, intent(out) :: n !< The value of J=n, if segment_str begins with I=, or the value of I=n
1135  character(len=*), intent(out) :: action_str(:) !< The "string" part of segment_str
1136  logical, intent(in) :: reentrant !< is domain reentrant in relevant direction?
1137  ! Local variables
1138  character(len=24) :: word1, word2, m_word, n_word !< Words delineated by commas in a string in form of
1139  !! "I=%,J=%:%,string"
1140  integer :: l_max !< Either ni_global or nj_global, depending on whether segment_str begins with "I=" or "J="
1141  integer :: mn_max !< Either nj_global or ni_global, depending on whether segment_str begins with "I=" or "J="
1142  integer :: j
1143  integer, parameter :: halo = 10
1144 
1145  ! Process first word which will started with either 'I=' or 'J='
1146  word1 = extract_word(segment_str,',',1)
1147  word2 = extract_word(segment_str,',',2)
1148  if (word1(1:2)=='I=') then
1149  l_max = ni_global
1150  mn_max = nj_global
1151  if (.not. (word2(1:2)=='J=')) call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
1152  "Second word of string '"//trim(segment_str)//"' must start with 'J='.")
1153  elseif (word1(1:2)=='J=') then ! Note that the file_parser uniformaly expands "=" to " = "
1154  l_max = nj_global
1155  mn_max = ni_global
1156  if (.not. (word2(1:2)=='I=')) call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
1157  "Second word of string '"//trim(segment_str)//"' must start with 'I='.")
1158  else
1159  call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str"//&
1160  "String '"//segment_str//"' must start with 'I=' or 'J='.")
1161  endif
1162 
1163  ! Read l
1164  l = interpret_int_expr( word1(3:24), l_max )
1165  if (l<0 .or. l>l_max) then
1166  call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
1167  "First value from string '"//trim(segment_str)//"' is outside of the physical domain.")
1168  endif
1169 
1170  ! Read m
1171  m_word = extract_word(word2(3:24),':',1)
1172  m = interpret_int_expr( m_word, mn_max )
1173  if (reentrant) then
1174  if (m<-halo .or. m>mn_max+halo) then
1175  call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
1176  "Beginning of range in string '"//trim(segment_str)//"' is outside of the physical domain.")
1177  endif
1178  else
1179  if (m<-1 .or. m>mn_max+1) then
1180  call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
1181  "Beginning of range in string '"//trim(segment_str)//"' is outside of the physical domain.")
1182  endif
1183  endif
1184 
1185  ! Read n
1186  n_word = extract_word(word2(3:24),':',2)
1187  n = interpret_int_expr( n_word, mn_max )
1188  if (reentrant) then
1189  if (n<-halo .or. n>mn_max+halo) then
1190  call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
1191  "End of range in string '"//trim(segment_str)//"' is outside of the physical domain.")
1192  endif
1193  else
1194  if (n<-1 .or. n>mn_max+1) then
1195  call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
1196  "End of range in string '"//trim(segment_str)//"' is outside of the physical domain.")
1197  endif
1198  endif
1199 
1200  if (abs(n-m)==0) then
1201  call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str: "//&
1202  "Range in string '"//trim(segment_str)//"' must span one cell.")
1203  endif
1204 
1205  ! Type of open boundary condition
1206  do j = 1, size(action_str)
1207  action_str(j) = extract_word(segment_str,',',2+j)
1208  enddo
1209 
1210  contains
1211 
1212  ! Returns integer value interpreted from string in form of %I, N or N+-%I
1213  integer function interpret_int_expr(string, imax)
1214  character(len=*), intent(in) :: string !< Integer in form or %I, N or N-%I
1215  integer, intent(in) :: imax !< Value to replace 'N' with
1216  ! Local variables
1217  integer slen
1218 
1219  slen = len_trim(string)
1220  if (slen==0) call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str"//&
1221  "Parsed string was empty!")
1222  if (len_trim(string)==1 .and. string(1:1)=='N') then
1223  interpret_int_expr = imax
1224  elseif (string(1:1)=='N') then
1225  if (string(2:2)=='+') then
1226  read(string(3:slen),*,err=911) interpret_int_expr
1228  elseif (string(2:2)=='-') then
1229  read(string(3:slen),*,err=911) interpret_int_expr
1231  endif
1232  else
1233  read(string(1:slen),*,err=911) interpret_int_expr
1234  endif
1235  return
1236  911 call mom_error(fatal, "MOM_open_boundary.F90, parse_segment_str"//&
1237  "Problem reading value from string '"//trim(string)//"'.")
1238  end function interpret_int_expr

References mom_string_functions::extract_word(), interpret_int_expr(), and mom_error_handler::mom_error().

Referenced by setup_u_point_obc(), and setup_v_point_obc().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ radiation_open_bdry_conds()

subroutine, public mom_open_boundary::radiation_open_bdry_conds ( type(ocean_obc_type), pointer  OBC,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout)  u_new,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  u_old,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)  v_new,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v_old,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
real, intent(in)  dt 
)

Apply radiation conditions to 3D u,v at open boundaries.

Parameters
[in,out]gOcean grid structure
obcOpen boundary control structure
[in,out]u_newOn exit, new u values on open boundaries On entry, the old time-level v but including barotropic accelerations [L T-1 ~> m s-1].
[in]u_oldOriginal unadjusted u [L T-1 ~> m s-1]
[in,out]v_newOn exit, new v values on open boundaries. On entry, the old time-level v but including barotropic accelerations [L T-1 ~> m s-1].
[in]v_oldOriginal unadjusted v [L T-1 ~> m s-1]
[in]usA dimensional unit scaling type
[in]dtAppropriate timestep [T ~> s]

Definition at line 1788 of file MOM_open_boundary.F90.

1788  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1789  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
1790  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new !< On exit, new u values on open boundaries
1791  !! On entry, the old time-level v but including
1792  !! barotropic accelerations [L T-1 ~> m s-1].
1793  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u_old !< Original unadjusted u [L T-1 ~> m s-1]
1794  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new !< On exit, new v values on open boundaries.
1795  !! On entry, the old time-level v but including
1796  !! barotropic accelerations [L T-1 ~> m s-1].
1797  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v_old !< Original unadjusted v [L T-1 ~> m s-1]
1798  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1799  real, intent(in) :: dt !< Appropriate timestep [T ~> s]
1800  ! Local variables
1801  real :: dhdt, dhdx, dhdy ! One-point differences in time or space [L T-1 ~> m s-1]
1802  real :: gamma_u, gamma_2 ! Fractional weightings of new values [nondim]
1803  real :: tau ! A local nudging timescale [T ~> s]
1804  real :: rx_max, ry_max ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2]
1805  real :: rx_new, rx_avg ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2]
1806  real :: ry_new, ry_avg ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2]
1807  real :: cff_new, cff_avg ! denominator in oblique [L2 T-2 ~> m2 s-2]
1808  real, allocatable, dimension(:,:,:) :: &
1809  rx_tang_rad, & ! The phase speed at u-points for tangential oblique OBCs
1810  ! in units of grid points per timestep [nondim]
1811  ry_tang_rad, & ! The phase speed at v-points for tangential oblique OBCs
1812  ! in units of grid points per timestep [nondim]
1813  rx_tang_obl, & ! The x-coefficient for tangential oblique OBCs [L2 T-2 ~> m2 s-2]
1814  ry_tang_obl, & ! The y-coefficient for tangential oblique OBCs [L2 T-2 ~> m2 s-2]
1815  cff_tangential ! The denominator for tangential oblique OBCs [L2 T-2 ~> m2 s-2]
1816  real :: eps ! A small velocity squared [L2 T-2 ~> m2 s-2]
1817  type(OBC_segment_type), pointer :: segment => null()
1818  integer :: i, j, k, is, ie, js, je, m, nz, n
1819  integer :: is_obc, ie_obc, js_obc, je_obc
1820 
1821  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1822 
1823  if (.not.associated(obc)) return
1824 
1825  if (.not.(obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
1826  return
1827 
1828  eps = 1.0e-20*us%m_s_to_L_T**2
1829 
1830  !! Copy previously calculated phase velocity from global arrays into segments
1831  !! This is terribly inefficient and temporary solution for continuity across restarts
1832  !! and needs to be revisited in the future.
1833  if (obc%gamma_uv < 1.0) then
1834  do n=1,obc%number_of_segments
1835  segment=>obc%segment(n)
1836  if (.not. segment%on_pe) cycle
1837  if (segment%is_E_or_W .and. segment%radiation) then
1838  do k=1,g%ke
1839  i=segment%HI%IsdB
1840  do j=segment%HI%jsd,segment%HI%jed
1841  segment%rx_norm_rad(i,j,k) = obc%rx_normal(i,j,k)
1842  enddo
1843  enddo
1844  elseif (segment%is_N_or_S .and. segment%radiation) then
1845  do k=1,g%ke
1846  j=segment%HI%JsdB
1847  do i=segment%HI%isd,segment%HI%ied
1848  segment%ry_norm_rad(i,j,k) = obc%ry_normal(i,j,k)
1849  enddo
1850  enddo
1851  endif
1852  if (segment%is_E_or_W .and. segment%oblique) then
1853  do k=1,g%ke
1854  i=segment%HI%IsdB
1855  do j=segment%HI%jsd,segment%HI%jed
1856  segment%rx_norm_obl(i,j,k) = obc%rx_oblique(i,j,k)
1857  segment%ry_norm_obl(i,j,k) = obc%ry_oblique(i,j,k)
1858  segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
1859  enddo
1860  enddo
1861  elseif (segment%is_N_or_S .and. segment%oblique) then
1862  do k=1,g%ke
1863  j=segment%HI%JsdB
1864  do i=segment%HI%isd,segment%HI%ied
1865  segment%rx_norm_obl(i,j,k) = obc%rx_oblique(i,j,k)
1866  segment%ry_norm_obl(i,j,k) = obc%ry_oblique(i,j,k)
1867  segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
1868  enddo
1869  enddo
1870  endif
1871  enddo
1872  endif
1873 
1874  ! Now tracers (if any)
1875  do n=1,obc%number_of_segments
1876  segment=>obc%segment(n)
1877  if (associated(segment%tr_Reg)) then
1878  if (segment%is_E_or_W) then
1879  i = segment%HI%IsdB
1880  do m=1,obc%ntr
1881  if (associated(segment%tr_Reg%Tr(m)%tres)) then
1882  do k=1,g%ke
1883  do j=segment%HI%jsd,segment%HI%jed
1884  segment%tr_Reg%Tr(m)%tres(i,j,k) = obc%tres_x(i,j,k,m)
1885  enddo
1886  enddo
1887  endif
1888  enddo
1889  else
1890  j = segment%HI%JsdB
1891  do m=1,obc%ntr
1892  if (associated(segment%tr_Reg%Tr(m)%tres)) then
1893  do k=1,g%ke
1894  do i=segment%HI%isd,segment%HI%ied
1895  segment%tr_Reg%Tr(m)%tres(i,j,k) = obc%tres_y(i,j,k,m)
1896  enddo
1897  enddo
1898  endif
1899  enddo
1900  endif
1901  endif
1902  enddo
1903 
1904  gamma_u = obc%gamma_uv
1905  rx_max = obc%rx_max ; ry_max = obc%rx_max
1906  do n=1,obc%number_of_segments
1907  segment=>obc%segment(n)
1908  if (.not. segment%on_pe) cycle
1909  if (segment%oblique) call gradient_at_q_points(g, segment, u_new(:,:,:), v_new(:,:,:))
1910  if (segment%direction == obc_direction_e) then
1911  i=segment%HI%IsdB
1912  if (i<g%HI%IscB) cycle
1913  do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
1914  if (segment%radiation) then
1915  dhdt = (u_old(i-1,j,k) - u_new(i-1,j,k)) !old-new
1916  dhdx = (u_new(i-1,j,k) - u_new(i-2,j,k)) !in new time backward sasha for I-1
1917  rx_new = 0.0
1918  if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max) ! outward phase speed
1919  if (gamma_u < 1.0) then
1920  rx_avg = (1.0-gamma_u)*segment%rx_norm_rad(i,j,k) + gamma_u*rx_new
1921  else
1922  rx_avg = rx_new
1923  endif
1924  segment%rx_norm_rad(i,j,k) = rx_avg
1925  ! The new boundary value is interpolated between future interior
1926  ! value, u_new(I-1) and past boundary value but with barotropic
1927  ! accelerations, u_new(I).
1928  segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) / (1.0+rx_avg)
1929  ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
1930  ! implemented as a work-around to limitations in restart capability
1931  if (gamma_u < 1.0) then
1932  obc%rx_normal(i,j,k) = segment%rx_norm_rad(i,j,k)
1933  endif
1934  elseif (segment%oblique) then
1935  dhdt = (u_old(i-1,j,k) - u_new(i-1,j,k)) !old-new
1936  dhdx = (u_new(i-1,j,k) - u_new(i-2,j,k)) !in new time backward sasha for I-1
1937  if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0) then
1938  dhdy = segment%grad_normal(j-1,1,k)
1939  elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0) then
1940  dhdy = 0.0
1941  else
1942  dhdy = segment%grad_normal(j,1,k)
1943  endif
1944  if (dhdt*dhdx < 0.0) dhdt = 0.0
1945  rx_new = dhdt*dhdx
1946  cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
1947  ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
1948  if (gamma_u < 1.0) then
1949  rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
1950  ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
1951  cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
1952  else
1953  rx_avg = rx_new
1954  ry_avg = ry_new
1955  cff_avg = cff_new
1956  endif
1957  segment%rx_norm_obl(i,j,k) = rx_avg
1958  segment%ry_norm_obl(i,j,k) = ry_avg
1959  segment%cff_normal(i,j,k) = cff_avg
1960  segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) - &
1961  (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + &
1962  min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
1963  (cff_avg + rx_avg)
1964  if (gamma_u < 1.0) then
1965  ! Copy restart fields into 3-d arrays. This is an inefficient and temporary
1966  ! implementation as a work-around to limitations in restart capability
1967  obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
1968  obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
1969  obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
1970  endif
1971  elseif (segment%gradient) then
1972  segment%normal_vel(i,j,k) = u_new(i-1,j,k)
1973  endif
1974  if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then
1975  ! dhdt gets set to 0 on inflow in oblique case
1976  if (dhdt*dhdx <= 0.0) then
1977  tau = segment%Velocity_nudging_timescale_in
1978  else
1979  tau = segment%Velocity_nudging_timescale_out
1980  endif
1981  gamma_2 = dt / (tau + dt)
1982  segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
1983  gamma_2 * segment%nudged_normal_vel(i,j,k)
1984  endif
1985  enddo ; enddo
1986  if (segment%radiation_tan .or. segment%radiation_grad) then
1987  i=segment%HI%IsdB
1988  allocate(rx_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1989  do k=1,nz
1990  if (gamma_u < 1.0) then
1991  rx_tang_rad(i,segment%HI%JsdB,k) = segment%rx_norm_rad(i,segment%HI%jsd,k)
1992  rx_tang_rad(i,segment%HI%JedB,k) = segment%rx_norm_rad(i,segment%HI%jed,k)
1993  do j=segment%HI%JsdB+1,segment%HI%JedB-1
1994  rx_tang_rad(i,j,k) = 0.5*(segment%rx_norm_rad(i,j,k) + segment%rx_norm_rad(i,j+1,k))
1995  enddo
1996  else
1997  do j=segment%HI%JsdB,segment%HI%JedB
1998  dhdt = v_old(i,j,k)-v_new(i,j,k) !old-new
1999  dhdx = v_new(i,j,k)-v_new(i-1,j,k) !in new time backward sasha for I-1
2000  rx_tang_rad(i,j,k) = 0.0
2001  if (dhdt*dhdx > 0.0) rx_tang_rad(i,j,k) = min( (dhdt/dhdx), rx_max) ! outward phase speed
2002  enddo
2003  endif
2004  enddo
2005  if (segment%radiation_tan) then
2006  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2007  rx_avg = rx_tang_rad(i,j,k)
2008  segment%tangential_vel(i,j,k) = (v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) / (1.0+rx_avg)
2009  enddo ; enddo
2010  endif
2011  if (segment%nudged_tan) then
2012  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2013  ! dhdt gets set to 0 on inflow in oblique case
2014  if (rx_tang_rad(i,j,k) <= 0.0) then
2015  tau = segment%Velocity_nudging_timescale_in
2016  else
2017  tau = segment%Velocity_nudging_timescale_out
2018  endif
2019  gamma_2 = dt / (tau + dt)
2020  segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2021  gamma_2 * segment%nudged_tangential_vel(i,j,k)
2022  enddo ; enddo
2023  endif
2024  if (segment%radiation_grad) then
2025  js_obc = max(segment%HI%JsdB,g%jsd+1)
2026  je_obc = min(segment%HI%JedB,g%jed-1)
2027  do k=1,nz ; do j=js_obc,je_obc
2028  rx_avg = rx_tang_rad(i,j,k)
2029 ! if (G%mask2dCu(I-1,j) > 0.0 .and. G%mask2dCu(I-1,j+1) > 0.0) then
2030 ! rx_avg = 0.5*(u_new(I-1,j,k) + u_new(I-1,j+1,k)) * dt * G%IdxBu(I-1,J)
2031 ! elseif (G%mask2dCu(I-1,j) > 0.0) then
2032 ! rx_avg = u_new(I-1,j,k) * dt * G%IdxBu(I-1,J)
2033 ! elseif (G%mask2dCu(I-1,j+1) > 0.0) then
2034 ! rx_avg = u_new(I-1,j+1,k) * dt * G%IdxBu(I-1,J)
2035 ! else
2036 ! rx_avg = 0.0
2037 ! endif
2038  segment%tangential_grad(i,j,k) = ((v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
2039  rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) / (1.0+rx_avg)
2040  enddo ; enddo
2041  endif
2042  if (segment%nudged_grad) then
2043  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2044  ! dhdt gets set to 0 on inflow in oblique case
2045  if (rx_tang_rad(i,j,k) <= 0.0) then
2046  tau = segment%Velocity_nudging_timescale_in
2047  else
2048  tau = segment%Velocity_nudging_timescale_out
2049  endif
2050  gamma_2 = dt / (tau + dt)
2051  segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2052  gamma_2 * segment%nudged_tangential_grad(i,j,k)
2053  enddo ; enddo
2054  endif
2055  deallocate(rx_tang_rad)
2056  endif
2057  if (segment%oblique_tan .or. segment%oblique_grad) then
2058  i=segment%HI%IsdB
2059  allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2060  allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2061  allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2062  do k=1,nz
2063  if (gamma_u < 1.0) then
2064  rx_tang_obl(i,segment%HI%JsdB,k) = segment%rx_norm_obl(i,segment%HI%jsd,k)
2065  rx_tang_obl(i,segment%HI%JedB,k) = segment%rx_norm_obl(i,segment%HI%jed,k)
2066  ry_tang_obl(i,segment%HI%JsdB,k) = segment%ry_norm_obl(i,segment%HI%jsd,k)
2067  ry_tang_obl(i,segment%HI%JedB,k) = segment%ry_norm_obl(i,segment%HI%jed,k)
2068  cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
2069  cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
2070  do j=segment%HI%JsdB+1,segment%HI%JedB-1
2071  rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i,j+1,k))
2072  ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i,j+1,k))
2073  cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
2074  enddo
2075  else
2076  do j=segment%HI%JsdB,segment%HI%JedB
2077  dhdt = v_old(i,j,k)-v_new(i,j,k) !old-new
2078  dhdx = v_new(i,j,k)-v_new(i-1,j,k) !in new time backward sasha for I-1
2079  if (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) > 0.0) then
2080  dhdy = segment%grad_tan(j,1,k)
2081  elseif (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) == 0.0) then
2082  dhdy = 0.0
2083  else
2084  dhdy = segment%grad_tan(j+1,1,k)
2085  endif
2086  if (dhdt*dhdx < 0.0) dhdt = 0.0
2087  rx_new = dhdt*dhdx
2088  cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2089  ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2090  rx_tang_obl(i,j,k) = rx_new
2091  ry_tang_obl(i,j,k) = ry_new
2092  cff_tangential(i,j,k) = cff_new
2093  enddo
2094  endif
2095  enddo
2096  if (segment%oblique_tan) then
2097  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2098  rx_avg = rx_tang_obl(i,j,k)
2099  ry_avg = ry_tang_obl(i,j,k)
2100  cff_avg = cff_tangential(i,j,k)
2101  segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) - &
2102  (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + &
2103  min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
2104  (cff_avg + rx_avg)
2105  enddo ; enddo
2106  endif
2107  if (segment%nudged_tan) then
2108  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2109  ! dhdt gets set to 0 on inflow in oblique case
2110  if (rx_tang_obl(i,j,k) <= 0.0) then
2111  tau = segment%Velocity_nudging_timescale_in
2112  else
2113  tau = segment%Velocity_nudging_timescale_out
2114  endif
2115  gamma_2 = dt / (tau + dt)
2116  segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2117  gamma_2 * segment%nudged_tangential_vel(i,j,k)
2118  enddo ; enddo
2119  endif
2120  if (segment%oblique_grad) then
2121  js_obc = max(segment%HI%JsdB,g%jsd+1)
2122  je_obc = min(segment%HI%JedB,g%jed-1)
2123  do k=1,nz ; do j=segment%HI%JsdB+1,segment%HI%JedB-1
2124  rx_avg = rx_tang_obl(i,j,k)
2125  ry_avg = ry_tang_obl(i,j,k)
2126  cff_avg = cff_tangential(i,j,k)
2127  segment%tangential_grad(i,j,k) = &
2128  ((cff_avg*(v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
2129  rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) - &
2130  (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + &
2131  min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k)) ) / &
2132  (cff_avg + rx_avg)
2133  enddo ; enddo
2134  endif
2135  if (segment%nudged_grad) then
2136  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2137  ! dhdt gets set to 0 on inflow in oblique case
2138  if (rx_tang_obl(i,j,k) <= 0.0) then
2139  tau = segment%Velocity_nudging_timescale_in
2140  else
2141  tau = segment%Velocity_nudging_timescale_out
2142  endif
2143  gamma_2 = dt / (tau + dt)
2144  segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2145  gamma_2 * segment%nudged_tangential_grad(i,j,k)
2146  enddo ; enddo
2147  endif
2148  deallocate(rx_tang_obl)
2149  deallocate(ry_tang_obl)
2150  deallocate(cff_tangential)
2151  endif
2152  endif
2153 
2154  if (segment%direction == obc_direction_w) then
2155  i=segment%HI%IsdB
2156  if (i>g%HI%IecB) cycle
2157  do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
2158  if (segment%radiation) then
2159  dhdt = (u_old(i+1,j,k) - u_new(i+1,j,k)) !old-new
2160  dhdx = (u_new(i+1,j,k) - u_new(i+2,j,k)) !in new time forward sasha for I+1
2161  rx_new = 0.0
2162  if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
2163  if (gamma_u < 1.0) then
2164  rx_avg = (1.0-gamma_u)*segment%rx_norm_rad(i,j,k) + gamma_u*rx_new
2165  else
2166  rx_avg = rx_new
2167  endif
2168  segment%rx_norm_rad(i,j,k) = rx_avg
2169  ! The new boundary value is interpolated between future interior
2170  ! value, u_new(I+1) and past boundary value but with barotropic
2171  ! accelerations, u_new(I).
2172  segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) / (1.0+rx_avg)
2173  if (gamma_u < 1.0) then
2174  ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
2175  ! implemented as a work-around to limitations in restart capability
2176  obc%rx_normal(i,j,k) = segment%rx_norm_rad(i,j,k)
2177  endif
2178  elseif (segment%oblique) then
2179  dhdt = (u_old(i+1,j,k) - u_new(i+1,j,k)) !old-new
2180  dhdx = (u_new(i+1,j,k) - u_new(i+2,j,k)) !in new time forward sasha for I+1
2181  if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0) then
2182  dhdy = segment%grad_normal(j-1,1,k)
2183  elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0) then
2184  dhdy = 0.0
2185  else
2186  dhdy = segment%grad_normal(j,1,k)
2187  endif
2188  if (dhdt*dhdx < 0.0) dhdt = 0.0
2189 
2190  rx_new = dhdt*dhdx
2191  cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2192  ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2193  if (gamma_u < 1.0) then
2194  rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2195  ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2196  cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2197  else
2198  rx_avg = rx_new
2199  ry_avg = ry_new
2200  cff_avg = cff_new
2201  endif
2202  segment%rx_norm_obl(i,j,k) = rx_avg
2203  segment%ry_norm_obl(i,j,k) = ry_avg
2204  segment%cff_normal(i,j,k) = cff_avg
2205  segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) - &
2206  (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + &
2207  min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
2208  (cff_avg + rx_avg)
2209  if (gamma_u < 1.0) then
2210  ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
2211  ! implemented as a work-around to limitations in restart capability
2212  obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2213  obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2214  obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2215  endif
2216  elseif (segment%gradient) then
2217  segment%normal_vel(i,j,k) = u_new(i+1,j,k)
2218  endif
2219  if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then
2220  ! dhdt gets set to 0. on inflow in oblique case
2221  if (dhdt*dhdx <= 0.0) then
2222  tau = segment%Velocity_nudging_timescale_in
2223  else
2224  tau = segment%Velocity_nudging_timescale_out
2225  endif
2226  gamma_2 = dt / (tau + dt)
2227  segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2228  gamma_2 * segment%nudged_normal_vel(i,j,k)
2229  endif
2230  enddo ; enddo
2231  if (segment%radiation_tan .or. segment%radiation_grad) then
2232  i=segment%HI%IsdB
2233  allocate(rx_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2234  do k=1,nz
2235  if (gamma_u < 1.0) then
2236  rx_tang_rad(i,segment%HI%JsdB,k) = segment%rx_norm_rad(i,segment%HI%jsd,k)
2237  rx_tang_rad(i,segment%HI%JedB,k) = segment%rx_norm_rad(i,segment%HI%jed,k)
2238  do j=segment%HI%JsdB+1,segment%HI%JedB-1
2239  rx_tang_rad(i,j,k) = 0.5*(segment%rx_norm_rad(i,j,k) + segment%rx_norm_rad(i,j+1,k))
2240  enddo
2241  else
2242  do j=segment%HI%JsdB,segment%HI%JedB
2243  dhdt = v_old(i+1,j,k)-v_new(i+1,j,k) !old-new
2244  dhdx = v_new(i+1,j,k)-v_new(i+2,j,k) !in new time backward sasha for I-1
2245  rx_tang_rad(i,j,k) = 0.0
2246  if (dhdt*dhdx > 0.0) rx_tang_rad(i,j,k) = min( (dhdt/dhdx), rx_max) ! outward phase speed
2247  enddo
2248  endif
2249  enddo
2250  if (segment%radiation_tan) then
2251  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2252  rx_avg = rx_tang_rad(i,j,k)
2253  segment%tangential_vel(i,j,k) = (v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) / (1.0+rx_avg)
2254  enddo ; enddo
2255  endif
2256  if (segment%nudged_tan) then
2257  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2258  ! dhdt gets set to 0 on inflow in oblique case
2259  if (rx_tang_rad(i,j,k) <= 0.0) then
2260  tau = segment%Velocity_nudging_timescale_in
2261  else
2262  tau = segment%Velocity_nudging_timescale_out
2263  endif
2264  gamma_2 = dt / (tau + dt)
2265  segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2266  gamma_2 * segment%nudged_tangential_vel(i,j,k)
2267  enddo ; enddo
2268  endif
2269  if (segment%radiation_grad) then
2270  js_obc = max(segment%HI%JsdB,g%jsd+1)
2271  je_obc = min(segment%HI%JedB,g%jed-1)
2272  do k=1,nz ; do j=js_obc,je_obc
2273  rx_avg = rx_tang_rad(i,j,k)
2274 ! if (G%mask2dCu(I+1,j) > 0.0 .and. G%mask2dCu(I+1,j+1) > 0.0) then
2275 ! rx_avg = 0.5*(u_new(I+1,j,k) + u_new(I+1,j+1,k)) * dt * G%IdxBu(I+1,J)
2276 ! elseif (G%mask2dCu(I+1,j) > 0.0) then
2277 ! rx_avg = u_new(I+1,j,k) * dt * G%IdxBu(I+1,J)
2278 ! elseif (G%mask2dCu(I+1,j+1) > 0.0) then
2279 ! rx_avg = u_new(I+1,j+1,k) * dt * G%IdxBu(I+1,J)
2280 ! else
2281 ! rx_avg = 0.0
2282 ! endif
2283  segment%tangential_grad(i,j,k) = ((v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
2284  rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) / (1.0+rx_avg)
2285  enddo ; enddo
2286  endif
2287  if (segment%nudged_grad) then
2288  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2289  ! dhdt gets set to 0 on inflow in oblique case
2290  if (rx_tang_rad(i,j,k) <= 0.0) then
2291  tau = segment%Velocity_nudging_timescale_in
2292  else
2293  tau = segment%Velocity_nudging_timescale_out
2294  endif
2295  gamma_2 = dt / (tau + dt)
2296  segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2297  gamma_2 * segment%nudged_tangential_grad(i,j,k)
2298  enddo ; enddo
2299  endif
2300  deallocate(rx_tang_rad)
2301  endif
2302  if (segment%oblique_tan .or. segment%oblique_grad) then
2303  i=segment%HI%IsdB
2304  allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2305  allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2306  allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2307  do k=1,nz
2308  if (gamma_u < 1.0) then
2309  rx_tang_obl(i,segment%HI%JsdB,k) = segment%rx_norm_obl(i,segment%HI%jsd,k)
2310  rx_tang_obl(i,segment%HI%JedB,k) = segment%rx_norm_obl(i,segment%HI%jed,k)
2311  ry_tang_obl(i,segment%HI%JsdB,k) = segment%ry_norm_obl(i,segment%HI%jsd,k)
2312  ry_tang_obl(i,segment%HI%JedB,k) = segment%ry_norm_obl(i,segment%HI%jed,k)
2313  cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
2314  cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
2315  do j=segment%HI%JsdB+1,segment%HI%JedB-1
2316  rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i,j+1,k))
2317  ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i,j+1,k))
2318  cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
2319  enddo
2320  else
2321  do j=segment%HI%JsdB,segment%HI%JedB
2322  dhdt = v_old(i+1,j,k)-v_new(i+1,j,k) !old-new
2323  dhdx = v_new(i+1,j,k)-v_new(i+2,j,k) !in new time backward sasha for I-1
2324  if (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) > 0.0) then
2325  dhdy = segment%grad_tan(j,1,k)
2326  elseif (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) == 0.0) then
2327  dhdy = 0.0
2328  else
2329  dhdy = segment%grad_tan(j+1,1,k)
2330  endif
2331  if (dhdt*dhdx < 0.0) dhdt = 0.0
2332  rx_new = dhdt*dhdx
2333  cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2334  ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2335  rx_tang_obl(i,j,k) = rx_new
2336  ry_tang_obl(i,j,k) = ry_new
2337  cff_tangential(i,j,k) = cff_new
2338  enddo
2339  endif
2340  enddo
2341  if (segment%oblique_tan) then
2342  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2343  rx_avg = rx_tang_obl(i,j,k)
2344  ry_avg = ry_tang_obl(i,j,k)
2345  cff_avg = cff_tangential(i,j,k)
2346  segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) - &
2347  (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + &
2348  min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
2349  (cff_avg + rx_avg)
2350  enddo ; enddo
2351  endif
2352  if (segment%nudged_tan) then
2353  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2354  ! dhdt gets set to 0 on inflow in oblique case
2355  if (rx_tang_obl(i,j,k) <= 0.0) then
2356  tau = segment%Velocity_nudging_timescale_in
2357  else
2358  tau = segment%Velocity_nudging_timescale_out
2359  endif
2360  gamma_2 = dt / (tau + dt)
2361  segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2362  gamma_2 * segment%nudged_tangential_vel(i,j,k)
2363  enddo ; enddo
2364  endif
2365  if (segment%oblique_grad) then
2366  js_obc = max(segment%HI%JsdB,g%jsd+1)
2367  je_obc = min(segment%HI%JedB,g%jed-1)
2368  do k=1,nz ; do j=segment%HI%JsdB+1,segment%HI%JedB-1
2369  rx_avg = rx_tang_obl(i,j,k)
2370  ry_avg = ry_tang_obl(i,j,k)
2371  cff_avg = cff_tangential(i,j,k)
2372  segment%tangential_grad(i,j,k) = &
2373  ((cff_avg*(v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
2374  rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) - &
2375  (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + &
2376  min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k))) / &
2377  (cff_avg + rx_avg)
2378  enddo ; enddo
2379  endif
2380  if (segment%nudged_grad) then
2381  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2382  ! dhdt gets set to 0 on inflow in oblique case
2383  if (rx_tang_obl(i,j,k) <= 0.0) then
2384  tau = segment%Velocity_nudging_timescale_in
2385  else
2386  tau = segment%Velocity_nudging_timescale_out
2387  endif
2388  gamma_2 = dt / (tau + dt)
2389  segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2390  gamma_2 * segment%nudged_tangential_grad(i,j,k)
2391  enddo ; enddo
2392  endif
2393  deallocate(rx_tang_obl)
2394  deallocate(ry_tang_obl)
2395  deallocate(cff_tangential)
2396  endif
2397  endif
2398 
2399  if (segment%direction == obc_direction_n) then
2400  j=segment%HI%JsdB
2401  if (j<g%HI%JscB) cycle
2402  do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
2403  if (segment%radiation) then
2404  dhdt = (v_old(i,j-1,k) - v_new(i,j-1,k)) !old-new
2405  dhdy = (v_new(i,j-1,k) - v_new(i,j-2,k)) !in new time backward sasha for J-1
2406  ry_new = 0.0
2407  if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2408  if (gamma_u < 1.0) then
2409  ry_avg = (1.0-gamma_u)*segment%ry_norm_rad(i,j,k) + gamma_u*ry_new
2410  else
2411  ry_avg = ry_new
2412  endif
2413  segment%ry_norm_rad(i,j,k) = ry_avg
2414  ! The new boundary value is interpolated between future interior
2415  ! value, v_new(J-1) and past boundary value but with barotropic
2416  ! accelerations, v_new(J).
2417  segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) / (1.0+ry_avg)
2418  if (gamma_u < 1.0) then
2419  ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
2420  ! implemented as a work-around to limitations in restart capability
2421  obc%ry_normal(i,j,k) = segment%ry_norm_rad(i,j,k)
2422  endif
2423  elseif (segment%oblique) then
2424  dhdt = (v_old(i,j-1,k) - v_new(i,j-1,k)) !old-new
2425  dhdy = (v_new(i,j-1,k) - v_new(i,j-2,k)) !in new time backward sasha for J-1
2426  if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0) then
2427  dhdx = segment%grad_normal(i-1,1,k)
2428  elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0) then
2429  dhdx = 0.0
2430  else
2431  dhdx = segment%grad_normal(i,1,k)
2432  endif
2433  if (dhdt*dhdy < 0.0) dhdt = 0.0
2434  ry_new = dhdt*dhdy
2435  cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2436  rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
2437  if (gamma_u < 1.0) then
2438  rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2439  ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2440  cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2441  else
2442  rx_avg = rx_new
2443  ry_avg = ry_new
2444  cff_avg = cff_new
2445  endif
2446  segment%rx_norm_obl(i,j,k) = rx_avg
2447  segment%ry_norm_obl(i,j,k) = ry_avg
2448  segment%cff_normal(i,j,k) = cff_avg
2449  segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) - &
2450  (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) +&
2451  min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
2452  (cff_avg + ry_avg)
2453  if (gamma_u < 1.0) then
2454  ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
2455  ! implemented as a work-around to limitations in restart capability
2456  obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2457  obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2458  obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2459  endif
2460  elseif (segment%gradient) then
2461  segment%normal_vel(i,j,k) = v_new(i,j-1,k)
2462  endif
2463  if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then
2464  ! dhdt gets set to 0 on inflow in oblique case
2465  if (dhdt*dhdy <= 0.0) then
2466  tau = segment%Velocity_nudging_timescale_in
2467  else
2468  tau = segment%Velocity_nudging_timescale_out
2469  endif
2470  gamma_2 = dt / (tau + dt)
2471  segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2472  gamma_2 * segment%nudged_normal_vel(i,j,k)
2473  endif
2474  enddo ; enddo
2475  if (segment%radiation_tan .or. segment%radiation_grad) then
2476  j=segment%HI%JsdB
2477  allocate(ry_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2478  do k=1,nz
2479  if (gamma_u < 1.0) then
2480  ry_tang_rad(segment%HI%IsdB,j,k) = segment%ry_norm_rad(segment%HI%isd,j,k)
2481  ry_tang_rad(segment%HI%IedB,j,k) = segment%ry_norm_rad(segment%HI%ied,j,k)
2482  do i=segment%HI%IsdB+1,segment%HI%IedB-1
2483  ry_tang_rad(i,j,k) = 0.5*(segment%ry_norm_rad(i,j,k) + segment%ry_norm_rad(i+1,j,k))
2484  enddo
2485  else
2486  do i=segment%HI%IsdB,segment%HI%IedB
2487  dhdt = u_old(i,j-1,k)-u_new(i,j-1,k) !old-new
2488  dhdy = u_new(i,j-1,k)-u_new(i,j-2,k) !in new time backward sasha for I-1
2489  ry_tang_rad(i,j,k) = 0.0
2490  if (dhdt*dhdy > 0.0) ry_tang_rad(i,j,k) = min( (dhdt/dhdy), rx_max) ! outward phase speed
2491  enddo
2492  endif
2493  enddo
2494  if (segment%radiation_tan) then
2495  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2496  ry_avg = ry_tang_rad(i,j,k)
2497  segment%tangential_vel(i,j,k) = (u_new(i,j,k) + ry_avg*u_new(i,j-1,k)) / (1.0+ry_avg)
2498  enddo ; enddo
2499  endif
2500  if (segment%nudged_tan) then
2501  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2502  ! dhdt gets set to 0 on inflow in oblique case
2503  if (ry_tang_rad(i,j,k) <= 0.0) then
2504  tau = segment%Velocity_nudging_timescale_in
2505  else
2506  tau = segment%Velocity_nudging_timescale_out
2507  endif
2508  gamma_2 = dt / (tau + dt)
2509  segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2510  gamma_2 * segment%nudged_tangential_vel(i,j,k)
2511  enddo ; enddo
2512  endif
2513  if (segment%radiation_grad) then
2514  is_obc = max(segment%HI%IsdB,g%isd+1)
2515  ie_obc = min(segment%HI%IedB,g%ied-1)
2516  do k=1,nz ; do i=is_obc,ie_obc
2517  ry_avg = ry_tang_rad(i,j,k)
2518 ! if (G%mask2dCv(i,J-1) > 0.0 .and. G%mask2dCv(i+1,J-1) > 0.0) then
2519 ! ry_avg = 0.5*(v_new(i,J-1,k) + v_new(i+1,J-1,k) * dt * G%IdyBu(I,J-1))
2520 ! elseif (G%mask2dCv(i,J-1) > 0.0) then
2521 ! ry_avg = v_new(i,J-1,k) * dt *G%IdyBu(I,J-1)
2522 ! elseif (G%mask2dCv(i+1,J-1) > 0.0) then
2523 ! ry_avg = v_new(i+1,J-1,k) * dt *G%IdyBu(I,J-1)
2524 ! else
2525 ! ry_avg = 0.0
2526 ! endif
2527  segment%tangential_grad(i,j,k) = ((u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
2528  ry_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) / (1.0+ry_avg)
2529  enddo ; enddo
2530  endif
2531  if (segment%nudged_grad) then
2532  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2533  ! dhdt gets set to 0 on inflow in oblique case
2534  if (ry_tang_rad(i,j,k) <= 0.0) then
2535  tau = segment%Velocity_nudging_timescale_in
2536  else
2537  tau = segment%Velocity_nudging_timescale_out
2538  endif
2539  gamma_2 = dt / (tau + dt)
2540  segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2541  gamma_2 * segment%nudged_tangential_grad(i,j,k)
2542  enddo ; enddo
2543  endif
2544  deallocate(ry_tang_rad)
2545  endif
2546  if (segment%oblique_tan .or. segment%oblique_grad) then
2547  j=segment%HI%JsdB
2548  allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2549  allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2550  allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2551  do k=1,nz
2552  if (gamma_u < 1.0) then
2553  rx_tang_obl(segment%HI%IsdB,j,k) = segment%rx_norm_obl(segment%HI%isd,j,k)
2554  rx_tang_obl(segment%HI%IedB,j,k) = segment%rx_norm_obl(segment%HI%ied,j,k)
2555  ry_tang_obl(segment%HI%IsdB,j,k) = segment%ry_norm_obl(segment%HI%isd,j,k)
2556  ry_tang_obl(segment%HI%IedB,j,k) = segment%ry_norm_obl(segment%HI%ied,j,k)
2557  cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
2558  cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
2559  do i=segment%HI%IsdB+1,segment%HI%IedB-1
2560  rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i+1,j,k))
2561  ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i+1,j,k))
2562  cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
2563  enddo
2564  else
2565  do i=segment%HI%IsdB,segment%HI%IedB
2566  dhdt = u_old(i,j,k)-u_new(i,j,k) !old-new
2567  dhdy = u_new(i,j,k)-u_new(i,j-1,k) !in new time backward sasha for I-1
2568  if (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) > 0.0) then
2569  dhdx = segment%grad_tan(i,1,k)
2570  elseif (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) == 0.0) then
2571  dhdx = 0.0
2572  else
2573  dhdx = segment%grad_tan(i+1,1,k)
2574  endif
2575  if (dhdt*dhdy < 0.0) dhdt = 0.0
2576  ry_new = dhdt*dhdy
2577  cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2578  rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
2579  rx_tang_obl(i,j,k) = rx_new
2580  ry_tang_obl(i,j,k) = ry_new
2581  cff_tangential(i,j,k) = cff_new
2582  enddo
2583  endif
2584  enddo
2585  if (segment%oblique_tan) then
2586  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2587  rx_avg = rx_tang_obl(i,j,k)
2588  ry_avg = ry_tang_obl(i,j,k)
2589  cff_avg = cff_tangential(i,j,k)
2590  segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + ry_avg*u_new(i,j-1,k)) - &
2591  (max(rx_avg,0.0)*segment%grad_tan(i,2,k) + &
2592  min(rx_avg,0.0)*segment%grad_tan(i+1,2,k))) / &
2593  (cff_avg + ry_avg)
2594  enddo ; enddo
2595  endif
2596  if (segment%nudged_tan) then
2597  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2598  ! dhdt gets set to 0 on inflow in oblique case
2599  if (ry_tang_obl(i,j,k) <= 0.0) then
2600  tau = segment%Velocity_nudging_timescale_in
2601  else
2602  tau = segment%Velocity_nudging_timescale_out
2603  endif
2604  gamma_2 = dt / (tau + dt)
2605  segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2606  gamma_2 * segment%nudged_tangential_vel(i,j,k)
2607  enddo ; enddo
2608  endif
2609  if (segment%oblique_grad) then
2610  is_obc = max(segment%HI%IsdB,g%isd+1)
2611  ie_obc = min(segment%HI%IedB,g%ied-1)
2612  do k=1,nz ; do i=segment%HI%IsdB+1,segment%HI%IedB-1
2613  rx_avg = rx_tang_obl(i,j,k)
2614  ry_avg = ry_tang_obl(i,j,k)
2615  cff_avg = cff_tangential(i,j,k)
2616  segment%tangential_grad(i,j,k) = &
2617  ((cff_avg*(u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
2618  ry_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) - &
2619  (max(rx_avg,0.0)*segment%grad_gradient(i,2,k) + &
2620  min(rx_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
2621  (cff_avg + ry_avg)
2622  enddo ; enddo
2623  endif
2624  if (segment%nudged_grad) then
2625  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2626  ! dhdt gets set to 0 on inflow in oblique case
2627  if (ry_tang_obl(i,j,k) <= 0.0) then
2628  tau = segment%Velocity_nudging_timescale_in
2629  else
2630  tau = segment%Velocity_nudging_timescale_out
2631  endif
2632  gamma_2 = dt / (tau + dt)
2633  segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2634  gamma_2 * segment%nudged_tangential_grad(i,j,k)
2635  enddo ; enddo
2636  endif
2637  deallocate(rx_tang_obl)
2638  deallocate(ry_tang_obl)
2639  deallocate(cff_tangential)
2640  endif
2641  endif
2642 
2643  if (segment%direction == obc_direction_s) then
2644  j=segment%HI%JsdB
2645  if (j>g%HI%JecB) cycle
2646  do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
2647  if (segment%radiation) then
2648  dhdt = (v_old(i,j+1,k) - v_new(i,j+1,k)) !old-new
2649  dhdy = (v_new(i,j+1,k) - v_new(i,j+2,k)) !in new time backward sasha for J-1
2650  ry_new = 0.0
2651  if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2652  if (gamma_u < 1.0) then
2653  ry_avg = (1.0-gamma_u)*segment%ry_norm_rad(i,j,k) + gamma_u*ry_new
2654  else
2655  ry_avg = ry_new
2656  endif
2657  segment%ry_norm_rad(i,j,k) = ry_avg
2658  ! The new boundary value is interpolated between future interior
2659  ! value, v_new(J+1) and past boundary value but with barotropic
2660  ! accelerations, v_new(J).
2661  segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) / (1.0+ry_avg)
2662  if (gamma_u < 1.0) then
2663  ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
2664  ! implemented as a work-around to limitations in restart capability
2665  obc%ry_normal(i,j,k) = segment%ry_norm_rad(i,j,k)
2666  endif
2667  elseif (segment%oblique) then
2668  dhdt = (v_old(i,j+1,k) - v_new(i,j+1,k)) !old-new
2669  dhdy = (v_new(i,j+1,k) - v_new(i,j+2,k)) !in new time backward sasha for J-1
2670  if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0) then
2671  dhdx = segment%grad_normal(i-1,1,k)
2672  elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0) then
2673  dhdx = 0.0
2674  else
2675  dhdx = segment%grad_normal(i,1,k)
2676  endif
2677  if (dhdt*dhdy < 0.0) dhdt = 0.0
2678 
2679  ry_new = dhdt*dhdy
2680  cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2681  rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
2682  if (gamma_u < 1.0) then
2683  rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2684  ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2685  cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2686  else
2687  rx_avg = rx_new
2688  ry_avg = ry_new
2689  cff_avg = cff_new
2690  endif
2691  segment%rx_norm_obl(i,j,k) = rx_avg
2692  segment%ry_norm_obl(i,j,k) = ry_avg
2693  segment%cff_normal(i,j,k) = cff_avg
2694  segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) - &
2695  (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) + &
2696  min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
2697  (cff_avg + ry_avg)
2698  if (gamma_u < 1.0) then
2699  ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues
2700  ! implemented as a work-around to limitations in restart capability
2701  obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2702  obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2703  obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2704  endif
2705  elseif (segment%gradient) then
2706  segment%normal_vel(i,j,k) = v_new(i,j+1,k)
2707  endif
2708  if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then
2709  ! dhdt gets set to 0 on inflow in oblique case
2710  if (dhdt*dhdy <= 0.0) then
2711  tau = segment%Velocity_nudging_timescale_in
2712  else
2713  tau = segment%Velocity_nudging_timescale_out
2714  endif
2715  gamma_2 = dt / (tau + dt)
2716  segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2717  gamma_2 * segment%nudged_normal_vel(i,j,k)
2718  endif
2719  enddo ; enddo
2720  if (segment%radiation_tan .or. segment%radiation_grad) then
2721  j=segment%HI%JsdB
2722  allocate(ry_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2723  do k=1,nz
2724  if (gamma_u < 1.0) then
2725  ry_tang_rad(segment%HI%IsdB,j,k) = segment%ry_norm_rad(segment%HI%isd,j,k)
2726  ry_tang_rad(segment%HI%IedB,j,k) = segment%ry_norm_rad(segment%HI%ied,j,k)
2727  do i=segment%HI%IsdB+1,segment%HI%IedB-1
2728  ry_tang_rad(i,j,k) = 0.5*(segment%ry_norm_rad(i,j,k) + segment%ry_norm_rad(i+1,j,k))
2729  enddo
2730  else
2731  do i=segment%HI%IsdB,segment%HI%IedB
2732  dhdt = u_old(i,j+1,k)-u_new(i,j+1,k) !old-new
2733  dhdy = u_new(i,j+1,k)-u_new(i,j+2,k) !in new time backward sasha for I-1
2734  ry_tang_rad(i,j,k) = 0.0
2735  if (dhdt*dhdy > 0.0) ry_tang_rad(i,j,k) = min( (dhdt/dhdy), rx_max) ! outward phase speed
2736  enddo
2737  endif
2738  enddo
2739  if (segment%radiation_tan) then
2740  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2741  ry_avg = ry_tang_rad(i,j,k)
2742  segment%tangential_vel(i,j,k) = (u_new(i,j+1,k) + ry_avg*u_new(i,j+2,k)) / (1.0+ry_avg)
2743  enddo ; enddo
2744  endif
2745  if (segment%nudged_tan) then
2746  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2747  ! dhdt gets set to 0 on inflow in oblique case
2748  if (ry_tang_rad(i,j,k) <= 0.0) then
2749  tau = segment%Velocity_nudging_timescale_in
2750  else
2751  tau = segment%Velocity_nudging_timescale_out
2752  endif
2753  gamma_2 = dt / (tau + dt)
2754  segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2755  gamma_2 * segment%nudged_tangential_vel(i,j,k)
2756  enddo ; enddo
2757  endif
2758  if (segment%radiation_grad) then
2759  is_obc = max(segment%HI%IsdB,g%isd+1)
2760  ie_obc = min(segment%HI%IedB,g%ied-1)
2761  do k=1,nz ; do i=is_obc,ie_obc
2762  ry_avg = ry_tang_rad(i,j,k)
2763 ! if (G%mask2dCv(i,J+1) > 0.0 .and. G%mask2dCv(i+1,J+1) > 0.0) then
2764 ! ry_avg = 0.5*(v_new(i,J+1,k) + v_new(i+1,J+1,k)) * dt * G%IdyBu(I,J+1)
2765 ! elseif (G%mask2dCv(i,J+1) > 0.0) then
2766 ! ry_avg = v_new(i,J+1,k) * dt * G%IdyBu(I,J+1)
2767 ! elseif (G%mask2dCv(i+1,J+1) > 0.0) then
2768 ! ry_avg = v_new(i+1,J+1,k) * dt * G%IdyBu(I,J+1)
2769 ! else
2770 ! ry_avg = 0.0
2771 ! endif
2772  segment%tangential_grad(i,j,k) = ((u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
2773  ry_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) / (1.0+ry_avg)
2774  enddo ; enddo
2775  endif
2776  if (segment%nudged_grad) then
2777  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2778  ! dhdt gets set to 0 on inflow in oblique case
2779  if (ry_tang_rad(i,j,k) <= 0.0) then
2780  tau = segment%Velocity_nudging_timescale_in
2781  else
2782  tau = segment%Velocity_nudging_timescale_out
2783  endif
2784  gamma_2 = dt / (tau + dt)
2785  segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2786  gamma_2 * segment%nudged_tangential_grad(i,j,k)
2787  enddo ; enddo
2788  endif
2789  deallocate(ry_tang_rad)
2790  endif
2791  if (segment%oblique_tan .or. segment%oblique_grad) then
2792  j=segment%HI%JsdB
2793  allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2794  allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2795  allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2796  do k=1,nz
2797  if (gamma_u < 1.0) then
2798  rx_tang_obl(segment%HI%IsdB,j,k) = segment%rx_norm_obl(segment%HI%isd,j,k)
2799  rx_tang_obl(segment%HI%IedB,j,k) = segment%rx_norm_obl(segment%HI%ied,j,k)
2800  ry_tang_obl(segment%HI%IsdB,j,k) = segment%ry_norm_obl(segment%HI%isd,j,k)
2801  ry_tang_obl(segment%HI%IedB,j,k) = segment%ry_norm_obl(segment%HI%ied,j,k)
2802  cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
2803  cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
2804  do i=segment%HI%IsdB+1,segment%HI%IedB-1
2805  rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i+1,j,k))
2806  ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i+1,j,k))
2807  cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
2808  enddo
2809  else
2810  do i=segment%HI%IsdB,segment%HI%IedB
2811  dhdt = u_old(i,j+1,k)-u_new(i,j+1,k) !old-new
2812  dhdy = u_new(i,j+1,k)-u_new(i,j+2,k) !in new time backward sasha for I-1
2813  if (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) > 0.0) then
2814  dhdx = segment%grad_tan(i,1,k)
2815  elseif (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) == 0.0) then
2816  dhdx = 0.0
2817  else
2818  dhdx = segment%grad_tan(i+1,1,k)
2819  endif
2820  if (dhdt*dhdy < 0.0) dhdt = 0.0
2821  ry_new = dhdt*dhdy
2822  cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2823  rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
2824  rx_tang_obl(i,j,k) = rx_new
2825  ry_tang_obl(i,j,k) = ry_new
2826  cff_tangential(i,j,k) = cff_new
2827  enddo
2828  endif
2829  enddo
2830  if (segment%oblique_tan) then
2831  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2832  rx_avg = rx_tang_obl(i,j,k)
2833  ry_avg = ry_tang_obl(i,j,k)
2834  cff_avg = cff_tangential(i,j,k)
2835  segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j+1,k) + ry_avg*u_new(i,j+2,k)) - &
2836  (max(rx_avg,0.0)*segment%grad_tan(i,2,k) + &
2837  min(rx_avg,0.0)*segment%grad_tan(i+1,2,k)) ) / &
2838  (cff_avg + ry_avg)
2839  enddo ; enddo
2840  endif
2841  if (segment%nudged_tan) then
2842  do k=1,nz ; do i=segment%HI%IsdB,segment%HI%IedB
2843  ! dhdt gets set to 0 on inflow in oblique case
2844  if (ry_tang_obl(i,j,k) <= 0.0) then
2845  tau = segment%Velocity_nudging_timescale_in
2846  else
2847  tau = segment%Velocity_nudging_timescale_out
2848  endif
2849  gamma_2 = dt / (tau + dt)
2850  segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2851  gamma_2 * segment%nudged_tangential_vel(i,j,k)
2852  enddo ; enddo
2853  endif
2854  if (segment%oblique_grad) then
2855  is_obc = max(segment%HI%IsdB,g%isd+1)
2856  ie_obc = min(segment%HI%IedB,g%ied-1)
2857  do k=1,nz ; do i=segment%HI%IsdB+1,segment%HI%IedB-1
2858  rx_avg = rx_tang_obl(i,j,k)
2859  ry_avg = ry_tang_obl(i,j,k)
2860  cff_avg = cff_tangential(i,j,k)
2861  segment%tangential_grad(i,j,k) = &
2862  ((cff_avg*(u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
2863  ry_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) - &
2864  (max(rx_avg,0.0)*segment%grad_gradient(i,2,k) + &
2865  min(rx_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
2866  (cff_avg + ry_avg)
2867  enddo ; enddo
2868  endif
2869  if (segment%nudged_grad) then
2870  do k=1,nz ; do j=segment%HI%JsdB,segment%HI%JedB
2871  ! dhdt gets set to 0 on inflow in oblique case
2872  if (ry_tang_obl(i,j,k) <= 0.0) then
2873  tau = segment%Velocity_nudging_timescale_in
2874  else
2875  tau = segment%Velocity_nudging_timescale_out
2876  endif
2877  gamma_2 = dt / (tau + dt)
2878  segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2879  gamma_2 * segment%nudged_tangential_grad(i,j,k)
2880  enddo ; enddo
2881  endif
2882  deallocate(rx_tang_obl)
2883  deallocate(ry_tang_obl)
2884  deallocate(cff_tangential)
2885  endif
2886  endif
2887  enddo
2888 
2889  ! Actually update u_new, v_new
2890  call open_boundary_apply_normal_flow(obc, g, u_new, v_new)
2891 
2892  call pass_vector(u_new, v_new, g%Domain, clock=id_clock_pass)
2893 

References gradient_at_q_points(), id_clock_pass, obc_direction_e, obc_direction_n, obc_direction_s, obc_direction_w, and open_boundary_apply_normal_flow().

Here is the call graph for this function:

◆ register_file_obc()

logical function, public mom_open_boundary::register_file_obc ( type(param_file_type), intent(in)  param_file,
type(file_obc_cs), pointer  CS,
type(obc_registry_type), pointer  OBC_Reg 
)

Add file to OBC registry.

Parameters
[in]param_fileparameter file.
csfile control structure.
obc_regOBC registry.

Definition at line 3951 of file MOM_open_boundary.F90.

3951  type(param_file_type), intent(in) :: param_file !< parameter file.
3952  type(file_OBC_CS), pointer :: CS !< file control structure.
3953  type(OBC_registry_type), pointer :: OBC_Reg !< OBC registry.
3954  logical :: register_file_OBC
3955  character(len=32) :: casename = "OBC file" !< This case's name.
3956 
3957  if (associated(cs)) then
3958  call mom_error(warning, "register_file_OBC called with an "// &
3959  "associated control structure.")
3960  return
3961  endif
3962  allocate(cs)
3963 
3964  ! Register the file for boundary updates.
3965  call register_obc(casename, param_file, obc_reg)
3966  register_file_obc = .true.
3967 

References register_obc().

Referenced by mom_boundary_update::call_obc_register().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ register_obc()

subroutine, public mom_open_boundary::register_obc ( character(len=32), intent(in)  name,
type(param_file_type), intent(in)  param_file,
type(obc_registry_type), pointer  Reg 
)

register open boundary objects for boundary updates.

Parameters
[in]nameOBC name used for error messages
[in]param_filefile to parse for model parameter values
regpointer to the tracer registry

Definition at line 3899 of file MOM_open_boundary.F90.

3899  character(len=32), intent(in) :: name !< OBC name used for error messages
3900  type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values
3901  type(OBC_registry_type), pointer :: Reg !< pointer to the tracer registry
3902  integer :: nobc
3903  character(len=256) :: mesg ! Message for error messages.
3904 
3905  if (.not. associated(reg)) call obc_registry_init(param_file, reg)
3906 
3907  if (reg%nobc>=max_fields_) then
3908  write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for &
3909  &all the open boundaries being registered via register_OBC.")') reg%nobc+1
3910  call mom_error(fatal,"MOM register_tracer: "//mesg)
3911  endif
3912  reg%nobc = reg%nobc + 1
3913  nobc = reg%nobc
3914 
3915  reg%OB(nobc)%name = name
3916 
3917  if (reg%locked) call mom_error(fatal, &
3918  "MOM register_tracer was called for variable "//trim(reg%OB(nobc)%name)//&
3919  " with a locked tracer registry.")
3920 

References obc_registry_init().

Referenced by dyed_channel_initialization::register_dyed_channel_obc(), and register_file_obc().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ register_segment_tracer()

subroutine, public mom_open_boundary::register_segment_tracer ( type(tracer_type), target  tr_ptr,
type(param_file_type), intent(in)  param_file,
type(verticalgrid_type), intent(in)  GV,
type(obc_segment_type), intent(inout)  segment,
real, intent(in), optional  OBC_scalar,
logical, intent(in), optional  OBC_array 
)
Parameters
[in]gvocean vertical grid structure
tr_ptrA 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.
[in]param_filefile to parse for model parameter values
[in,out]segmentcurrent segment data structure
[in]obc_scalarIf present, use scalar value for segment tracer inflow concentration.
[in]obc_arrayIf true, use array values for segment tracer inflow concentration.

Definition at line 4013 of file MOM_open_boundary.F90.

4013  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
4014  type(tracer_type), target :: tr_ptr !< A target that can be used to set a pointer to the
4015  !! stored value of tr. This target must be
4016  !! an enduring part of the control structure,
4017  !! because the tracer registry will use this memory,
4018  !! but it also means that any updates to this
4019  !! structure in the calling module will be
4020  !! available subsequently to the tracer registry.
4021  type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values
4022  type(OBC_segment_type), intent(inout) :: segment !< current segment data structure
4023  real, optional, intent(in) :: OBC_scalar !< If present, use scalar value for segment tracer
4024  !! inflow concentration.
4025  logical, optional, intent(in) :: OBC_array !< If true, use array values for segment tracer
4026  !! inflow concentration.
4027 
4028 
4029 ! Local variables
4030  integer :: ntseg
4031  integer :: isd, ied, jsd, jed
4032  integer :: IsdB, IedB, JsdB, JedB
4033  character(len=256) :: mesg ! Message for error messages.
4034 
4035  call segment_tracer_registry_init(param_file, segment)
4036 
4037  if (segment%tr_Reg%ntseg>=max_fields_) then
4038  write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for &
4039  &all the tracers being registered via register_tracer.")') segment%tr_Reg%ntseg+1
4040  call mom_error(fatal,"MOM register_tracer: "//mesg)
4041  endif
4042  segment%tr_Reg%ntseg = segment%tr_Reg%ntseg + 1
4043  ntseg = segment%tr_Reg%ntseg
4044 
4045  isd = segment%HI%isd ; ied = segment%HI%ied
4046  jsd = segment%HI%jsd ; jed = segment%HI%jed
4047  isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
4048  jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
4049 
4050  segment%tr_Reg%Tr(ntseg)%Tr => tr_ptr
4051  segment%tr_Reg%Tr(ntseg)%name = tr_ptr%name
4052 
4053  if (segment%tr_Reg%locked) call mom_error(fatal, &
4054  "MOM register_tracer was called for variable "//trim(segment%tr_Reg%Tr(ntseg)%name)//&
4055  " with a locked tracer registry.")
4056 
4057  if (present(obc_scalar)) segment%tr_Reg%Tr(ntseg)%OBC_inflow_conc = obc_scalar ! initialize tracer value later
4058  if (present(obc_array)) then
4059  if (segment%is_E_or_W) then
4060  allocate(segment%tr_Reg%Tr(ntseg)%t(isdb:iedb,jsd:jed,1:gv%ke));segment%tr_Reg%Tr(ntseg)%t(:,:,:)=0.0
4061  allocate(segment%tr_Reg%Tr(ntseg)%tres(isdb:iedb,jsd:jed,1:gv%ke));segment%tr_Reg%Tr(ntseg)%tres(:,:,:)=0.0
4062  segment%tr_Reg%Tr(ntseg)%is_initialized=.false.
4063  elseif (segment%is_N_or_S) then
4064  allocate(segment%tr_Reg%Tr(ntseg)%t(isd:ied,jsdb:jedb,1:gv%ke));segment%tr_Reg%Tr(ntseg)%t(:,:,:)=0.0
4065  allocate(segment%tr_Reg%Tr(ntseg)%tres(isd:ied,jsdb:jedb,1:gv%ke));segment%tr_Reg%Tr(ntseg)%tres(:,:,:)=0.0
4066  segment%tr_Reg%Tr(ntseg)%is_initialized=.false.
4067  endif
4068  endif
4069 

References segment_tracer_registry_init().

Referenced by dome_initialization::dome_set_obc_data(), dyed_obcs_initialization::dyed_obcs_set_obc_data(), and register_temp_salt_segments().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ register_temp_salt_segments()

subroutine, public mom_open_boundary::register_temp_salt_segments ( type(verticalgrid_type), intent(in)  GV,
type(ocean_obc_type), pointer  OBC,
type(tracer_registry_type), pointer  tr_Reg,
type(param_file_type), intent(in)  param_file 
)
Parameters
[in]gvocean vertical grid structure
obcOpen boundary structure
tr_regTracer registry
[in]param_filefile to parse for model parameter values

Definition at line 4088 of file MOM_open_boundary.F90.

4088  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
4089  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
4090  type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry
4091  type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values
4092 
4093 ! Local variables
4094  integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, nz, nf
4095  integer :: i, j, k, n
4096  character(len=32) :: name
4097  type(OBC_segment_type), pointer :: segment => null() ! pointer to segment type list
4098  type(tracer_type), pointer :: tr_ptr => null()
4099 
4100  if (.not. associated(obc)) return
4101 
4102  do n=1, obc%number_of_segments
4103  segment=>obc%segment(n)
4104  if (.not. segment%on_pe) cycle
4105 
4106  if (associated(segment%tr_Reg)) &
4107  call mom_error(fatal,"register_temp_salt_segments: tracer array was previously allocated")
4108 
4109  name = 'temp'
4110  call tracer_name_lookup(tr_reg, tr_ptr, name)
4111  call register_segment_tracer(tr_ptr, param_file, gv, segment, &
4112  obc_array=segment%temp_segment_data_exists)
4113  name = 'salt'
4114  call tracer_name_lookup(tr_reg, tr_ptr, name)
4115  call register_segment_tracer(tr_ptr, param_file, gv, segment, &
4116  obc_array=segment%salt_segment_data_exists)
4117  enddo
4118 

References register_segment_tracer().

Here is the call graph for this function:

◆ segment_tracer_registry_end()

subroutine, public mom_open_boundary::segment_tracer_registry_end ( type(segment_tracer_registry_type), pointer  Reg)

Clean up the segment tracer registry.

Parameters
regpointer to tracer registry

Definition at line 4074 of file MOM_open_boundary.F90.

4074  type(segment_tracer_registry_type), pointer :: Reg !< pointer to tracer registry
4075 
4076 ! Local variables
4077  integer n
4078 
4079  if (associated(reg)) then
4080  do n = 1, reg%ntseg
4081  if (associated(reg%Tr(n)%t)) deallocate(reg%Tr(n)%t)
4082  enddo
4083  deallocate(reg)
4084  endif

Referenced by deallocate_obc_segment_data().

Here is the caller graph for this function:

◆ segment_tracer_registry_init()

subroutine, public mom_open_boundary::segment_tracer_registry_init ( type(param_file_type), intent(in)  param_file,
type(obc_segment_type), intent(inout)  segment 
)

Initialize the segment tracer registry.

Parameters
[in]param_fileopen file to parse for model parameters
[in,out]segmentthe segment

Definition at line 3981 of file MOM_open_boundary.F90.

3981  type(param_file_type), intent(in) :: param_file !< open file to parse for model parameters
3982  type(OBC_segment_type), intent(inout) :: segment !< the segment
3983 
3984  integer, save :: init_calls = 0
3985 
3986 ! This include declares and sets the variable "version".
3987 #include "version_variable.h"
3988  character(len=40) :: mdl = "segment_tracer_registry_init" ! This routine's name.
3989  character(len=256) :: mesg ! Message for error messages.
3990 
3991  if (.not.associated(segment%tr_Reg)) then
3992  allocate(segment%tr_Reg)
3993  else
3994  return
3995  endif
3996 
3997  init_calls = init_calls + 1
3998 
3999  ! Read all relevant parameters and write them to the model log.
4000  if (init_calls == 1) call log_version(param_file, mdl, version, "")
4001 
4002 ! Need to call once per segment with tracers...
4003 ! if (init_calls > 1) then
4004 ! write(mesg,'("segment_tracer_registry_init called ",I3, &
4005 ! &" times with different registry pointers.")') init_calls
4006 ! if (is_root_pe()) call MOM_error(WARNING,"MOM_tracer"//mesg)
4007 ! endif
4008 

References mdl.

Referenced by register_segment_tracer().

Here is the caller graph for this function:

◆ set_tracer_data()

subroutine, public mom_open_boundary::set_tracer_data ( type(ocean_obc_type), pointer  OBC,
type(thermo_var_ptrs), intent(inout)  tv,
real, dimension(szi_(g),szj_(g), szk_(g)), intent(inout)  h,
type(ocean_grid_type), intent(inout)  G,
type(param_file_type), intent(in)  PF,
type(tracer_registry_type), pointer  tracer_Reg 
)

Sets the initial values of the tracer open boundary conditions. Redoing this elsewhere.

Parameters
[in,out]gOcean grid structure
obcOpen boundary structure
[in,out]tvThermodynamics structure
[in,out]hThickness
[in]pfParameter file handle
tracer_regTracer registry

Definition at line 3090 of file MOM_open_boundary.F90.

3090  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
3091  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
3092  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
3093  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h !< Thickness
3094  type(param_file_type), intent(in) :: PF !< Parameter file handle
3095  type(tracer_registry_type), pointer :: tracer_Reg !< Tracer registry
3096  ! Local variables
3097  integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, n
3098  integer :: isd_off, jsd_off
3099  integer :: IsdB, IedB, JsdB, JedB
3100  type(OBC_segment_type), pointer :: segment => null() ! pointer to segment type list
3101  character(len=40) :: mdl = "set_tracer_data" ! This subroutine's name.
3102  character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path
3103 
3104  real :: temp_u(G%domain%niglobal+1,G%domain%njglobal)
3105  real :: temp_v(G%domain%niglobal,G%domain%njglobal+1)
3106 
3107  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3108  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3109  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3110 
3111  ! For now, there are no radiation conditions applied to the thicknesses, since
3112  ! the thicknesses might not be physically motivated. Instead, sponges should be
3113  ! used to enforce the near-boundary layer structure.
3114 
3115  if (associated(tv%T)) then
3116 
3117  call pass_var(tv%T, g%Domain)
3118  call pass_var(tv%S, g%Domain)
3119 
3120  do n=1,obc%number_of_segments
3121  segment => obc%segment(n)
3122  if (.not. segment%on_pe) cycle
3123 
3124  if (segment%direction == obc_direction_e) then
3125  i=segment%HI%IsdB
3126  do k=1,g%ke ; do j=segment%HI%jsd,segment%HI%jed
3127  tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k)
3128  enddo ; enddo
3129  elseif (segment%direction == obc_direction_w) then
3130  i=segment%HI%IsdB
3131  do k=1,g%ke ; do j=segment%HI%jsd,segment%HI%jed
3132  tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k)
3133  enddo ; enddo
3134  elseif (segment%direction == obc_direction_n) then
3135  j=segment%HI%JsdB
3136  do k=1,g%ke ; do i=segment%HI%isd,segment%HI%ied
3137  tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k)
3138  enddo ; enddo
3139  elseif (segment%direction == obc_direction_s) then
3140  j=segment%HI%JsdB
3141  do k=1,g%ke ; do i=segment%HI%isd,segment%HI%ied
3142  tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k)
3143  enddo ; enddo
3144  endif
3145  enddo
3146  endif
3147 

References obc_direction_e, obc_direction_n, obc_direction_s, and obc_direction_w.

◆ setup_obc_tracer_reservoirs()

subroutine mom_open_boundary::setup_obc_tracer_reservoirs ( type(ocean_grid_type), intent(inout)  G,
type(ocean_obc_type), pointer  OBC 
)
private

Make sure the OBC tracer reservoirs are initialized.

Parameters
[in,out]gOcean grid structure
obcOpen boundary control structure

Definition at line 1749 of file MOM_open_boundary.F90.

1749  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1750  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
1751  ! Local variables
1752  type(OBC_segment_type), pointer :: segment => null()
1753  integer :: i, j, k, m, n
1754 
1755  do n=1,obc%number_of_segments
1756  segment=>obc%segment(n)
1757  if (associated(segment%tr_Reg)) then
1758  if (segment%is_E_or_W) then
1759  i = segment%HI%IsdB
1760  do m=1,obc%ntr
1761  if (associated(segment%tr_Reg%Tr(m)%tres)) then
1762  do k=1,g%ke
1763  do j=segment%HI%jsd,segment%HI%jed
1764  obc%tres_x(i,j,k,m) = segment%tr_Reg%Tr(m)%t(i,j,k)
1765  enddo
1766  enddo
1767  endif
1768  enddo
1769  else
1770  j = segment%HI%JsdB
1771  do m=1,obc%ntr
1772  if (associated(segment%tr_Reg%Tr(m)%tres)) then
1773  do k=1,g%ke
1774  do i=segment%HI%isd,segment%HI%ied
1775  obc%tres_y(i,j,k,m) = segment%tr_Reg%Tr(m)%t(i,j,k)
1776  enddo
1777  enddo
1778  endif
1779  enddo
1780  endif
1781  endif
1782  enddo
1783 

Referenced by fill_temp_salt_segments().

Here is the caller graph for this function:

◆ setup_segment_indices()

subroutine mom_open_boundary::setup_segment_indices ( type(dyn_horgrid_type), intent(in)  G,
type(obc_segment_type), intent(inout)  seg,
integer, intent(in)  Is_obc,
integer, intent(in)  Ie_obc,
integer, intent(in)  Js_obc,
integer, intent(in)  Je_obc 
)

Define indices for segment and store in hor_index_type using global segment bounds corresponding to q-points.

Parameters
[in]ggrid type
[in,out]segOpen boundary segment
[in]is_obcQ-point global i-index of start of segment
[in]ie_obcQ-point global i-index of end of segment
[in]js_obcQ-point global j-index of start of segment
[in]je_obcQ-point global j-index of end of segment

Definition at line 799 of file MOM_open_boundary.F90.

799  type(dyn_horgrid_type), intent(in) :: G !< grid type
800  type(OBC_segment_type), intent(inout) :: seg !< Open boundary segment
801  integer, intent(in) :: Is_obc !< Q-point global i-index of start of segment
802  integer, intent(in) :: Ie_obc !< Q-point global i-index of end of segment
803  integer, intent(in) :: Js_obc !< Q-point global j-index of start of segment
804  integer, intent(in) :: Je_obc !< Q-point global j-index of end of segment
805  ! Local variables
806  integer :: Isg,Ieg,Jsg,Jeg
807 
808  ! Isg, Ieg will be I*_obc in global space
809  if (ie_obc<is_obc) then
810  isg=ie_obc;ieg=is_obc
811  else
812  isg=is_obc;ieg=ie_obc
813  endif
814  if (je_obc<js_obc) then
815  jsg=je_obc;jeg=js_obc
816  else
817  jsg=js_obc;jeg=je_obc
818  endif
819 
820  ! Global space I*_obc but sorted
821  seg%HI%IsgB = isg ; seg%HI%IegB = ieg
822  seg%HI%isg = isg+1 ; seg%HI%ieg = ieg
823  seg%HI%JsgB = jsg ; seg%HI%JegB = jeg
824  seg%HI%jsg = jsg+1 ; seg%HI%Jeg = jeg
825 
826  ! Move into local index space
827  isg = isg - g%idg_offset
828  jsg = jsg - g%jdg_offset
829  ieg = ieg - g%idg_offset
830  jeg = jeg - g%jdg_offset
831 
832  ! This is the i-extent of the segment on this PE.
833  ! The values are nonsense if the segment is not on this PE.
834  seg%HI%IsdB = min( max(isg, g%HI%IsdB), g%HI%IedB)
835  seg%HI%IedB = min( max(ieg, g%HI%IsdB), g%HI%IedB)
836  seg%HI%isd = min( max(isg+1, g%HI%isd), g%HI%ied)
837  seg%HI%ied = min( max(ieg, g%HI%isd), g%HI%ied)
838  seg%HI%IscB = min( max(isg, g%HI%IscB), g%HI%IecB)
839  seg%HI%IecB = min( max(ieg, g%HI%IscB), g%HI%IecB)
840  seg%HI%isc = min( max(isg+1, g%HI%isc), g%HI%iec)
841  seg%HI%iec = min( max(ieg, g%HI%isc), g%HI%iec)
842 
843  ! This is the j-extent of the segment on this PE.
844  ! The values are nonsense if the segment is not on this PE.
845  seg%HI%JsdB = min( max(jsg, g%HI%JsdB), g%HI%JedB)
846  seg%HI%JedB = min( max(jeg, g%HI%JsdB), g%HI%JedB)
847  seg%HI%jsd = min( max(jsg+1, g%HI%jsd), g%HI%jed)
848  seg%HI%jed = min( max(jeg, g%HI%jsd), g%HI%jed)
849  seg%HI%JscB = min( max(jsg, g%HI%JscB), g%HI%JecB)
850  seg%HI%JecB = min( max(jeg, g%HI%JscB), g%HI%JecB)
851  seg%HI%jsc = min( max(jsg+1, g%HI%jsc), g%HI%jec)
852  seg%HI%jec = min( max(jeg, g%HI%jsc), g%HI%jec)
853 

Referenced by setup_u_point_obc(), and setup_v_point_obc().

Here is the caller graph for this function:

◆ setup_u_point_obc()

subroutine mom_open_boundary::setup_u_point_obc ( type(ocean_obc_type), pointer  OBC,
type(dyn_horgrid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
character(len=*), intent(in)  segment_str,
integer, intent(in)  l_seg,
type(param_file_type), intent(in)  PF,
logical, intent(in)  reentrant_y 
)
private

Parse an OBC_SEGMENT_%%% string starting with "I=" and configure placement and type of OBC accordingly.

Parameters
obcOpen boundary control structure
[in]gOcean grid structure
[in]usA dimensional unit scaling type
[in]segment_strA string in form of "I=%,J=%:%,string"
[in]l_segwhich segment is this?
[in]pfParameter file handle
[in]reentrant_yis the domain reentrant in y?

Definition at line 858 of file MOM_open_boundary.F90.

858  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
859  type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure
860  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
861  character(len=*), intent(in) :: segment_str !< A string in form of "I=%,J=%:%,string"
862  integer, intent(in) :: l_seg !< which segment is this?
863  type(param_file_type), intent(in) :: PF !< Parameter file handle
864  logical, intent(in) :: reentrant_y !< is the domain reentrant in y?
865  ! Local variables
866  integer :: I_obc, Js_obc, Je_obc ! Position of segment in global index space
867  integer :: j, a_loop
868  character(len=32) :: action_str(8)
869  character(len=128) :: segment_param_str
870  real, allocatable, dimension(:) :: tnudge
871  ! This returns the global indices for the segment
872  call parse_segment_str(g%ieg, g%jeg, segment_str, i_obc, js_obc, je_obc, action_str, reentrant_y)
873 
874  call setup_segment_indices(g, obc%segment(l_seg),i_obc,i_obc,js_obc,je_obc)
875 
876  i_obc = i_obc - g%idg_offset ! Convert to local tile indices on this tile
877  js_obc = js_obc - g%jdg_offset ! Convert to local tile indices on this tile
878  je_obc = je_obc - g%jdg_offset ! Convert to local tile indices on this tile
879 
880  if (je_obc>js_obc) then
881  obc%segment(l_seg)%direction = obc_direction_e
882  elseif (je_obc<js_obc) then
883  obc%segment(l_seg)%direction = obc_direction_w
884  j=js_obc;js_obc=je_obc;je_obc=j
885  endif
886 
887  obc%segment(l_seg)%on_pe = .false.
888 
889  do a_loop = 1,8 ! up to 8 options available
890  if (len_trim(action_str(a_loop)) == 0) then
891  cycle
892  elseif (trim(action_str(a_loop)) == 'FLATHER') then
893  obc%segment(l_seg)%Flather = .true.
894  obc%segment(l_seg)%open = .true.
895  obc%Flather_u_BCs_exist_globally = .true.
896  obc%open_u_BCs_exist_globally = .true.
897  obc%segment%z_values_needed = .true.
898  obc%segment%u_values_needed = .true.
899  elseif (trim(action_str(a_loop)) == 'ORLANSKI') then
900  obc%segment(l_seg)%radiation = .true.
901  obc%segment(l_seg)%open = .true.
902  obc%open_u_BCs_exist_globally = .true.
903  obc%radiation_BCs_exist_globally = .true.
904  elseif (trim(action_str(a_loop)) == 'ORLANSKI_TAN') then
905  obc%segment(l_seg)%radiation = .true.
906  obc%segment(l_seg)%radiation_tan = .true.
907  obc%radiation_BCs_exist_globally = .true.
908  elseif (trim(action_str(a_loop)) == 'ORLANSKI_GRAD') then
909  obc%segment(l_seg)%radiation = .true.
910  obc%segment(l_seg)%radiation_grad = .true.
911  elseif (trim(action_str(a_loop)) == 'OBLIQUE') then
912  obc%segment(l_seg)%oblique = .true.
913  obc%segment(l_seg)%open = .true.
914  obc%oblique_BCs_exist_globally = .true.
915  obc%open_u_BCs_exist_globally = .true.
916  elseif (trim(action_str(a_loop)) == 'OBLIQUE_TAN') then
917  obc%segment(l_seg)%oblique = .true.
918  obc%segment(l_seg)%oblique_tan = .true.
919  obc%oblique_BCs_exist_globally = .true.
920  elseif (trim(action_str(a_loop)) == 'OBLIQUE_GRAD') then
921  obc%segment(l_seg)%oblique = .true.
922  obc%segment(l_seg)%oblique_grad = .true.
923  elseif (trim(action_str(a_loop)) == 'NUDGED') then
924  obc%segment(l_seg)%nudged = .true.
925  obc%nudged_u_BCs_exist_globally = .true.
926  obc%segment%u_values_needed = .true.
927  elseif (trim(action_str(a_loop)) == 'NUDGED_TAN') then
928  obc%segment(l_seg)%nudged_tan = .true.
929  obc%nudged_u_BCs_exist_globally = .true.
930  obc%segment%v_values_needed = .true.
931  elseif (trim(action_str(a_loop)) == 'NUDGED_GRAD') then
932  obc%segment(l_seg)%nudged_grad = .true.
933  obc%segment%g_values_needed = .true.
934  elseif (trim(action_str(a_loop)) == 'GRADIENT') then
935  obc%segment(l_seg)%gradient = .true.
936  obc%segment(l_seg)%open = .true.
937  obc%open_u_BCs_exist_globally = .true.
938  elseif (trim(action_str(a_loop)) == 'SIMPLE') then
939  obc%segment(l_seg)%specified = .true.
940  obc%specified_u_BCs_exist_globally = .true. ! This avoids deallocation
941  obc%segment%u_values_needed = .true.
942  elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then
943  obc%segment(l_seg)%specified_tan = .true.
944  else
945  call mom_error(fatal, "MOM_open_boundary.F90, setup_u_point_obc: "//&
946  "String '"//trim(action_str(a_loop))//"' not understood.")
947  endif
948  if (obc%segment(l_seg)%nudged .or. obc%segment(l_seg)%nudged_tan) then
949  write(segment_param_str(1:43),"('OBC_SEGMENT_',i3.3,'_VELOCITY_NUDGING_TIMESCALES')") l_seg
950  allocate(tnudge(2))
951  call get_param(pf, mdl, segment_param_str(1:43), tnudge, &
952  "Timescales in days for nudging along a segment, "//&
953  "for inflow, then outflow. Setting both to zero should "//&
954  "behave like SIMPLE obcs for the baroclinic velocities.", &
955  fail_if_missing=.true., default=0., units="days", scale=86400.0*us%s_to_T)
956  obc%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1)
957  obc%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2)
958  deallocate(tnudge)
959  endif
960 
961  enddo ! a_loop
962 
963  obc%segment(l_seg)%is_E_or_W_2 = .true.
964 
965  if (i_obc<=g%HI%IsdB+1 .or. i_obc>=g%HI%IedB-1) return ! Boundary is not on tile
966  if (je_obc<=g%HI%JsdB .or. js_obc>=g%HI%JedB) return ! Segment is not on tile
967 
968  obc%segment(l_seg)%on_pe = .true.
969  obc%segment(l_seg)%is_E_or_W = .true.
970 
971  do j=g%HI%jsd, g%HI%jed
972  if (j>js_obc .and. j<=je_obc) then
973  obc%segnum_u(i_obc,j) = l_seg
974  endif
975  enddo
976  obc%segment(l_seg)%Is_obc = i_obc
977  obc%segment(l_seg)%Ie_obc = i_obc
978  obc%segment(l_seg)%Js_obc = js_obc
979  obc%segment(l_seg)%Je_obc = je_obc
980  call allocate_obc_segment_data(obc, obc%segment(l_seg))
981 
982  if (obc%segment(l_seg)%oblique .and. obc%segment(l_seg)%radiation) &
983  call mom_error(fatal, "MOM_open_boundary.F90, setup_u_point_obc: \n"//&
984  "Orlanski and Oblique OBC options cannot be used together on one segment.")
985 
986  if (obc%segment(l_seg)%u_values_needed .or. obc%segment(l_seg)%v_values_needed .or. &
987  obc%segment(l_seg)%t_values_needed .or. obc%segment(l_seg)%s_values_needed .or. &
988  obc%segment(l_seg)%z_values_needed .or. obc%segment(l_seg)%g_values_needed) &
989  obc%segment(l_seg)%values_needed = .true.

References allocate_obc_segment_data(), mdl, mom_error_handler::mom_error(), obc_direction_e, obc_direction_w, parse_segment_str(), and setup_segment_indices().

Referenced by open_boundary_config().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ setup_v_point_obc()

subroutine mom_open_boundary::setup_v_point_obc ( type(ocean_obc_type), pointer  OBC,
type(dyn_horgrid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
character(len=*), intent(in)  segment_str,
integer, intent(in)  l_seg,
type(param_file_type), intent(in)  PF,
logical, intent(in)  reentrant_x 
)
private

Parse an OBC_SEGMENT_%%% string starting with "J=" and configure placement and type of OBC accordingly.

Parameters
obcOpen boundary control structure
[in]gOcean grid structure
[in]usA dimensional unit scaling type
[in]segment_strA string in form of "J=%,I=%:%,string"
[in]l_segwhich segment is this?
[in]pfParameter file handle
[in]reentrant_xis the domain reentrant in x?

Definition at line 994 of file MOM_open_boundary.F90.

994  type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
995  type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure
996  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
997  character(len=*), intent(in) :: segment_str !< A string in form of "J=%,I=%:%,string"
998  integer, intent(in) :: l_seg !< which segment is this?
999  type(param_file_type), intent(in) :: PF !< Parameter file handle
1000  logical, intent(in) :: reentrant_x !< is the domain reentrant in x?
1001  ! Local variables
1002  integer :: J_obc, Is_obc, Ie_obc ! Position of segment in global index space
1003  integer :: i, a_loop
1004  character(len=32) :: action_str(8)
1005  character(len=128) :: segment_param_str
1006  real, allocatable, dimension(:) :: tnudge
1007 
1008  ! This returns the global indices for the segment
1009  call parse_segment_str(g%ieg, g%jeg, segment_str, j_obc, is_obc, ie_obc, action_str, reentrant_x)
1010 
1011  call setup_segment_indices(g, obc%segment(l_seg),is_obc,ie_obc,j_obc,j_obc)
1012 
1013  j_obc = j_obc - g%jdg_offset ! Convert to local tile indices on this tile
1014  is_obc = is_obc - g%idg_offset ! Convert to local tile indices on this tile
1015  ie_obc = ie_obc - g%idg_offset ! Convert to local tile indices on this tile
1016 
1017  if (ie_obc>is_obc) then
1018  obc%segment(l_seg)%direction = obc_direction_s
1019  elseif (ie_obc<is_obc) then
1020  obc%segment(l_seg)%direction = obc_direction_n
1021  i=is_obc;is_obc=ie_obc;ie_obc=i
1022  endif
1023 
1024  obc%segment(l_seg)%on_pe = .false.
1025 
1026  do a_loop = 1,8
1027  if (len_trim(action_str(a_loop)) == 0) then
1028  cycle
1029  elseif (trim(action_str(a_loop)) == 'FLATHER') then
1030  obc%segment(l_seg)%Flather = .true.
1031  obc%segment(l_seg)%open = .true.
1032  obc%Flather_v_BCs_exist_globally = .true.
1033  obc%open_v_BCs_exist_globally = .true.
1034  obc%segment%z_values_needed = .true.
1035  obc%segment%v_values_needed = .true.
1036  elseif (trim(action_str(a_loop)) == 'ORLANSKI') then
1037  obc%segment(l_seg)%radiation = .true.
1038  obc%segment(l_seg)%open = .true.
1039  obc%open_v_BCs_exist_globally = .true.
1040  obc%radiation_BCs_exist_globally = .true.
1041  elseif (trim(action_str(a_loop)) == 'ORLANSKI_TAN') then
1042  obc%segment(l_seg)%radiation = .true.
1043  obc%segment(l_seg)%radiation_tan = .true.
1044  obc%radiation_BCs_exist_globally = .true.
1045  elseif (trim(action_str(a_loop)) == 'ORLANSKI_GRAD') then
1046  obc%segment(l_seg)%radiation = .true.
1047  obc%segment(l_seg)%radiation_grad = .true.
1048  elseif (trim(action_str(a_loop)) == 'OBLIQUE') then
1049  obc%segment(l_seg)%oblique = .true.
1050  obc%segment(l_seg)%open = .true.
1051  obc%oblique_BCs_exist_globally = .true.
1052  obc%open_v_BCs_exist_globally = .true.
1053  elseif (trim(action_str(a_loop)) == 'OBLIQUE_TAN') then
1054  obc%segment(l_seg)%oblique = .true.
1055  obc%segment(l_seg)%oblique_tan = .true.
1056  obc%oblique_BCs_exist_globally = .true.
1057  elseif (trim(action_str(a_loop)) == 'OBLIQUE_GRAD') then
1058  obc%segment(l_seg)%oblique = .true.
1059  obc%segment(l_seg)%oblique_grad = .true.
1060  elseif (trim(action_str(a_loop)) == 'NUDGED') then
1061  obc%segment(l_seg)%nudged = .true.
1062  obc%nudged_v_BCs_exist_globally = .true.
1063  obc%segment%v_values_needed = .true.
1064  elseif (trim(action_str(a_loop)) == 'NUDGED_TAN') then
1065  obc%segment(l_seg)%nudged_tan = .true.
1066  obc%nudged_v_BCs_exist_globally = .true.
1067  obc%segment%u_values_needed = .true.
1068  elseif (trim(action_str(a_loop)) == 'NUDGED_GRAD') then
1069  obc%segment(l_seg)%nudged_grad = .true.
1070  obc%segment%g_values_needed = .true.
1071  elseif (trim(action_str(a_loop)) == 'GRADIENT') then
1072  obc%segment(l_seg)%gradient = .true.
1073  obc%segment(l_seg)%open = .true.
1074  obc%open_v_BCs_exist_globally = .true.
1075  elseif (trim(action_str(a_loop)) == 'SIMPLE') then
1076  obc%segment(l_seg)%specified = .true.
1077  obc%specified_v_BCs_exist_globally = .true. ! This avoids deallocation
1078  obc%segment%v_values_needed = .true.
1079  elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then
1080  obc%segment(l_seg)%specified_tan = .true.
1081  else
1082  call mom_error(fatal, "MOM_open_boundary.F90, setup_v_point_obc: "//&
1083  "String '"//trim(action_str(a_loop))//"' not understood.")
1084  endif
1085  if (obc%segment(l_seg)%nudged .or. obc%segment(l_seg)%nudged_tan) then
1086  write(segment_param_str(1:43),"('OBC_SEGMENT_',i3.3,'_VELOCITY_NUDGING_TIMESCALES')") l_seg
1087  allocate(tnudge(2))
1088  call get_param(pf, mdl, segment_param_str(1:43), tnudge, &
1089  "Timescales in days for nudging along a segment, "//&
1090  "for inflow, then outflow. Setting both to zero should "//&
1091  "behave like SIMPLE obcs for the baroclinic velocities.", &
1092  fail_if_missing=.true., default=0., units="days", scale=86400.0*us%s_to_T)
1093  obc%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1)
1094  obc%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2)
1095  deallocate(tnudge)
1096  endif
1097 
1098  enddo ! a_loop
1099 
1100  if (j_obc<=g%HI%JsdB+1 .or. j_obc>=g%HI%JedB-1) return ! Boundary is not on tile
1101  if (ie_obc<=g%HI%IsdB .or. is_obc>=g%HI%IedB) return ! Segment is not on tile
1102 
1103  obc%segment(l_seg)%on_pe = .true.
1104  obc%segment(l_seg)%is_N_or_S = .true.
1105 
1106  do i=g%HI%isd, g%HI%ied
1107  if (i>is_obc .and. i<=ie_obc) then
1108  obc%segnum_v(i,j_obc) = l_seg
1109  endif
1110  enddo
1111  obc%segment(l_seg)%Is_obc = is_obc
1112  obc%segment(l_seg)%Ie_obc = ie_obc
1113  obc%segment(l_seg)%Js_obc = j_obc
1114  obc%segment(l_seg)%Je_obc = j_obc
1115  call allocate_obc_segment_data(obc, obc%segment(l_seg))
1116 
1117  if (obc%segment(l_seg)%oblique .and. obc%segment(l_seg)%radiation) &
1118  call mom_error(fatal, "MOM_open_boundary.F90, setup_v_point_obc: \n"//&
1119  "Orlanski and Oblique OBC options cannot be used together on one segment.")
1120 
1121  if (obc%segment(l_seg)%u_values_needed .or. obc%segment(l_seg)%v_values_needed .or. &
1122  obc%segment(l_seg)%t_values_needed .or. obc%segment(l_seg)%s_values_needed .or. &
1123  obc%segment(l_seg)%z_values_needed .or. obc%segment(l_seg)%g_values_needed) &
1124  obc%segment(l_seg)%values_needed = .true.

References allocate_obc_segment_data(), mdl, mom_error_handler::mom_error(), obc_direction_n, obc_direction_s, parse_segment_str(), and setup_segment_indices().

Referenced by open_boundary_config().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ update_obc_segment_data()

subroutine, public mom_open_boundary::update_obc_segment_data ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(ocean_obc_type), pointer  OBC,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h,
type(time_type), intent(in)  Time 
)

Update the OBC values on the segments.

Parameters
[in]gOcean grid structure
[in]gvOcean vertical grid structure
[in]usA dimensional unit scaling type
obcOpen boundary structure
[in]tvThermodynamics structure
[in,out]hThickness [m]
[in]timeModel time

Definition at line 3401 of file MOM_open_boundary.F90.

3401  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
3402  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
3403  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3404  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
3405  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
3406  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h !< Thickness [m]
3407  type(time_type), intent(in) :: Time !< Model time
3408  ! Local variables
3409  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed
3410  integer :: IsdB, IedB, JsdB, JedB, n, m, nz
3411  character(len=40) :: mdl = "set_OBC_segment_data" ! This subroutine's name.
3412  character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path
3413  type(OBC_segment_type), pointer :: segment => null()
3414  integer, dimension(4) :: siz,siz2
3415  real :: sumh ! column sum of thicknesses [m]
3416  integer :: ni_seg, nj_seg ! number of src gridpoints along the segments
3417  integer :: i2, j2 ! indices for referencing local domain array
3418  integer :: is_obc, ie_obc, js_obc, je_obc ! segment indices within local domain
3419  integer :: ishift, jshift ! offsets for staggered locations
3420  real, dimension(:,:), pointer :: seg_vel => null() ! pointer to segment velocity array
3421  real, dimension(:,:), pointer :: seg_trans => null() ! pointer to segment transport array
3422  real, dimension(:,:,:), allocatable :: tmp_buffer
3423  real, dimension(:), allocatable :: h_stack
3424  integer :: is_obc2, js_obc2
3425  real :: net_H_src, net_H_int, scl_fac
3426  real, pointer, dimension(:,:) :: normal_trans_bt=>null() ! barotropic transport
3427 
3428  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3429  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3430  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3431  nz=g%ke
3432 
3433  if (.not. associated(obc)) return
3434 
3435  do n = 1, obc%number_of_segments
3436  segment => obc%segment(n)
3437 
3438  if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain
3439 
3440  ni_seg = segment%ie_obc-segment%is_obc+1
3441  nj_seg = segment%je_obc-segment%js_obc+1
3442  is_obc = max(segment%is_obc,isd-1)
3443  ie_obc = min(segment%ie_obc,ied)
3444  js_obc = max(segment%js_obc,jsd-1)
3445  je_obc = min(segment%je_obc,jed)
3446 
3447 ! Calculate auxiliary fields at staggered locations.
3448 ! Segment indices are on q points:
3449 !
3450 ! |-----------|------------|-----------|-----------| J_obc
3451 ! Is_obc Ie_obc
3452 !
3453 ! i2 has to start at Is_obc+1 and end at Ie_obc.
3454 ! j2 is J_obc and jshift has to be +1 at both the north and south.
3455 
3456  ! calculate auxiliary fields at staggered locations
3457  ishift=0;jshift=0
3458  if (segment%is_E_or_W) then
3459  allocate(normal_trans_bt(segment%HI%IsdB:segment%HI%IedB,segment%HI%jsd:segment%HI%jed))
3460  normal_trans_bt(:,:)=0.0
3461  if (segment%direction == obc_direction_w) ishift=1
3462  i=segment%HI%IsdB
3463  do j=segment%HI%jsd,segment%HI%jed
3464  segment%Cg(i,j) = sqrt(gv%g_prime(1)*g%bathyT(i+ishift,j))
3465  segment%Htot(i,j)=0.0
3466  do k=1,g%ke
3467  segment%h(i,j,k) = h(i+ishift,j,k)
3468  segment%Htot(i,j)=segment%Htot(i,j)+segment%h(i,j,k)
3469  enddo
3470  enddo
3471  else! (segment%direction == OBC_DIRECTION_N .or. segment%direction == OBC_DIRECTION_S)
3472  allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB))
3473  normal_trans_bt(:,:)=0.0
3474  if (segment%direction == obc_direction_s) jshift=1
3475  j=segment%HI%JsdB
3476  do i=segment%HI%isd,segment%HI%ied
3477  segment%Cg(i,j) = sqrt(gv%g_prime(1)*g%bathyT(i,j+jshift))
3478  segment%Htot(i,j)=0.0
3479  do k=1,g%ke
3480  segment%h(i,j,k) = h(i,j+jshift,k)
3481  segment%Htot(i,j)=segment%Htot(i,j)+segment%h(i,j,k)
3482  enddo
3483  enddo
3484  endif
3485 
3486  allocate(h_stack(g%ke))
3487  h_stack(:) = 0.0
3488  do m = 1,segment%num_fields
3489  if (segment%field(m)%fid > 0) then
3490  siz(1)=size(segment%field(m)%buffer_src,1)
3491  siz(2)=size(segment%field(m)%buffer_src,2)
3492  siz(3)=size(segment%field(m)%buffer_src,3)
3493  if (.not.associated(segment%field(m)%buffer_dst)) then
3494  if (siz(3) /= segment%field(m)%nk_src) call mom_error(fatal,'nk_src inconsistency')
3495  if (segment%field(m)%nk_src > 1) then
3496  if (segment%is_E_or_W) then
3497  if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
3498  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3499  else
3500  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
3501  endif
3502  if (segment%field(m)%name == 'U') then
3503  allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
3504  segment%field(m)%bt_vel(:,:)=0.0
3505  endif
3506  else
3507  if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
3508  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3509  else
3510  allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
3511  endif
3512  if (segment%field(m)%name == 'V') then
3513  allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
3514  segment%field(m)%bt_vel(:,:)=0.0
3515  endif
3516  endif
3517  else
3518  if (segment%is_E_or_W) then
3519  if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
3520  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
3521  else
3522  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,1))
3523  endif
3524  if (segment%field(m)%name == 'U') then
3525  allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
3526  segment%field(m)%bt_vel(:,:)=0.0
3527  endif
3528  else
3529  if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
3530  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
3531  else
3532  allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,1))
3533  endif
3534  if (segment%field(m)%name == 'V') then
3535  allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
3536  segment%field(m)%bt_vel(:,:)=0.0
3537  endif
3538  endif
3539  endif
3540  segment%field(m)%buffer_dst(:,:,:)=0.0
3541  endif
3542  ! read source data interpolated to the current model time
3543  if (siz(1)==1) then
3544  if (obc%brushcutter_mode) then
3545  allocate(tmp_buffer(1,nj_seg*2-1,segment%field(m)%nk_src)) ! segment data is currrently on supergrid
3546  else
3547  allocate(tmp_buffer(1,nj_seg,segment%field(m)%nk_src)) ! segment data is currrently on supergrid
3548  endif
3549  else
3550  if (obc%brushcutter_mode) then
3551  allocate(tmp_buffer(ni_seg*2-1,1,segment%field(m)%nk_src)) ! segment data is currrently on supergrid
3552  else
3553  allocate(tmp_buffer(ni_seg,1,segment%field(m)%nk_src)) ! segment data is currrently on supergrid
3554  endif
3555  endif
3556 
3557  call time_interp_external(segment%field(m)%fid,time, tmp_buffer)
3558  if (obc%brushcutter_mode) then
3559  if (segment%is_E_or_W) then
3560  if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
3561  segment%field(m)%buffer_src(is_obc,:,:) = &
3562  tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset)+1:2,:)
3563  else
3564  segment%field(m)%buffer_src(is_obc,:,:) = &
3565  tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset):2,:)
3566  endif
3567  else
3568  if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
3569  segment%field(m)%buffer_src(:,js_obc,:) = &
3570  tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset)+1:2,1,:)
3571  else
3572  segment%field(m)%buffer_src(:,js_obc,:) = &
3573  tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset):2,1,:)
3574  endif
3575  endif
3576  else
3577  if (segment%is_E_or_W) then
3578  if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
3579  segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset+1,:)
3580  else
3581  segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset,:)
3582  endif
3583  else
3584  if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
3585  segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset+1,1,:)
3586  else
3587  segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset,1,:)
3588  endif
3589  endif
3590  endif
3591  if (segment%field(m)%nk_src > 1) then
3592  call time_interp_external(segment%field(m)%fid_dz,time, tmp_buffer)
3593  if (obc%brushcutter_mode) then
3594  if (segment%is_E_or_W) then
3595  if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
3596  segment%field(m)%dz_src(is_obc,:,:) = &
3597  tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset)+1:2,:)
3598  else
3599  segment%field(m)%dz_src(is_obc,:,:) = &
3600  tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset):2,:)
3601  endif
3602  else
3603  if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
3604  segment%field(m)%dz_src(:,js_obc,:) = &
3605  tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset)+1:2,1,:)
3606  else
3607  segment%field(m)%dz_src(:,js_obc,:) = &
3608  tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset):2,1,:)
3609  endif
3610  endif
3611  else
3612  if (segment%is_E_or_W) then
3613  if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
3614  segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset+1,:)
3615  else
3616  segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset,:)
3617  endif
3618  else
3619  if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
3620  segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset+1,1,:)
3621  else
3622  segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset,1,:)
3623  endif
3624  endif
3625  endif
3626 
3627  call adjustsegmentetatofitbathymetry(g,gv,us,segment,m)
3628 
3629  if (segment%is_E_or_W) then
3630  ishift=1
3631  if (segment%direction == obc_direction_e) ishift=0
3632  i=is_obc
3633  if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then
3634  ! Do q points for the whole segment
3635  do j=max(js_obc,jsd),min(je_obc,jed-1)
3636  ! Using the h remapping approach
3637  ! Pretty sure we need to check for source/target grid consistency here
3638  segment%field(m)%buffer_dst(i,j,:)=0.0 ! initialize remap destination buffer
3639  if (g%mask2dCu(i,j)>0. .and. g%mask2dCu(i,j+1)>0.) then
3640  h_stack(:) = 0.5*(h(i+ishift,j,:) + h(i+ishift,j+1,:))
3641  call remapping_core_h(obc%remap_CS, &
3642  segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3643  segment%field(m)%buffer_src(i,j,:), &
3644  g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3645  elseif (g%mask2dCu(i,j)>0.) then
3646  h_stack(:) = h(i+ishift,j,:)
3647  call remapping_core_h(obc%remap_CS, &
3648  segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3649  segment%field(m)%buffer_src(i,j,:), &
3650  g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3651  elseif (g%mask2dCu(i,j+1)>0.) then
3652  h_stack(:) = h(i+ishift,j+1,:)
3653  call remapping_core_h(obc%remap_CS, &
3654  segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3655  segment%field(m)%buffer_src(i,j,:), &
3656  g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3657  endif
3658  enddo
3659  else
3660  do j=js_obc+1,je_obc
3661  ! Using the h remapping approach
3662  ! Pretty sure we need to check for source/target grid consistency here
3663  segment%field(m)%buffer_dst(i,j,:)=0.0 ! initialize remap destination buffer
3664  if (g%mask2dCu(i,j)>0.) then
3665  net_h_src = sum( segment%field(m)%dz_src(i,j,:) )
3666  net_h_int = sum( h(i+ishift,j,:) )
3667  scl_fac = net_h_int / net_h_src
3668  call remapping_core_h(obc%remap_CS, &
3669  segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,j,:), &
3670  segment%field(m)%buffer_src(i,j,:), &
3671  g%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(i,j,:))
3672  endif
3673  enddo
3674  endif
3675  else
3676  jshift=1
3677  if (segment%direction == obc_direction_n) jshift=0
3678  j=js_obc
3679  if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then
3680  ! Do q points for the whole segment
3681  do i=max(is_obc,isd),min(ie_obc,ied-1)
3682  segment%field(m)%buffer_dst(i,j,:)=0.0 ! initialize remap destination buffer
3683  if (g%mask2dCv(i,j)>0. .and. g%mask2dCv(i+1,j)>0.) then
3684  ! Using the h remapping approach
3685  ! Pretty sure we need to check for source/target grid consistency here
3686  h_stack(:) = 0.5*(h(i,j+jshift,:) + h(i+1,j+jshift,:))
3687  call remapping_core_h(obc%remap_CS, &
3688  segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3689  segment%field(m)%buffer_src(i,j,:), &
3690  g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3691  elseif (g%mask2dCv(i,j)>0.) then
3692  h_stack(:) = h(i,j+jshift,:)
3693  call remapping_core_h(obc%remap_CS, &
3694  segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3695  segment%field(m)%buffer_src(i,j,:), &
3696  g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3697  elseif (g%mask2dCv(i+1,j)>0.) then
3698  h_stack(:) = h(i+1,j+jshift,:)
3699  call remapping_core_h(obc%remap_CS, &
3700  segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3701  segment%field(m)%buffer_src(i,j,:), &
3702  g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3703  endif
3704  enddo
3705  else
3706  do i=is_obc+1,ie_obc
3707  ! Using the h remapping approach
3708  ! Pretty sure we need to check for source/target grid consistency here
3709  segment%field(m)%buffer_dst(i,j,:)=0.0 ! initialize remap destination buffer
3710  if (g%mask2dCv(i,j)>0.) then
3711  net_h_src = sum( segment%field(m)%dz_src(i,j,:) )
3712  net_h_int = sum( h(i,j+jshift,:) )
3713  scl_fac = net_h_int / net_h_src
3714  call remapping_core_h(obc%remap_CS, &
3715  segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3716  segment%field(m)%buffer_src(i,j,:), &
3717  g%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,j,:))
3718  endif
3719  enddo
3720  endif
3721  endif
3722  else ! 2d data
3723  segment%field(m)%buffer_dst(:,:,1) = segment%field(m)%buffer_src(:,:,1) ! initialize remap destination buffer
3724  endif
3725  deallocate(tmp_buffer)
3726  else ! fid <= 0 (Uniform value)
3727  if (.not. associated(segment%field(m)%buffer_dst)) then
3728  if (segment%is_E_or_W) then
3729  if (segment%field(m)%name == 'V') then
3730  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3731  allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc))
3732  elseif (segment%field(m)%name == 'U') then
3733  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
3734  allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
3735  elseif (segment%field(m)%name == 'DVDX') then
3736  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3737  elseif (segment%field(m)%name == 'SSH') then
3738  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
3739  else
3740  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
3741  endif
3742  else
3743  if (segment%field(m)%name == 'U') then
3744  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3745  allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc))
3746  elseif (segment%field(m)%name == 'V') then
3747  allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
3748  allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
3749  elseif (segment%field(m)%name == 'DUDY') then
3750  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3751  elseif (segment%field(m)%name == 'SSH') then
3752  allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
3753  else
3754  allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
3755  endif
3756  endif
3757  segment%field(m)%buffer_dst(:,:,:) = segment%field(m)%value
3758  if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then
3759  segment%field(m)%bt_vel(:,:) = segment%field(m)%value
3760  endif
3761  endif
3762  endif
3763 
3764  if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then
3765  if (segment%field(m)%fid>0) then ! calculate external BT velocity and transport if needed
3766  if (trim(segment%field(m)%name) == 'U' .and. segment%is_E_or_W) then
3767  i=is_obc
3768  do j=js_obc+1,je_obc
3769  normal_trans_bt(i,j) = 0.0
3770  do k=1,g%ke
3771  segment%normal_vel(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)
3772  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) * &
3773  g%dyCu(i,j)
3774  normal_trans_bt(i,j) = normal_trans_bt(i,j) + segment%normal_trans(i,j,k)
3775  enddo
3776  segment%normal_vel_bt(i,j) = normal_trans_bt(i,j) / (max(segment%Htot(i,j),1.e-12) * g%dyCu(i,j))
3777  if (associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,j,:) = segment%normal_vel(i,j,:)
3778  enddo
3779  elseif (trim(segment%field(m)%name) == 'V' .and. segment%is_N_or_S) then
3780  j=js_obc
3781  do i=is_obc+1,ie_obc
3782  normal_trans_bt(i,j) = 0.0
3783  do k=1,g%ke
3784  segment%normal_vel(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)
3785  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) * &
3786  g%dxCv(i,j)
3787  normal_trans_bt(i,j) = normal_trans_bt(i,j) + segment%normal_trans(i,j,k)
3788  enddo
3789  segment%normal_vel_bt(i,j) = normal_trans_bt(i,j) / (max(segment%Htot(i,j),1.e-12) * g%dxCv(i,j))
3790  if (associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,j,:) = segment%normal_vel(i,j,:)
3791  enddo
3792  elseif (trim(segment%field(m)%name) == 'V' .and. segment%is_E_or_W .and. &
3793  associated(segment%tangential_vel)) then
3794  i=is_obc
3795  do j=js_obc,je_obc
3796  do k=1,g%ke
3797  segment%tangential_vel(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)
3798  enddo
3799  if (associated(segment%nudged_tangential_vel)) &
3800  segment%nudged_tangential_vel(i,j,:) = segment%tangential_vel(i,j,:)
3801  enddo
3802  elseif (trim(segment%field(m)%name) == 'U' .and. segment%is_N_or_S .and. &
3803  associated(segment%tangential_vel)) then
3804  j=js_obc
3805  do i=is_obc,ie_obc
3806  do k=1,g%ke
3807  segment%tangential_vel(i,j,k) = us%m_s_to_L_T*segment%field(m)%buffer_dst(i,j,k)
3808  enddo
3809  if (associated(segment%nudged_tangential_vel)) &
3810  segment%nudged_tangential_vel(i,j,:) = segment%tangential_vel(i,j,:)
3811  enddo
3812  elseif (trim(segment%field(m)%name) == 'DVDX' .and. segment%is_E_or_W .and. &
3813  associated(segment%tangential_grad)) then
3814  i=is_obc
3815  do j=js_obc,je_obc
3816  do k=1,g%ke
3817  segment%tangential_grad(i,j,k) = us%T_to_s*segment%field(m)%buffer_dst(i,j,k)
3818  enddo
3819  enddo
3820  elseif (trim(segment%field(m)%name) == 'DUDY' .and. segment%is_N_or_S .and. &
3821  associated(segment%tangential_grad)) then
3822  j=js_obc
3823  do i=is_obc,ie_obc
3824  do k=1,g%ke
3825  segment%tangential_grad(i,j,k) = us%T_to_s*segment%field(m)%buffer_dst(i,j,k)
3826  enddo
3827  enddo
3828  endif
3829  endif
3830  endif
3831 
3832  ! from this point on, data are entirely on segments - will
3833  ! write all segment loops as 2d loops.
3834  if (segment%is_E_or_W) then
3835  js_obc2 = js_obc+1
3836  is_obc2 = is_obc
3837  else
3838  js_obc2 = js_obc
3839  is_obc2 = is_obc+1
3840  endif
3841  if (segment%is_N_or_S) then
3842  is_obc2 = is_obc+1
3843  js_obc2 = js_obc
3844  else
3845  is_obc2 = is_obc
3846  js_obc2 = js_obc+1
3847  endif
3848 
3849  if (trim(segment%field(m)%name) == 'SSH') then
3850  do j=js_obc2,je_obc
3851  do i=is_obc2,ie_obc
3852  segment%eta(i,j) = segment%field(m)%buffer_dst(i,j,1)
3853  enddo
3854  enddo
3855  endif
3856 
3857  if (trim(segment%field(m)%name) == 'TEMP') then
3858  if (associated(segment%field(m)%buffer_dst)) then
3859  do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc
3860  segment%tr_Reg%Tr(1)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3861  enddo ; enddo ; enddo
3862  if (.not. segment%tr_Reg%Tr(1)%is_initialized) then
3863  ! if the tracer reservoir has not yet been initialized, then set to external value.
3864  do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc
3865  segment%tr_Reg%Tr(1)%tres(i,j,k) = segment%tr_Reg%Tr(1)%t(i,j,k)
3866  enddo ; enddo ; enddo
3867  segment%tr_Reg%Tr(1)%is_initialized=.true.
3868  endif
3869  else
3870  segment%tr_Reg%Tr(1)%OBC_inflow_conc = segment%field(m)%value
3871  endif
3872  elseif (trim(segment%field(m)%name) == 'SALT') then
3873  if (associated(segment%field(m)%buffer_dst)) then
3874  do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc
3875  segment%tr_Reg%Tr(2)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3876  enddo ; enddo ; enddo
3877  if (.not. segment%tr_Reg%Tr(2)%is_initialized) then
3878  !if the tracer reservoir has not yet been initialized, then set to external value.
3879  do k=1,nz; do j=js_obc2, je_obc; do i=is_obc2,ie_obc
3880  segment%tr_Reg%Tr(2)%tres(i,j,k) = segment%tr_Reg%Tr(2)%t(i,j,k)
3881  enddo ; enddo ; enddo
3882  segment%tr_Reg%Tr(2)%is_initialized=.true.
3883  endif
3884  else
3885  segment%tr_Reg%Tr(2)%OBC_inflow_conc = segment%field(m)%value
3886  endif
3887  endif
3888 
3889  enddo ! end field loop
3890  deallocate(h_stack)
3891  deallocate(normal_trans_bt)
3892 
3893  enddo ! end segment loop
3894 

References adjustsegmentetatofitbathymetry(), obc_direction_e, obc_direction_n, obc_direction_s, and obc_direction_w.

Referenced by mom_state_initialization::mom_initialize_state().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ update_segment_tracer_reservoirs()

subroutine, public mom_open_boundary::update_segment_tracer_reservoirs ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  uhr,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  vhr,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(ocean_obc_type), pointer  OBC,
real, intent(in)  dt,
type(tracer_registry_type), pointer  Reg 
)

Update the OBC tracer reservoirs after the tracers have been updated.

Parameters
[in]gThe ocean's grid structure
[in]gvOcean vertical grid structure
[in]uhraccumulated volume/mass flux through the zonal face [H L2 ~> m3 or kg]
[in]vhraccumulated volume/mass flux through the meridional face [H L2 ~> m3 or kg]
[in]hlayer thickness after advection [H ~> m or kg m-2]
obcOpen boundary structure
[in]dttime increment [T ~> s]
regpointer to tracer registry

Definition at line 4513 of file MOM_open_boundary.F90.

4513  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
4514  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
4515  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uhr !< accumulated volume/mass flux through
4516  !! the zonal face [H L2 ~> m3 or kg]
4517  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vhr !< accumulated volume/mass flux through
4518  !! the meridional face [H L2 ~> m3 or kg]
4519  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness after advection
4520  !! [H ~> m or kg m-2]
4521  type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
4522  real, intent(in) :: dt !< time increment [T ~> s]
4523  type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry
4524  ! Local variables
4525  type(OBC_segment_type), pointer :: segment=>null()
4526  real :: u_L_in, u_L_out ! The zonal distance moved in or out of a cell [L ~> m]
4527  real :: v_L_in, v_L_out ! The meridional distance moved in or out of a cell [L ~> m]
4528  real :: fac1 ! The denominator of the expression for tracer updates [nondim]
4529  integer :: i, j, k, m, n, ntr, nz
4530  integer :: ishift, idir, jshift, jdir
4531 
4532  nz = gv%ke
4533  ntr = reg%ntr
4534  if (associated(obc)) then ; if (obc%OBC_pe) then ; do n=1,obc%number_of_segments
4535  segment=>obc%segment(n)
4536  if (.not. associated(segment%tr_Reg)) cycle
4537  if (segment%is_E_or_W) then
4538  do j=segment%HI%jsd,segment%HI%jed
4539  i = segment%HI%IsdB
4540  ! ishift+I corresponds to the nearest interior tracer cell index
4541  ! idir switches the sign of the flow so that positive is into the reservoir
4542  if (segment%direction == obc_direction_w) then
4543  ishift = 1 ; idir = -1
4544  else
4545  ishift = 0 ; idir = 1
4546  endif
4547  ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
4548  do m=1,ntr ; if (associated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
4549  u_l_out = max(0.0, (idir*uhr(i,j,k))*segment%Tr_InvLscale_out / (h(i+ishift,j,k)*g%dyCu(i,j)))
4550  u_l_in = min(0.0, (idir*uhr(i,j,k))*segment%Tr_InvLscale_in / (h(i+ishift,j,k)*g%dyCu(i,j)))
4551  fac1 = 1.0 + (u_l_out-u_l_in)
4552  segment%tr_Reg%Tr(m)%tres(i,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,j,k) + &
4553  (u_l_out*reg%Tr(m)%t(i+ishift,j,k) - &
4554  u_l_in*segment%tr_Reg%Tr(m)%t(i,j,k)))
4555  if (associated(obc%tres_x)) obc%tres_x(i,j,k,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
4556  enddo ; endif ; enddo
4557  enddo
4558  else
4559  do i=segment%HI%isd,segment%HI%ied
4560  j = segment%HI%JsdB
4561  ! jshift+J corresponds to the nearest interior tracer cell index
4562  ! jdir switches the sign of the flow so that positive is into the reservoir
4563  if (segment%direction == obc_direction_s) then
4564  jshift = 1 ; jdir = -1
4565  else
4566  jshift = 0 ; jdir = 1
4567  endif
4568  ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
4569  do m=1,ntr ; if (associated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
4570  v_l_out = max(0.0, (jdir*vhr(i,j,k))*segment%Tr_InvLscale_out / (h(i,j+jshift,k)*g%dxCv(i,j)))
4571  v_l_in = min(0.0, (jdir*vhr(i,j,k))*segment%Tr_InvLscale_in / (h(i,j+jshift,k)*g%dxCv(i,j)))
4572  fac1 = 1.0 + (v_l_out-v_l_in)
4573  segment%tr_Reg%Tr(m)%tres(i,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,j,k) + &
4574  (v_l_out*reg%Tr(m)%t(i,j+jshift,k) - &
4575  v_l_in*segment%tr_Reg%Tr(m)%t(i,j,k)))
4576  if (associated(obc%tres_y)) obc%tres_y(i,j,k,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
4577  enddo ; endif ; enddo
4578  enddo
4579  endif
4580  enddo ; endif ; endif
4581 

References obc_direction_s, and obc_direction_w.

Referenced by mom::step_mom_tracer_dyn().

Here is the caller graph for this function:

Variable Documentation

◆ id_clock_pass

integer mom_open_boundary::id_clock_pass

A CPU time clock.

Definition at line 301 of file MOM_open_boundary.F90.

301 integer :: id_clock_pass !< A CPU time clock

Referenced by open_boundary_init(), and radiation_open_bdry_conds().

◆ max_obc_fields

integer, parameter mom_open_boundary::max_obc_fields = 100
private

Maximum number of data fields needed for OBC segments.

Definition at line 67 of file MOM_open_boundary.F90.

67 integer, parameter :: MAX_OBC_FIELDS = 100 !< Maximum number of data fields needed for OBC segments

◆ mdl

character(len=40) mom_open_boundary::mdl = "MOM_open_boundary"
private

This module's name.

Definition at line 303 of file MOM_open_boundary.F90.

303 character(len=40) :: mdl = "MOM_open_boundary" !< This module's name.

Referenced by initialize_segment_data(), mask_outside_obcs(), open_boundary_config(), parse_for_tracer_reservoirs(), segment_tracer_registry_init(), setup_u_point_obc(), and setup_v_point_obc().

◆ obc_direction_e

integer, parameter, public mom_open_boundary::obc_direction_e = 300

◆ obc_direction_n

◆ obc_direction_s

◆ obc_direction_w

integer, parameter, public mom_open_boundary::obc_direction_w = 400

◆ obc_flather

integer, parameter, public mom_open_boundary::obc_flather = 3

Indicates the use of a Flather open boundary.

Definition at line 61 of file MOM_open_boundary.F90.

61 integer, parameter, public :: OBC_FLATHER = 3 !< Indicates the use of a Flather open boundary

◆ obc_none

integer, parameter, public mom_open_boundary::obc_none = 0

Indicates the use of no open boundary.

Definition at line 58 of file MOM_open_boundary.F90.

58 integer, parameter, public :: OBC_NONE = 0 !< Indicates the use of no open boundary

Referenced by open_boundary_config(), and open_boundary_impose_land_mask().

◆ obc_radiation

integer, parameter, public mom_open_boundary::obc_radiation = 4

Indicates the use of a radiation open boundary.

Definition at line 62 of file MOM_open_boundary.F90.

62 integer, parameter, public :: OBC_RADIATION = 4 !< Indicates the use of a radiation open boundary

◆ obc_simple

integer, parameter, public mom_open_boundary::obc_simple = 1

Indicates the use of a simple inflow open boundary.

Definition at line 59 of file MOM_open_boundary.F90.

59 integer, parameter, public :: OBC_SIMPLE = 1 !< Indicates the use of a simple inflow open boundary

◆ obc_wall

integer, parameter, public mom_open_boundary::obc_wall = 2

Indicates the use of a closed wall.

Definition at line 60 of file MOM_open_boundary.F90.

60 integer, parameter, public :: OBC_WALL = 2 !< Indicates the use of a closed wall
interpret_int_expr
integer function interpret_int_expr(string, imax)
Definition: MOM_open_boundary.F90:1214