9 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
12 use mom_domains,
only : to_all, scalar_pair, cgrid_ne, corner
18 use mom_io,
only : east_face, north_face
19 use mom_io,
only : slasher, read_data, field_size, single_file
23 use mom_obsolete_params,
only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char
26 use mom_time_manager,
only : set_date, time_type, time_type_to_real,
operator(-)
28 use time_interp_external_mod,
only : init_external_field, time_interp_external
29 use time_interp_external_mod,
only : time_interp_external_init
31 use mom_remapping,
only : initialize_remapping, remapping_core_h, end_remapping
37 implicit none ;
private 39 #include <MOM_memory.h> 41 public open_boundary_apply_normal_flow
42 public open_boundary_config
43 public open_boundary_init
44 public open_boundary_query
45 public open_boundary_end
46 public open_boundary_impose_normal_slope
47 public open_boundary_impose_land_mask
48 public radiation_open_bdry_conds
49 public set_tracer_data
50 public update_obc_segment_data
51 public open_boundary_test_extern_uv
52 public open_boundary_test_extern_h
53 public open_boundary_zero_normal_flow
54 public register_obc, obc_registry_init
55 public register_file_obc, file_obc_end
56 public segment_tracer_registry_init
57 public segment_tracer_registry_end
58 public register_segment_tracer
59 public register_temp_salt_segments
60 public fill_temp_salt_segments
61 public open_boundary_register_restarts
62 public update_segment_tracer_reservoirs
63 public update_obc_ramp
64 public rotate_obc_config
65 public rotate_obc_init
66 public initialize_segment_data
68 integer,
parameter,
public :: obc_none = 0
69 integer,
parameter,
public :: obc_simple = 1
70 integer,
parameter,
public :: obc_wall = 2
71 integer,
parameter,
public :: obc_flather = 3
72 integer,
parameter,
public :: obc_radiation = 4
73 integer,
parameter,
public :: obc_direction_n = 100
74 integer,
parameter,
public :: obc_direction_s = 200
75 integer,
parameter,
public :: obc_direction_e = 300
76 integer,
parameter,
public :: obc_direction_w = 400
77 integer,
parameter :: max_obc_fields = 100
83 character(len=8) :: name
84 real,
dimension(:,:,:),
allocatable :: buffer_src
87 real,
dimension(:,:,:),
allocatable :: dz_src
89 real,
dimension(:,:,:),
pointer :: buffer_dst=>null()
95 real,
dimension(:,:,:),
pointer :: t => null()
96 real :: obc_inflow_conc = 0.0
97 character(len=32) :: name
99 real,
dimension(:,:,:),
pointer :: tres => null()
100 logical :: is_initialized
107 logical :: locked = .false.
117 logical :: radiation_tan
119 logical :: radiation_grad
122 logical :: oblique_tan
124 logical :: oblique_grad
127 logical :: nudged_tan
128 logical :: nudged_grad
130 logical :: specified_tan
131 logical :: specified_grad
134 logical :: values_needed
135 logical :: u_values_needed
136 logical :: uamp_values_needed
137 logical :: uphase_values_needed
138 logical :: v_values_needed
139 logical :: vamp_values_needed
140 logical :: vphase_values_needed
141 logical :: t_values_needed
142 logical :: s_values_needed
143 logical :: z_values_needed
144 logical :: zamp_values_needed
145 logical :: zphase_values_needed
146 logical :: g_values_needed
150 logical :: is_e_or_w_2
152 integer :: num_fields
153 character(len=32),
pointer,
dimension(:) :: field_names=>null()
158 integer :: uamp_index
159 integer :: uphase_index
160 integer :: vamp_index
161 integer :: vphase_index
162 integer :: zamp_index
163 integer :: zphase_index
164 real :: velocity_nudging_timescale_in
165 real :: velocity_nudging_timescale_out
167 logical :: temp_segment_data_exists
168 logical :: salt_segment_data_exists
169 real,
pointer,
dimension(:,:) :: cg=>null()
171 real,
pointer,
dimension(:,:) :: htot=>null()
172 real,
pointer,
dimension(:,:,:) :: h=>null()
173 real,
pointer,
dimension(:,:,:) :: normal_vel=>null()
175 real,
pointer,
dimension(:,:,:) :: tangential_vel=>null()
177 real,
pointer,
dimension(:,:,:) :: tangential_grad=>null()
179 real,
pointer,
dimension(:,:,:) :: normal_trans=>null()
181 real,
pointer,
dimension(:,:) :: normal_vel_bt=>null()
183 real,
pointer,
dimension(:,:) :: eta=>null()
184 real,
pointer,
dimension(:,:,:) :: grad_normal=>null()
186 real,
pointer,
dimension(:,:,:) :: grad_tan=>null()
188 real,
pointer,
dimension(:,:,:) :: grad_gradient=>null()
190 real,
pointer,
dimension(:,:,:) :: rx_norm_rad=>null()
192 real,
pointer,
dimension(:,:,:) :: ry_norm_rad=>null()
194 real,
pointer,
dimension(:,:,:) :: rx_norm_obl=>null()
196 real,
pointer,
dimension(:,:,:) :: ry_norm_obl=>null()
198 real,
pointer,
dimension(:,:,:) :: cff_normal=>null()
200 real,
pointer,
dimension(:,:,:) :: nudged_normal_vel=>null()
202 real,
pointer,
dimension(:,:,:) :: nudged_tangential_vel=>null()
204 real,
pointer,
dimension(:,:,:) :: nudged_tangential_grad=>null()
207 type(hor_index_type) :: hi
208 real :: tr_invlscale_out
212 real :: tr_invlscale_in
219 integer :: number_of_segments = 0
221 logical :: open_u_bcs_exist_globally = .false.
223 logical :: open_v_bcs_exist_globally = .false.
225 logical :: flather_u_bcs_exist_globally = .false.
227 logical :: flather_v_bcs_exist_globally = .false.
229 logical :: oblique_bcs_exist_globally = .false.
231 logical :: nudged_u_bcs_exist_globally = .false.
233 logical :: nudged_v_bcs_exist_globally = .false.
235 logical :: specified_u_bcs_exist_globally = .false.
237 logical :: specified_v_bcs_exist_globally = .false.
239 logical :: radiation_bcs_exist_globally = .false.
240 logical :: user_bcs_set_globally = .false.
242 logical :: update_obc = .false.
243 logical :: needs_io_for_data = .false.
244 logical :: zero_vorticity = .false.
245 logical :: freeslip_vorticity = .false.
247 logical :: computed_vorticity = .false.
249 logical :: specified_vorticity = .false.
251 logical :: zero_strain = .false.
252 logical :: freeslip_strain = .false.
254 logical :: computed_strain = .false.
256 logical :: specified_strain = .false.
258 logical :: zero_biharmonic = .false.
260 logical :: brushcutter_mode = .false.
261 logical,
pointer,
dimension(:) :: &
262 tracer_x_reservoirs_used => null()
264 logical,
pointer,
dimension(:) :: &
265 tracer_y_reservoirs_used => null()
268 integer :: n_tide_constituents = 0
269 logical :: add_tide_constituents = .false.
271 character(len=2),
allocatable,
dimension(:) :: tide_names
272 real,
allocatable,
dimension(:) :: tide_frequencies
273 real,
allocatable,
dimension(:) :: tide_eq_phases
274 real,
allocatable,
dimension(:) :: tide_fn
275 real,
allocatable,
dimension(:) :: tide_un
276 logical :: add_eq_phase = .false.
278 logical :: add_nodal_terms = .false.
280 type(time_type) :: time_ref
286 integer,
pointer,
dimension(:,:) :: &
287 segnum_u => null(), &
300 real,
pointer,
dimension(:,:,:) :: &
301 rx_normal => null(), &
303 ry_normal => null(), &
305 rx_oblique => null(), &
306 ry_oblique => null(), &
308 real,
pointer,
dimension(:,:,:,:) :: &
315 logical :: ramp = .false.
317 logical :: ramping_is_activated = .false.
318 real :: ramp_timescale
319 real :: trunc_ramp_time
322 type(time_type) :: ramp_start_time
328 real :: tide_flow = 3.0e6
333 character(len=32) :: name
340 logical :: locked = .false.
344 integer :: id_clock_pass
346 character(len=40) :: mdl =
"MOM_open_boundary" 348 #include "version_variable.h" 359 subroutine open_boundary_config(G, US, param_file, OBC)
366 logical :: debug_obc, debug, mask_outside, reentrant_x, reentrant_y
367 character(len=15) :: segment_param_str
368 character(len=1024) :: segment_str
369 character(len=200) :: config1
370 real :: lscale_in, lscale_out
371 character(len=128) :: inputdir
372 logical :: answers_2018, default_2018_answers
373 logical :: check_reconstruction, check_remapping, force_bounds_in_subcell
374 character(len=32) :: remappingscheme
377 call get_param(param_file, mdl,
"OBC_NUMBER_OF_SEGMENTS", obc%number_of_segments, &
378 default=0, do_not_log=.true.)
380 "Controls where open boundaries are located, what kind of boundary condition "//&
381 "to impose, and what data to apply, if any.", &
382 all_default=(obc%number_of_segments<=0))
383 call get_param(param_file, mdl,
"OBC_NUMBER_OF_SEGMENTS", obc%number_of_segments, &
384 "The number of open boundary segments.", &
386 call get_param(param_file, mdl,
"OBC_USER_CONFIG", config1, &
387 "A string that sets how the open boundary conditions are "//&
388 " configured: \n", default=
"none", do_not_log=.true.)
389 call get_param(param_file, mdl,
"NK", obc%ke, &
390 "The number of model layers", default=0, do_not_log=.true.)
392 if (config1 /=
"none" .and. config1 /=
"dyed_obcs") obc%user_BCs_set_globally = .true.
394 if (obc%number_of_segments > 0)
then 395 call get_param(param_file, mdl,
"OBC_ZERO_VORTICITY", obc%zero_vorticity, &
396 "If true, sets relative vorticity to zero on open boundaries.", &
398 call get_param(param_file, mdl,
"OBC_FREESLIP_VORTICITY", obc%freeslip_vorticity, &
399 "If true, sets the normal gradient of tangential velocity to "//&
400 "zero in the relative vorticity on open boundaries. This cannot "//&
401 "be true if another OBC_XXX_VORTICITY option is True.", default=.true.)
402 call get_param(param_file, mdl,
"OBC_COMPUTED_VORTICITY", obc%computed_vorticity, &
403 "If true, uses the external values of tangential velocity "//&
404 "in the relative vorticity on open boundaries. This cannot "//&
405 "be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
406 call get_param(param_file, mdl,
"OBC_SPECIFIED_VORTICITY", obc%specified_vorticity, &
407 "If true, uses the external values of tangential velocity "//&
408 "in the relative vorticity on open boundaries. This cannot "//&
409 "be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
410 if ((obc%zero_vorticity .and. obc%freeslip_vorticity) .or. &
411 (obc%zero_vorticity .and. obc%computed_vorticity) .or. &
412 (obc%zero_vorticity .and. obc%specified_vorticity) .or. &
413 (obc%freeslip_vorticity .and. obc%computed_vorticity) .or. &
414 (obc%freeslip_vorticity .and. obc%specified_vorticity) .or. &
415 (obc%computed_vorticity .and. obc%specified_vorticity)) &
416 call mom_error(fatal,
"MOM_open_boundary.F90, open_boundary_config:\n"//&
417 "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"//&
418 "and OBC_IMPORTED_VORTICITY can be True at once.")
419 call get_param(param_file, mdl,
"OBC_ZERO_STRAIN", obc%zero_strain, &
420 "If true, sets the strain used in the stress tensor to zero on open boundaries.", &
422 call get_param(param_file, mdl,
"OBC_FREESLIP_STRAIN", obc%freeslip_strain, &
423 "If true, sets the normal gradient of tangential velocity to "//&
424 "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
425 "be true if another OBC_XXX_STRAIN option is True.", default=.true.)
426 call get_param(param_file, mdl,
"OBC_COMPUTED_STRAIN", obc%computed_strain, &
427 "If true, sets the normal gradient of tangential velocity to "//&
428 "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
429 "be true if another OBC_XXX_STRAIN option is True.", default=.false.)
430 call get_param(param_file, mdl,
"OBC_SPECIFIED_STRAIN", obc%specified_strain, &
431 "If true, sets the normal gradient of tangential velocity to "//&
432 "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
433 "be true if another OBC_XXX_STRAIN option is True.", default=.false.)
434 if ((obc%zero_strain .and. obc%freeslip_strain) .or. &
435 (obc%zero_strain .and. obc%computed_strain) .or. &
436 (obc%zero_strain .and. obc%specified_strain) .or. &
437 (obc%freeslip_strain .and. obc%computed_strain) .or. &
438 (obc%freeslip_strain .and. obc%specified_strain) .or. &
439 (obc%computed_strain .and. obc%specified_strain)) &
440 call mom_error(fatal,
"MOM_open_boundary.F90, open_boundary_config: \n"//&
441 "Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN, OBC_COMPUTED_STRAIN \n"//&
442 "and OBC_IMPORTED_STRAIN can be True at once.")
443 call get_param(param_file, mdl,
"OBC_ZERO_BIHARMONIC", obc%zero_biharmonic, &
444 "If true, zeros the Laplacian of flow on open boundaries in the biharmonic "//&
445 "viscosity term.", default=.false.)
446 call get_param(param_file, mdl,
"MASK_OUTSIDE_OBCS", mask_outside, &
447 "If true, set the areas outside open boundaries to be land.", &
449 call get_param(param_file, mdl,
"RAMP_OBCS", obc%ramp, &
450 "If true, ramps from zero to the external values over time, with"//&
451 "a ramping timescale given by RAMP_TIMESCALE. Ramping SSH only so far", &
453 call get_param(param_file, mdl,
"OBC_RAMP_TIMESCALE", obc%ramp_timescale, &
454 "If RAMP_OBCS is true, this sets the ramping timescale.", &
455 units=
"days", default=1.0, scale=86400.0*us%s_to_T)
456 call get_param(param_file, mdl,
"OBC_TIDE_N_CONSTITUENTS", obc%n_tide_constituents, &
457 "Number of tidal constituents being added to the open boundary.", &
460 if (obc%n_tide_constituents > 0)
then 461 obc%add_tide_constituents = .true.
463 obc%add_tide_constituents = .false.
466 call get_param(param_file, mdl,
"DEBUG", debug, default=.false.)
467 call get_param(param_file, mdl,
"DEBUG_OBC", debug_obc, default=.false.)
468 if (debug_obc .or. debug) &
469 call log_param(param_file, mdl,
"DEBUG_OBC", debug_obc, &
470 "If true, do additional calls to help debug the performance "//&
471 "of the open boundary condition code.", default=.false., &
472 debuggingparam=.true.)
474 call get_param(param_file, mdl,
"OBC_SILLY_THICK", obc%silly_h, &
475 "A silly value of thicknesses used outside of open boundary "//&
476 "conditions for debugging.", units=
"m", default=0.0, scale=us%m_to_Z, &
477 do_not_log=.not.debug_obc, debuggingparam=.true.)
478 call get_param(param_file, mdl,
"OBC_SILLY_VEL", obc%silly_u, &
479 "A silly value of velocities used outside of open boundary "//&
480 "conditions for debugging.", units=
"m/s", default=0.0, scale=us%m_s_to_L_T, &
481 do_not_log=.not.debug_obc, debuggingparam=.true.)
482 reentrant_x = .false.
483 call get_param(param_file, mdl,
"REENTRANT_X", reentrant_x, default=.true.)
484 reentrant_y = .false.
485 call get_param(param_file, mdl,
"REENTRANT_Y", reentrant_y, default=.false.)
488 allocate(obc%segment(1:obc%number_of_segments))
489 do l=1,obc%number_of_segments
490 obc%segment(l)%Flather = .false.
491 obc%segment(l)%radiation = .false.
492 obc%segment(l)%radiation_tan = .false.
493 obc%segment(l)%radiation_grad = .false.
494 obc%segment(l)%oblique = .false.
495 obc%segment(l)%oblique_tan = .false.
496 obc%segment(l)%oblique_grad = .false.
497 obc%segment(l)%nudged = .false.
498 obc%segment(l)%nudged_tan = .false.
499 obc%segment(l)%nudged_grad = .false.
500 obc%segment(l)%specified = .false.
501 obc%segment(l)%specified_tan = .false.
502 obc%segment(l)%specified_grad = .false.
503 obc%segment(l)%open = .false.
504 obc%segment(l)%gradient = .false.
505 obc%segment(l)%values_needed = .false.
506 obc%segment(l)%u_values_needed = .false.
507 obc%segment(l)%uamp_values_needed = obc%add_tide_constituents
508 obc%segment(l)%uphase_values_needed = obc%add_tide_constituents
509 obc%segment(l)%v_values_needed = .false.
510 obc%segment(l)%vamp_values_needed = obc%add_tide_constituents
511 obc%segment(l)%vphase_values_needed = obc%add_tide_constituents
512 obc%segment(l)%t_values_needed = .false.
513 obc%segment(l)%s_values_needed = .false.
514 obc%segment(l)%z_values_needed = .false.
515 obc%segment(l)%zamp_values_needed = obc%add_tide_constituents
516 obc%segment(l)%zphase_values_needed = obc%add_tide_constituents
517 obc%segment(l)%g_values_needed = .false.
518 obc%segment(l)%direction = obc_none
519 obc%segment(l)%is_N_or_S = .false.
520 obc%segment(l)%is_E_or_W = .false.
521 obc%segment(l)%is_E_or_W_2 = .false.
522 obc%segment(l)%Velocity_nudging_timescale_in = 0.0
523 obc%segment(l)%Velocity_nudging_timescale_out = 0.0
524 obc%segment(l)%num_fields = 0
526 allocate(obc%segnum_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; obc%segnum_u(:,:) = obc_none
527 allocate(obc%segnum_v(g%isd:g%ied,g%JsdB:g%JedB)) ; obc%segnum_v(:,:) = obc_none
529 do l = 1, obc%number_of_segments
530 write(segment_param_str(1:15),
"('OBC_SEGMENT_',i3.3)") l
531 call get_param(param_file, mdl, segment_param_str, segment_str, &
532 "Documentation needs to be dynamic?????", &
533 fail_if_missing=.true.)
534 segment_str = remove_spaces(segment_str)
535 if (segment_str(1:2) ==
'I=')
then 536 call setup_u_point_obc(obc, g, us, segment_str, l, param_file, reentrant_y)
537 elseif (segment_str(1:2) ==
'J=')
then 538 call setup_v_point_obc(obc, g, us, segment_str, l, param_file, reentrant_x)
540 call mom_error(fatal,
"MOM_open_boundary.F90, open_boundary_config: "//&
541 "Unable to interpret "//segment_param_str//
" = "//trim(segment_str))
547 if (obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
548 obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)
then 550 call time_interp_external_init()
555 if (open_boundary_query(obc, apply_open_obc=.true.))
then 556 call get_param(param_file, mdl,
"OBC_RADIATION_MAX", obc%rx_max, &
557 "The maximum magnitude of the baroclinic radiation velocity (or speed of "//&
558 "characteristics), in gridpoints per timestep. This is only "//&
559 "used if one of the open boundary segments is using Orlanski.", &
560 units=
"nondim", default=1.0)
561 call get_param(param_file, mdl,
"OBC_RAD_VEL_WT", obc%gamma_uv, &
562 "The relative weighting for the baroclinic radiation "//&
563 "velocities (or speed of characteristics) at the new "//&
564 "time level (1) or the running mean (0) for velocities. "//&
565 "Valid values range from 0 to 1. This is only used if "//&
566 "one of the open boundary segments is using Orlanski.", &
567 units=
"nondim", default=0.3)
572 if (open_boundary_query(obc, apply_open_obc=.true.))
then 573 call get_param(param_file, mdl,
"OBC_TRACER_RESERVOIR_LENGTH_SCALE_OUT ", lscale_out, &
574 "An effective length scale for restoring the tracer concentration "//&
575 "at the boundaries to externally imposed values when the flow "//&
576 "is exiting the domain.", units=
"m", default=0.0, scale=us%m_to_L)
578 call get_param(param_file, mdl,
"OBC_TRACER_RESERVOIR_LENGTH_SCALE_IN ", lscale_in, &
579 "An effective length scale for restoring the tracer concentration "//&
580 "at the boundaries to values from the interior when the flow "//&
581 "is entering the domain.", units=
"m", default=0.0, scale=us%m_to_L)
584 if (mask_outside)
call mask_outside_obcs(g, us, param_file, obc)
589 do l = 1, obc%number_of_segments
590 obc%segment(l)%Tr_InvLscale_in = 0.0
591 if (lscale_in>0.) obc%segment(l)%Tr_InvLscale_in = 1.0/lscale_in
592 obc%segment(l)%Tr_InvLscale_out = 0.0
593 if (lscale_out>0.) obc%segment(l)%Tr_InvLscale_out = 1.0/lscale_out
596 call get_param(param_file, mdl,
"REMAPPING_SCHEME", remappingscheme, &
597 "This sets the reconstruction scheme used "//&
598 "for vertical remapping for all variables. "//&
599 "It can be one of the following schemes: \n"//&
600 trim(remappingschemesdoc), default=remappingdefaultscheme,do_not_log=.true.)
601 call get_param(param_file, mdl,
"FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
602 "If true, cell-by-cell reconstructions are checked for "//&
603 "consistency and if non-monotonicity or an inconsistency is "//&
604 "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
605 call get_param(param_file, mdl,
"FATAL_CHECK_REMAPPING", check_remapping, &
606 "If true, the results of remapping are checked for "//&
607 "conservation and new extrema and if an inconsistency is "//&
608 "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
609 call get_param(param_file, mdl,
"BRUSHCUTTER_MODE", obc%brushcutter_mode, &
610 "If true, read external OBC data on the supergrid.", &
612 call get_param(param_file, mdl,
"REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
613 "If true, the values on the intermediate grid used for remapping "//&
614 "are forced to be bounded, which might not be the case due to "//&
615 "round off.", default=.false.,do_not_log=.true.)
616 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
617 "This sets the default value for the various _2018_ANSWERS parameters.", &
619 call get_param(param_file, mdl,
"REMAPPING_2018_ANSWERS", answers_2018, &
620 "If true, use the order of arithmetic and expressions that recover the "//&
621 "answers from the end of 2018. Otherwise, use updated and more robust "//&
622 "forms of the same expressions.", default=default_2018_answers)
624 allocate(obc%remap_CS)
625 call initialize_remapping(obc%remap_CS, remappingscheme, boundary_extrapolation = .false., &
626 check_reconstruction=check_reconstruction, check_remapping=check_remapping, &
627 force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018)
632 if ((obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally) .and. &
633 .not.g%symmetric )
call mom_error(fatal, &
634 "MOM_open_boundary, open_boundary_config: "//&
635 "Symmetric memory must be used when using Flather OBCs.")
637 if (obc%add_tide_constituents)
then 638 call initialize_obc_tides(obc, param_file)
640 obc%update_OBC = .true.
643 if (.not.(obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
644 obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally))
then 646 call open_boundary_dealloc(obc)
649 end subroutine open_boundary_config
653 subroutine initialize_segment_data(G, OBC, PF)
654 use mpp_mod,
only : mpp_pe, mpp_set_current_pelist, mpp_get_current_pelist,mpp_npes
660 integer :: n,m,num_fields
661 character(len=1024) :: segstr
662 character(len=256) :: filename
663 character(len=20) :: segnam, suffix
664 character(len=32) :: varnam, fieldname
666 character(len=32),
dimension(MAX_OBC_FIELDS) :: fields
667 character(len=128) :: inputdir
669 character(len=256) :: mesg
670 integer,
dimension(4) :: siz,siz2
671 integer :: is, ie, js, je
672 integer :: isd, ied, jsd, jed
673 integer :: isdb, iedb, jsdb, jedb
674 integer,
dimension(:),
allocatable :: saved_pelist
675 integer :: current_pe
676 integer,
dimension(1) :: single_pelist
679 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
686 call get_param(pf, mdl,
"INPUTDIR", inputdir, default=
".")
687 inputdir = slasher(inputdir)
689 if (obc%user_BCs_set_globally)
return 692 do n=1, obc%number_of_segments
693 segment => obc%segment(n)
694 write(segnam,
"('OBC_SEGMENT_',i3.3,'_DATA')") n
695 call get_param(pf, mdl, segnam, segstr,
'OBC segment docs')
700 allocate(saved_pelist(0:mpp_npes()-1))
701 call mpp_get_current_pelist(saved_pelist)
702 current_pe = mpp_pe()
703 single_pelist(1) = current_pe
704 call mpp_set_current_pelist(single_pelist)
706 do n=1, obc%number_of_segments
707 segment => obc%segment(n)
708 if (.not. segment%values_needed) cycle
710 write(segnam,
"('OBC_SEGMENT_',i3.3,'_DATA')") n
711 write(suffix,
"('_segment_',i3.3)") n
718 if (segstr ==
'')
then 719 write(mesg,
'("No OBC_SEGMENT_XXX_DATA string for OBC segment ",I3)') n
720 call mom_error(fatal, mesg)
723 call parse_segment_manifest_str(trim(segstr), num_fields, fields)
724 if (num_fields == 0)
then 725 call mom_mesg(
'initialize_segment_data: num_fields = 0')
729 allocate(segment%field(num_fields))
730 segment%num_fields = num_fields
732 segment%temp_segment_data_exists=.false.
733 segment%salt_segment_data_exists=.false.
738 isd = segment%HI%isd ; ied = segment%HI%ied
739 jsd = segment%HI%jsd ; jed = segment%HI%jed
740 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
741 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
744 call parse_segment_data_str(trim(segstr), m, trim(fields(m)), &
745 value, filename, fieldname)
746 if (trim(filename) /=
'none')
then 747 obc%update_OBC = .true.
748 obc%needs_IO_for_data = .true.
750 segment%field(m)%name = trim(fields(m))
751 if (segment%field(m)%name ==
'TEMP')
then 752 segment%temp_segment_data_exists=.true.
753 segment%t_values_needed = .false.
755 if (segment%field(m)%name ==
'SALT')
then 756 segment%salt_segment_data_exists=.true.
757 segment%s_values_needed = .false.
759 filename = trim(inputdir)//trim(filename)
760 fieldname = trim(fieldname)//trim(suffix)
761 call field_size(filename,fieldname,siz,no_domain=.true.)
763 if (segment%on_pe)
then 764 if (obc%brushcutter_mode .and. (modulo(siz(1),2) == 0 .or. modulo(siz(2),2) == 0))
then 765 call mom_error(fatal,
'segment data are not on the supergrid')
770 if (obc%brushcutter_mode)
then 778 if (obc%brushcutter_mode)
then 786 if (segment%is_E_or_W)
then 787 if (segment%field(m)%name ==
'V')
then 788 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
789 segment%v_values_needed = .false.
790 elseif (segment%field(m)%name ==
'Vamp')
then 791 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
792 segment%vamp_values_needed = .false.
793 segment%vamp_index = m
794 elseif (segment%field(m)%name ==
'Vphase')
then 795 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
796 segment%vphase_values_needed = .false.
797 segment%vphase_index = m
798 elseif (segment%field(m)%name ==
'DVDX')
then 799 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
800 segment%g_values_needed = .false.
802 allocate(segment%field(m)%buffer_src(isdb:iedb,jsd:jed,siz2(3)))
803 if (segment%field(m)%name ==
'U')
then 804 segment%u_values_needed = .false.
805 elseif (segment%field(m)%name ==
'Uamp')
then 806 segment%uamp_values_needed = .false.
807 segment%uamp_index = m
808 elseif (segment%field(m)%name ==
'Uphase')
then 809 segment%uphase_values_needed = .false.
810 segment%uphase_index = m
811 elseif (segment%field(m)%name ==
'SSH')
then 812 segment%z_values_needed = .false.
813 elseif (segment%field(m)%name ==
'SSHamp')
then 814 segment%zamp_values_needed = .false.
815 segment%zamp_index = m
816 elseif (segment%field(m)%name ==
'SSHphase')
then 817 segment%zphase_values_needed = .false.
818 segment%zphase_index = m
819 elseif (segment%field(m)%name ==
'TEMP')
then 820 segment%t_values_needed = .false.
821 elseif (segment%field(m)%name ==
'SALT')
then 822 segment%s_values_needed = .false.
826 if (segment%field(m)%name ==
'U')
then 827 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
828 segment%u_values_needed = .false.
829 elseif (segment%field(m)%name ==
'DUDY')
then 830 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
831 segment%g_values_needed = .false.
832 elseif (segment%field(m)%name ==
'Uamp')
then 833 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
834 segment%uamp_values_needed = .false.
835 segment%uamp_index = m
836 elseif (segment%field(m)%name ==
'Uphase')
then 837 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
838 segment%uphase_values_needed = .false.
839 segment%uphase_index = m
841 allocate(segment%field(m)%buffer_src(isd:ied,jsdb:jedb,siz2(3)))
842 if (segment%field(m)%name ==
'V')
then 843 segment%v_values_needed = .false.
844 elseif (segment%field(m)%name ==
'Vamp')
then 845 segment%vamp_values_needed = .false.
846 segment%vamp_index = m
847 elseif (segment%field(m)%name ==
'Vphase')
then 848 segment%vphase_values_needed = .false.
849 segment%vphase_index = m
850 elseif (segment%field(m)%name ==
'SSH')
then 851 segment%z_values_needed = .false.
852 elseif (segment%field(m)%name ==
'SSHamp')
then 853 segment%zamp_values_needed = .false.
854 segment%zamp_index = m
855 elseif (segment%field(m)%name ==
'SSHphase')
then 856 segment%zphase_values_needed = .false.
857 segment%zphase_index = m
858 elseif (segment%field(m)%name ==
'TEMP')
then 859 segment%t_values_needed = .false.
860 elseif (segment%field(m)%name ==
'SALT')
then 861 segment%s_values_needed = .false.
865 segment%field(m)%buffer_src(:,:,:)=0.0
866 segment%field(m)%fid = init_external_field(trim(filename),&
867 trim(fieldname),ignore_axis_atts=.true.,threading=single_file)
869 if ((index(segment%field(m)%name,
'phase') > 0) .or. (index(segment%field(m)%name,
'amp') > 0))
then 871 call field_size(filename,
'constituent', siz, no_domain=.true.)
873 if (siz(3) .ne. obc%n_tide_constituents .and. obc%add_tide_constituents)
then 874 call mom_error(fatal,
'Number of constituents in input data is not '//&
875 'the same as the number specified')
877 segment%field(m)%nk_src=siz(3)
880 fieldname =
'dz_'//trim(fieldname)
881 call field_size(filename,fieldname,siz,no_domain=.true.)
882 if (segment%is_E_or_W)
then 883 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then 884 allocate(segment%field(m)%dz_src(isdb:iedb,jsdb:jedb,siz(3)))
886 allocate(segment%field(m)%dz_src(isdb:iedb,jsd:jed,siz(3)))
889 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then 890 allocate(segment%field(m)%dz_src(isdb:iedb,jsdb:jedb,siz(3)))
892 allocate(segment%field(m)%dz_src(isd:ied,jsdb:jedb,siz(3)))
895 segment%field(m)%dz_src(:,:,:)=0.0
896 segment%field(m)%nk_src=siz(3)
897 segment%field(m)%fid_dz = init_external_field(trim(filename),trim(fieldname),&
898 ignore_axis_atts=.true.,threading=single_file)
901 segment%field(m)%nk_src=1
905 segment%field(m)%fid = -1
906 segment%field(m)%value =
value 907 segment%field(m)%name = trim(fields(m))
910 if ((index(segment%field(m)%name,
'phase') > 0) .or. (index(segment%field(m)%name,
'amp') > 0))
then 911 if (obc%n_tide_constituents .gt. 1 .and. obc%add_tide_constituents)
then 912 call mom_error(fatal,
'Only one constituent is supported when specifying '//&
913 'tidal boundary conditions by value rather than file.')
916 if (segment%field(m)%name ==
'U')
then 917 segment%u_values_needed = .false.
918 elseif (segment%field(m)%name ==
'Uamp')
then 919 segment%uamp_values_needed = .false.
920 segment%uamp_index = m
921 elseif (segment%field(m)%name ==
'Uphase')
then 922 segment%uphase_values_needed = .false.
923 segment%uphase_index = m
924 elseif (segment%field(m)%name ==
'V')
then 925 segment%v_values_needed = .false.
926 elseif (segment%field(m)%name ==
'Vamp')
then 927 segment%vamp_values_needed = .false.
928 segment%vamp_index = m
929 elseif (segment%field(m)%name ==
'Vphase')
then 930 segment%vphase_values_needed = .false.
931 segment%vphase_index = m
932 elseif (segment%field(m)%name ==
'SSH')
then 933 segment%z_values_needed = .false.
934 elseif (segment%field(m)%name ==
'SSHamp')
then 935 segment%zamp_values_needed = .false.
936 segment%zamp_index = m
937 elseif (segment%field(m)%name ==
'SSHphase')
then 938 segment%zphase_values_needed = .false.
939 segment%zphase_index = m
940 elseif (segment%field(m)%name ==
'TEMP')
then 941 segment%t_values_needed = .false.
942 elseif (segment%field(m)%name ==
'SALT')
then 943 segment%s_values_needed = .false.
944 elseif (segment%field(m)%name ==
'DVDX' .or. segment%field(m)%name ==
'DUDY')
then 945 segment%g_values_needed = .false.
949 if (segment%u_values_needed .or. segment%uamp_values_needed .or. segment%uphase_values_needed .or. &
950 segment%v_values_needed .or. segment%vamp_values_needed .or. segment%vphase_values_needed .or. &
951 segment%t_values_needed .or. segment%s_values_needed .or. segment%g_values_needed .or. &
952 segment%z_values_needed .or. segment%zamp_values_needed .or. segment%zphase_values_needed )
then 953 write(mesg,
'("Values needed for OBC segment ",I3)') n
954 call mom_error(fatal, mesg)
958 call mpp_set_current_pelist(saved_pelist)
960 end subroutine initialize_segment_data
962 subroutine initialize_obc_tides(OBC, param_file)
965 integer,
dimension(3) :: tide_ref_date
966 integer,
dimension(3) :: nodal_ref_date
967 character(len=50) :: tide_constituent_str
969 type(time_type) :: nodal_time
972 call get_param(param_file, mdl,
"OBC_TIDE_CONSTITUENTS", tide_constituent_str, &
973 "Names of tidal constituents being added to the open boundaries.", &
974 fail_if_missing=.true.)
976 call get_param(param_file, mdl,
"OBC_TIDE_ADD_EQ_PHASE", obc%add_eq_phase, &
977 "If true, add the equilibrium phase argument to the specified tidal phases.", &
978 default=.false., fail_if_missing=.false.)
980 call get_param(param_file, mdl,
"OBC_TIDE_ADD_NODAL", obc%add_nodal_terms, &
981 "If true, include 18.6 year nodal modulation in the boundary tidal forcing.", &
984 call get_param(param_file, mdl,
"OBC_TIDE_REF_DATE", tide_ref_date, &
985 "Reference date to use for tidal calculations and equilibrium phase.", &
986 fail_if_missing=.true.)
988 call get_param(param_file, mdl,
"OBC_TIDE_NODAL_REF_DATE", nodal_ref_date, &
989 "Fixed reference date to use for nodal modulation of boundary tides.", &
990 fail_if_missing=.false., default=0)
992 if (.not. obc%add_eq_phase)
then 995 call mom_mesg(
'OBC tidal phases will *not* be corrected with equilibrium arguments.')
998 allocate(obc%tide_names(obc%n_tide_constituents))
999 read(tide_constituent_str, *) obc%tide_names
1002 obc%time_ref = set_date(tide_ref_date(1), tide_ref_date(2), tide_ref_date(3))
1005 if (obc%add_eq_phase)
call astro_longitudes_init(obc%time_ref, obc%tidal_longitudes)
1009 if (obc%add_nodal_terms)
then 1010 if (sum(nodal_ref_date) .ne. 0)
then 1012 nodal_time = set_date(nodal_ref_date(1), nodal_ref_date(2), nodal_ref_date(3))
1013 call astro_longitudes_init(nodal_time, nodal_longitudes)
1014 elseif (obc%add_eq_phase)
then 1017 nodal_longitudes = obc%tidal_longitudes
1020 call astro_longitudes_init(obc%time_ref, nodal_longitudes)
1024 allocate(obc%tide_frequencies(obc%n_tide_constituents))
1025 allocate(obc%tide_eq_phases(obc%n_tide_constituents))
1026 allocate(obc%tide_fn(obc%n_tide_constituents))
1027 allocate(obc%tide_un(obc%n_tide_constituents))
1029 do c=1,obc%n_tide_constituents
1032 call get_param(param_file, mdl,
"TIDE_"//trim(obc%tide_names(c))//
"_FREQ", obc%tide_frequencies(c), &
1033 "Frequency of the "//trim(obc%tide_names(c))//
" tidal constituent. "//&
1034 "This is only used if TIDES and TIDE_"//trim(obc%tide_names(c))// &
1035 " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(obc%tide_names(c))//&
1036 " is in OBC_TIDE_CONSTITUENTS.", units=
"s-1", default=tidal_frequency(trim(obc%tide_names(c))))
1039 if (obc%add_eq_phase)
then 1040 obc%tide_eq_phases(c) = eq_phase(trim(obc%tide_names(c)), obc%tidal_longitudes)
1042 obc%tide_eq_phases(c) = 0.0
1046 if (obc%add_nodal_terms)
then 1047 call nodal_fu(trim(obc%tide_names(c)), nodal_longitudes%N, obc%tide_fn(c), obc%tide_un(c))
1049 obc%tide_fn(c) = 1.0
1050 obc%tide_un(c) = 0.0
1053 end subroutine initialize_obc_tides
1057 subroutine setup_segment_indices(G, seg, Is_obc, Ie_obc, Js_obc, Je_obc)
1060 integer,
intent(in) :: Is_obc
1061 integer,
intent(in) :: Ie_obc
1062 integer,
intent(in) :: Js_obc
1063 integer,
intent(in) :: Je_obc
1065 integer :: IsgB, IegB, JsgB, JegB
1066 integer :: isg, ieg, jsg, jeg
1069 if (ie_obc < is_obc)
then 1077 if (je_obc < js_obc)
then 1100 if (is_obc > ie_obc)
then 1108 if (is_obc < ie_obc)
then 1116 if (js_obc < je_obc)
then 1124 if (js_obc > je_obc)
then 1144 isgb = isgb - g%idg_offset
1145 jsgb = jsgb - g%jdg_offset
1146 iegb = iegb - g%idg_offset
1147 jegb = jegb - g%jdg_offset
1149 isg = isg - g%idg_offset
1150 jsg = jsg - g%jdg_offset
1151 ieg = ieg - g%idg_offset
1152 jeg = jeg - g%jdg_offset
1156 seg%HI%IsdB = min(max(isgb, g%HI%IsdB), g%HI%IedB)
1157 seg%HI%IedB = min(max(iegb, g%HI%IsdB), g%HI%IedB)
1158 seg%HI%isd = min(max(isg, g%HI%isd), g%HI%ied)
1159 seg%HI%ied = min(max(ieg, g%HI%isd), g%HI%ied)
1160 seg%HI%IscB = min(max(isgb, g%HI%IscB), g%HI%IecB)
1161 seg%HI%IecB = min(max(iegb, g%HI%IscB), g%HI%IecB)
1162 seg%HI%isc = min(max(isg, g%HI%isc), g%HI%iec)
1163 seg%HI%iec = min(max(ieg, g%HI%isc), g%HI%iec)
1167 seg%HI%JsdB = min(max(jsgb, g%HI%JsdB), g%HI%JedB)
1168 seg%HI%JedB = min(max(jegb, g%HI%JsdB), g%HI%JedB)
1169 seg%HI%jsd = min(max(jsg, g%HI%jsd), g%HI%jed)
1170 seg%HI%jed = min(max(jeg, g%HI%jsd), g%HI%jed)
1171 seg%HI%JscB = min(max(jsgb, g%HI%JscB), g%HI%JecB)
1172 seg%HI%JecB = min(max(jegb, g%HI%JscB), g%HI%JecB)
1173 seg%HI%jsc = min(max(jsg, g%HI%jsc), g%HI%jec)
1174 seg%HI%jec = min(max(jeg, g%HI%jsc), g%HI%jec)
1176 end subroutine setup_segment_indices
1179 subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y)
1183 character(len=*),
intent(in) :: segment_str
1184 integer,
intent(in) :: l_seg
1186 logical,
intent(in) :: reentrant_y
1188 integer :: I_obc, Js_obc, Je_obc
1189 integer :: j, a_loop
1190 character(len=32) :: action_str(8)
1191 character(len=128) :: segment_param_str
1192 real,
allocatable,
dimension(:) :: tnudge
1194 call parse_segment_str(g%ieg, g%jeg, segment_str, i_obc, js_obc, je_obc, action_str, reentrant_y)
1196 call setup_segment_indices(g, obc%segment(l_seg),i_obc,i_obc,js_obc,je_obc)
1198 i_obc = i_obc - g%idg_offset
1199 js_obc = js_obc - g%jdg_offset
1200 je_obc = je_obc - g%jdg_offset
1202 if (je_obc>js_obc)
then 1203 obc%segment(l_seg)%direction = obc_direction_e
1204 elseif (je_obc<js_obc)
then 1205 obc%segment(l_seg)%direction = obc_direction_w
1206 j=js_obc;js_obc=je_obc;je_obc=j
1209 obc%segment(l_seg)%on_pe = .false.
1212 if (len_trim(action_str(a_loop)) == 0)
then 1214 elseif (trim(action_str(a_loop)) ==
'FLATHER')
then 1215 obc%segment(l_seg)%Flather = .true.
1216 obc%segment(l_seg)%open = .true.
1217 obc%Flather_u_BCs_exist_globally = .true.
1218 obc%open_u_BCs_exist_globally = .true.
1219 obc%segment(l_seg)%z_values_needed = .true.
1220 obc%segment(l_seg)%u_values_needed = .true.
1221 elseif (trim(action_str(a_loop)) ==
'ORLANSKI')
then 1222 obc%segment(l_seg)%radiation = .true.
1223 obc%segment(l_seg)%open = .true.
1224 obc%open_u_BCs_exist_globally = .true.
1225 obc%radiation_BCs_exist_globally = .true.
1226 elseif (trim(action_str(a_loop)) ==
'ORLANSKI_TAN')
then 1227 obc%segment(l_seg)%radiation = .true.
1228 obc%segment(l_seg)%radiation_tan = .true.
1229 obc%radiation_BCs_exist_globally = .true.
1230 elseif (trim(action_str(a_loop)) ==
'ORLANSKI_GRAD')
then 1231 obc%segment(l_seg)%radiation = .true.
1232 obc%segment(l_seg)%radiation_grad = .true.
1233 elseif (trim(action_str(a_loop)) ==
'OBLIQUE')
then 1234 obc%segment(l_seg)%oblique = .true.
1235 obc%segment(l_seg)%open = .true.
1236 obc%oblique_BCs_exist_globally = .true.
1237 obc%open_u_BCs_exist_globally = .true.
1238 elseif (trim(action_str(a_loop)) ==
'OBLIQUE_TAN')
then 1239 obc%segment(l_seg)%oblique = .true.
1240 obc%segment(l_seg)%oblique_tan = .true.
1241 obc%oblique_BCs_exist_globally = .true.
1242 elseif (trim(action_str(a_loop)) ==
'OBLIQUE_GRAD')
then 1243 obc%segment(l_seg)%oblique = .true.
1244 obc%segment(l_seg)%oblique_grad = .true.
1245 elseif (trim(action_str(a_loop)) ==
'NUDGED')
then 1246 obc%segment(l_seg)%nudged = .true.
1247 obc%nudged_u_BCs_exist_globally = .true.
1248 obc%segment(l_seg)%u_values_needed = .true.
1249 elseif (trim(action_str(a_loop)) ==
'NUDGED_TAN')
then 1250 obc%segment(l_seg)%nudged_tan = .true.
1251 obc%nudged_u_BCs_exist_globally = .true.
1252 obc%segment(l_seg)%v_values_needed = .true.
1253 elseif (trim(action_str(a_loop)) ==
'NUDGED_GRAD')
then 1254 obc%segment(l_seg)%nudged_grad = .true.
1255 obc%segment(l_seg)%g_values_needed = .true.
1256 elseif (trim(action_str(a_loop)) ==
'GRADIENT')
then 1257 obc%segment(l_seg)%gradient = .true.
1258 obc%segment(l_seg)%open = .true.
1259 obc%open_u_BCs_exist_globally = .true.
1260 elseif (trim(action_str(a_loop)) ==
'SIMPLE')
then 1261 obc%segment(l_seg)%specified = .true.
1262 obc%specified_u_BCs_exist_globally = .true.
1263 obc%segment(l_seg)%u_values_needed = .true.
1264 elseif (trim(action_str(a_loop)) ==
'SIMPLE_TAN')
then 1265 obc%segment(l_seg)%specified_tan = .true.
1266 obc%segment(l_seg)%v_values_needed = .true.
1267 elseif (trim(action_str(a_loop)) ==
'SIMPLE_GRAD')
then 1268 obc%segment(l_seg)%specified_grad = .true.
1269 obc%segment(l_seg)%g_values_needed = .true.
1271 call mom_error(fatal,
"MOM_open_boundary.F90, setup_u_point_obc: "//&
1272 "String '"//trim(action_str(a_loop))//
"' not understood.")
1274 if (obc%segment(l_seg)%nudged .or. obc%segment(l_seg)%nudged_tan)
then 1275 write(segment_param_str(1:43),
"('OBC_SEGMENT_',i3.3,'_VELOCITY_NUDGING_TIMESCALES')") l_seg
1277 call get_param(pf, mdl, segment_param_str(1:43), tnudge, &
1278 "Timescales in days for nudging along a segment, "//&
1279 "for inflow, then outflow. Setting both to zero should "//&
1280 "behave like SIMPLE obcs for the baroclinic velocities.", &
1281 fail_if_missing=.true., default=0., units=
"days", scale=86400.0*us%s_to_T)
1282 obc%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1)
1283 obc%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2)
1289 obc%segment(l_seg)%is_E_or_W_2 = .true.
1291 if (i_obc<=g%HI%IsdB+1 .or. i_obc>=g%HI%IedB-1)
return 1292 if (je_obc<=g%HI%JsdB .or. js_obc>=g%HI%JedB)
return 1294 obc%segment(l_seg)%on_pe = .true.
1295 obc%segment(l_seg)%is_E_or_W = .true.
1297 do j=g%HI%jsd, g%HI%jed
1298 if (j>js_obc .and. j<=je_obc)
then 1299 obc%segnum_u(i_obc,j) = l_seg
1302 obc%segment(l_seg)%Is_obc = i_obc
1303 obc%segment(l_seg)%Ie_obc = i_obc
1304 obc%segment(l_seg)%Js_obc = js_obc
1305 obc%segment(l_seg)%Je_obc = je_obc
1306 call allocate_obc_segment_data(obc, obc%segment(l_seg))
1308 if (obc%segment(l_seg)%oblique .and. obc%segment(l_seg)%radiation) &
1309 call mom_error(fatal,
"MOM_open_boundary.F90, setup_u_point_obc: \n"//&
1310 "Orlanski and Oblique OBC options cannot be used together on one segment.")
1312 if (obc%segment(l_seg)%u_values_needed .or. obc%segment(l_seg)%v_values_needed .or. &
1313 obc%segment(l_seg)%t_values_needed .or. obc%segment(l_seg)%s_values_needed .or. &
1314 obc%segment(l_seg)%z_values_needed .or. obc%segment(l_seg)%g_values_needed) &
1315 obc%segment(l_seg)%values_needed = .true.
1316 end subroutine setup_u_point_obc
1319 subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x)
1323 character(len=*),
intent(in) :: segment_str
1324 integer,
intent(in) :: l_seg
1326 logical,
intent(in) :: reentrant_x
1328 integer :: J_obc, Is_obc, Ie_obc
1329 integer :: i, a_loop
1330 character(len=32) :: action_str(8)
1331 character(len=128) :: segment_param_str
1332 real,
allocatable,
dimension(:) :: tnudge
1335 call parse_segment_str(g%ieg, g%jeg, segment_str, j_obc, is_obc, ie_obc, action_str, reentrant_x)
1337 call setup_segment_indices(g, obc%segment(l_seg),is_obc,ie_obc,j_obc,j_obc)
1339 j_obc = j_obc - g%jdg_offset
1340 is_obc = is_obc - g%idg_offset
1341 ie_obc = ie_obc - g%idg_offset
1343 if (ie_obc>is_obc)
then 1344 obc%segment(l_seg)%direction = obc_direction_s
1345 elseif (ie_obc<is_obc)
then 1346 obc%segment(l_seg)%direction = obc_direction_n
1347 i=is_obc;is_obc=ie_obc;ie_obc=i
1350 obc%segment(l_seg)%on_pe = .false.
1353 if (len_trim(action_str(a_loop)) == 0)
then 1355 elseif (trim(action_str(a_loop)) ==
'FLATHER')
then 1356 obc%segment(l_seg)%Flather = .true.
1357 obc%segment(l_seg)%open = .true.
1358 obc%Flather_v_BCs_exist_globally = .true.
1359 obc%open_v_BCs_exist_globally = .true.
1360 obc%segment(l_seg)%z_values_needed = .true.
1361 obc%segment(l_seg)%v_values_needed = .true.
1362 elseif (trim(action_str(a_loop)) ==
'ORLANSKI')
then 1363 obc%segment(l_seg)%radiation = .true.
1364 obc%segment(l_seg)%open = .true.
1365 obc%open_v_BCs_exist_globally = .true.
1366 obc%radiation_BCs_exist_globally = .true.
1367 elseif (trim(action_str(a_loop)) ==
'ORLANSKI_TAN')
then 1368 obc%segment(l_seg)%radiation = .true.
1369 obc%segment(l_seg)%radiation_tan = .true.
1370 obc%radiation_BCs_exist_globally = .true.
1371 elseif (trim(action_str(a_loop)) ==
'ORLANSKI_GRAD')
then 1372 obc%segment(l_seg)%radiation = .true.
1373 obc%segment(l_seg)%radiation_grad = .true.
1374 elseif (trim(action_str(a_loop)) ==
'OBLIQUE')
then 1375 obc%segment(l_seg)%oblique = .true.
1376 obc%segment(l_seg)%open = .true.
1377 obc%oblique_BCs_exist_globally = .true.
1378 obc%open_v_BCs_exist_globally = .true.
1379 elseif (trim(action_str(a_loop)) ==
'OBLIQUE_TAN')
then 1380 obc%segment(l_seg)%oblique = .true.
1381 obc%segment(l_seg)%oblique_tan = .true.
1382 obc%oblique_BCs_exist_globally = .true.
1383 elseif (trim(action_str(a_loop)) ==
'OBLIQUE_GRAD')
then 1384 obc%segment(l_seg)%oblique = .true.
1385 obc%segment(l_seg)%oblique_grad = .true.
1386 elseif (trim(action_str(a_loop)) ==
'NUDGED')
then 1387 obc%segment(l_seg)%nudged = .true.
1388 obc%nudged_v_BCs_exist_globally = .true.
1389 obc%segment(l_seg)%v_values_needed = .true.
1390 elseif (trim(action_str(a_loop)) ==
'NUDGED_TAN')
then 1391 obc%segment(l_seg)%nudged_tan = .true.
1392 obc%nudged_v_BCs_exist_globally = .true.
1393 obc%segment(l_seg)%u_values_needed = .true.
1394 elseif (trim(action_str(a_loop)) ==
'NUDGED_GRAD')
then 1395 obc%segment(l_seg)%nudged_grad = .true.
1396 obc%segment(l_seg)%g_values_needed = .true.
1397 elseif (trim(action_str(a_loop)) ==
'GRADIENT')
then 1398 obc%segment(l_seg)%gradient = .true.
1399 obc%segment(l_seg)%open = .true.
1400 obc%open_v_BCs_exist_globally = .true.
1401 elseif (trim(action_str(a_loop)) ==
'SIMPLE')
then 1402 obc%segment(l_seg)%specified = .true.
1403 obc%specified_v_BCs_exist_globally = .true.
1404 obc%segment(l_seg)%v_values_needed = .true.
1405 elseif (trim(action_str(a_loop)) ==
'SIMPLE_TAN')
then 1406 obc%segment(l_seg)%specified_tan = .true.
1407 obc%segment(l_seg)%u_values_needed = .true.
1408 elseif (trim(action_str(a_loop)) ==
'SIMPLE_GRAD')
then 1409 obc%segment(l_seg)%specified_grad = .true.
1410 obc%segment(l_seg)%g_values_needed = .true.
1412 call mom_error(fatal,
"MOM_open_boundary.F90, setup_v_point_obc: "//&
1413 "String '"//trim(action_str(a_loop))//
"' not understood.")
1415 if (obc%segment(l_seg)%nudged .or. obc%segment(l_seg)%nudged_tan)
then 1416 write(segment_param_str(1:43),
"('OBC_SEGMENT_',i3.3,'_VELOCITY_NUDGING_TIMESCALES')") l_seg
1418 call get_param(pf, mdl, segment_param_str(1:43), tnudge, &
1419 "Timescales in days for nudging along a segment, "//&
1420 "for inflow, then outflow. Setting both to zero should "//&
1421 "behave like SIMPLE obcs for the baroclinic velocities.", &
1422 fail_if_missing=.true., default=0., units=
"days", scale=86400.0*us%s_to_T)
1423 obc%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1)
1424 obc%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2)
1430 if (j_obc<=g%HI%JsdB+1 .or. j_obc>=g%HI%JedB-1)
return 1431 if (ie_obc<=g%HI%IsdB .or. is_obc>=g%HI%IedB)
return 1433 obc%segment(l_seg)%on_pe = .true.
1434 obc%segment(l_seg)%is_N_or_S = .true.
1436 do i=g%HI%isd, g%HI%ied
1437 if (i>is_obc .and. i<=ie_obc)
then 1438 obc%segnum_v(i,j_obc) = l_seg
1441 obc%segment(l_seg)%Is_obc = is_obc
1442 obc%segment(l_seg)%Ie_obc = ie_obc
1443 obc%segment(l_seg)%Js_obc = j_obc
1444 obc%segment(l_seg)%Je_obc = j_obc
1445 call allocate_obc_segment_data(obc, obc%segment(l_seg))
1447 if (obc%segment(l_seg)%oblique .and. obc%segment(l_seg)%radiation) &
1448 call mom_error(fatal,
"MOM_open_boundary.F90, setup_v_point_obc: \n"//&
1449 "Orlanski and Oblique OBC options cannot be used together on one segment.")
1451 if (obc%segment(l_seg)%u_values_needed .or. obc%segment(l_seg)%v_values_needed .or. &
1452 obc%segment(l_seg)%t_values_needed .or. obc%segment(l_seg)%s_values_needed .or. &
1453 obc%segment(l_seg)%z_values_needed .or. obc%segment(l_seg)%g_values_needed) &
1454 obc%segment(l_seg)%values_needed = .true.
1455 end subroutine setup_v_point_obc
1458 subroutine parse_segment_str(ni_global, nj_global, segment_str, l, m, n, action_str, reentrant)
1459 integer,
intent(in) :: ni_global
1460 integer,
intent(in) :: nj_global
1461 character(len=*),
intent(in) :: segment_str
1462 integer,
intent(out) :: l
1463 integer,
intent(out) :: m
1464 integer,
intent(out) :: n
1465 character(len=*),
intent(out) :: action_str(:)
1466 logical,
intent(in) :: reentrant
1468 character(len=24) :: word1, word2, m_word, n_word
1473 integer,
parameter :: halo = 10
1476 word1 = extract_word(segment_str,
',',1)
1477 word2 = extract_word(segment_str,
',',2)
1478 if (word1(1:2)==
'I=')
then 1481 if (.not. (word2(1:2)==
'J='))
call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1482 "Second word of string '"//trim(segment_str)//
"' must start with 'J='.")
1483 elseif (word1(1:2)==
'J=')
then 1486 if (.not. (word2(1:2)==
'I='))
call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1487 "Second word of string '"//trim(segment_str)//
"' must start with 'I='.")
1489 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str"//&
1490 "String '"//segment_str//
"' must start with 'I=' or 'J='.")
1494 l = interpret_int_expr( word1(3:24), l_max )
1495 if (l<0 .or. l>l_max)
then 1496 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1497 "First value from string '"//trim(segment_str)//
"' is outside of the physical domain.")
1501 m_word = extract_word(word2(3:24),
':',1)
1502 m = interpret_int_expr( m_word, mn_max )
1504 if (m<-halo .or. m>mn_max+halo)
then 1505 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1506 "Beginning of range in string '"//trim(segment_str)//
"' is outside of the physical domain.")
1509 if (m<-1 .or. m>mn_max+1)
then 1510 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1511 "Beginning of range in string '"//trim(segment_str)//
"' is outside of the physical domain.")
1516 n_word = extract_word(word2(3:24),
':',2)
1517 n = interpret_int_expr( n_word, mn_max )
1519 if (n<-halo .or. n>mn_max+halo)
then 1520 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1521 "End of range in string '"//trim(segment_str)//
"' is outside of the physical domain.")
1524 if (n<-1 .or. n>mn_max+1)
then 1525 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1526 "End of range in string '"//trim(segment_str)//
"' is outside of the physical domain.")
1530 if (abs(n-m)==0)
then 1531 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1532 "Range in string '"//trim(segment_str)//
"' must span one cell.")
1536 do j = 1,
size(action_str)
1537 action_str(j) = extract_word(segment_str,
',',2+j)
1543 integer function interpret_int_expr(string, imax)
1544 character(len=*),
intent(in) :: string
1545 integer,
intent(in) :: imax
1549 slen = len_trim(string)
1550 if (slen==0)
call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str"//&
1551 "Parsed string was empty!")
1552 if (len_trim(string)==1 .and. string(1:1)==
'N')
then 1553 interpret_int_expr = imax
1554 elseif (string(1:1)==
'N')
then 1555 if (string(2:2)==
'+')
then 1556 read(string(3:slen),*,err=911) interpret_int_expr
1557 interpret_int_expr = imax + interpret_int_expr
1558 elseif (string(2:2)==
'-')
then 1559 read(string(3:slen),*,err=911) interpret_int_expr
1560 interpret_int_expr = imax - interpret_int_expr
1563 read(string(1:slen),*,err=911) interpret_int_expr
1566 911
call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str"//&
1567 "Problem reading value from string '"//trim(string)//
"'.")
1568 end function interpret_int_expr
1569 end subroutine parse_segment_str
1573 subroutine parse_segment_manifest_str(segment_str, num_fields, fields)
1574 character(len=*),
intent(in) :: segment_str
1576 integer,
intent(out) :: num_fields
1577 character(len=*),
dimension(MAX_OBC_FIELDS),
intent(out) :: fields
1581 character(len=128) :: word1, word2
1585 word1 = extract_word(segment_str,
',', num_fields+1)
1586 if (trim(word1) ==
'')
exit 1587 num_fields = num_fields + 1
1588 word2 = extract_word(word1,
'=', 1)
1589 fields(num_fields) = trim(word2)
1591 end subroutine parse_segment_manifest_str
1595 subroutine parse_segment_data_str(segment_str, idx, var, value, filename, fieldname)
1596 character(len=*),
intent(in) :: segment_str
1598 integer,
intent(in) :: idx
1599 character(len=*),
intent(in) :: var
1600 character(len=*),
intent(out) :: filename
1601 character(len=*),
intent(out) :: fieldname
1603 real,
optional,
intent(out) ::
value 1606 character(len=128) :: word1, word2, word3, method
1610 word3 = extract_word(segment_str,
',', idx)
1611 word1 = extract_word(word3,
':', 1)
1613 word2 = extract_word(word1,
'=', 1)
1614 if (trim(word2) == trim(var))
then 1615 method = trim(extract_word(word1,
'=', 2))
1616 lword = len_trim(method)
1617 if (method(lword-3:lword) ==
'file')
then 1619 word1 = extract_word(word3,
':', 2)
1620 filename = extract_word(word1,
'(', 1)
1621 fieldname = extract_word(word1,
'(', 2)
1622 lword = len_trim(fieldname)
1623 fieldname = fieldname(1:lword-1)
1625 elseif (method(lword-4:lword) ==
'value')
then 1628 word1 = extract_word(word3,
':', 2)
1629 lword = len_trim(word1)
1630 read(word1(1:lword), *, end=986, err=987)
value 1635 986
call mom_error(fatal,
'End of record while parsing segment data specification! '//trim(segment_str))
1636 987
call mom_error(fatal,
'Error while parsing segment data specification! '//trim(segment_str))
1637 end subroutine parse_segment_data_str
1642 subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature)
1645 logical,
intent(in) :: use_temperature
1648 integer :: n,m,num_fields
1649 character(len=1024) :: segstr
1650 character(len=256) :: filename
1651 character(len=20) :: segnam, suffix
1652 character(len=32) :: varnam, fieldname
1654 character(len=32),
dimension(MAX_OBC_FIELDS) :: fields
1656 character(len=256) :: mesg
1658 do n=1, obc%number_of_segments
1659 segment => obc%segment(n)
1660 write(segnam,
"('OBC_SEGMENT_',i3.3,'_DATA')") n
1661 write(suffix,
"('_segment_',i3.3)") n
1665 if (segstr ==
'') cycle
1667 call parse_segment_manifest_str(trim(segstr), num_fields, fields)
1668 if (num_fields == 0) cycle
1672 call parse_segment_data_str(trim(segstr), m, trim(fields(m)), &
1673 value, filename, fieldname)
1674 if (trim(filename) /=
'none')
then 1675 if (fields(m) ==
'TEMP')
then 1676 if (segment%is_E_or_W_2)
then 1677 obc%tracer_x_reservoirs_used(1) = .true.
1679 obc%tracer_y_reservoirs_used(1) = .true.
1682 if (fields(m) ==
'SALT')
then 1683 if (segment%is_E_or_W_2)
then 1684 obc%tracer_x_reservoirs_used(2) = .true.
1686 obc%tracer_y_reservoirs_used(2) = .true.
1692 if (use_temperature)
then 1693 if (segment%is_E_or_W_2)
then 1694 obc%tracer_x_reservoirs_used(1) = .true.
1695 obc%tracer_x_reservoirs_used(2) = .true.
1697 obc%tracer_y_reservoirs_used(1) = .true.
1698 obc%tracer_y_reservoirs_used(2) = .true.
1705 end subroutine parse_for_tracer_reservoirs
1708 subroutine parse_segment_param_real(segment_str, var, param_value, debug )
1709 character(len=*),
intent(in) :: segment_str
1711 character(len=*),
intent(in) :: var
1712 real,
intent(out) :: param_value
1713 logical,
optional,
intent(in) :: debug
1715 character(len=128) :: word1, word2, word3, method
1716 integer :: lword, nfields, n, m
1717 logical :: continue,dbg
1718 character(len=32),
dimension(MAX_OBC_FIELDS) :: flds
1723 if (
PRESENT(debug)) dbg=debug
1726 word1 = extract_word(segment_str,
',',nfields+1)
1727 if (trim(word1) ==
'')
exit 1729 word2 = extract_word(word1,
'=',1)
1730 flds(nfields) = trim(word2)
1747 if (trim(var)==trim(flds(n)))
then 1757 word3 = extract_word(segment_str,
',',m)
1760 word2 = extract_word(word1,
'=',1)
1761 if (trim(word2) == trim(var))
then 1762 method=trim(extract_word(word1,
'=',2))
1763 lword=len_trim(method)
1764 read(method(1:lword),*,err=987) param_value
1784 986
call mom_error(fatal,
'End of record while parsing segment data specification! '//trim(segment_str))
1785 987
call mom_error(fatal,
'Error while parsing segment parameter specification! '//trim(segment_str))
1787 end subroutine parse_segment_param_real
1791 subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CSp)
1800 real :: vel2_rescale
1802 integer :: i, j, k, isd, ied, jsd, jed, nz, m
1803 integer :: isdb, iedb, jsdb, jedb
1804 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
1805 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1807 if (.not.
associated(obc))
return 1809 id_clock_pass = cpu_clock_id(
'(Ocean OBC halo updates)', grain=clock_routine)
1810 if (obc%radiation_BCs_exist_globally)
call pass_vector(obc%rx_normal, obc%ry_normal, g%Domain, &
1812 if (obc%oblique_BCs_exist_globally)
call pass_vector(obc%rx_oblique, obc%ry_oblique, g%Domain, &
1814 if (
associated(obc%cff_normal))
call pass_var(obc%cff_normal, g%Domain, position=corner)
1815 if (
associated(obc%tres_x) .and.
associated(obc%tres_y))
then 1817 call pass_vector(obc%tres_x(:,:,:,m), obc%tres_y(:,:,:,m), g%Domain, to_all+scalar_pair)
1819 elseif (
associated(obc%tres_x))
then 1821 call pass_var(obc%tres_x(:,:,:,m), g%Domain, position=east_face)
1823 elseif (
associated(obc%tres_y))
then 1825 call pass_var(obc%tres_y(:,:,:,m), g%Domain, position=north_face)
1848 if ( obc%oblique_BCs_exist_globally .and. (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1849 ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) )
then 1850 vel2_rescale = (us%m_to_L * us%s_to_T_restart)**2 / (us%m_to_L_restart * us%s_to_T)**2
1852 do k=1,nz ;
do j=jsd,jed ;
do i=isdb,iedb
1853 obc%rx_oblique(i,j,k) = vel2_rescale * obc%rx_oblique(i,j,k)
1854 enddo ;
enddo ;
enddo 1857 do k=1,nz ;
do j=jsdb,jedb ;
do i=isd,ied
1858 obc%ry_oblique(i,j,k) = vel2_rescale * obc%ry_oblique(i,j,k)
1859 enddo ;
enddo ;
enddo 1862 do k=1,nz ;
do j=jsdb,jedb ;
do i=isdb,iedb
1863 obc%cff_normal(i,j,k) = vel2_rescale * obc%cff_normal(i,j,k)
1864 enddo ;
enddo ;
enddo 1868 end subroutine open_boundary_init
1870 logical function open_boundary_query(OBC, apply_open_OBC, apply_specified_OBC, apply_Flather_OBC, &
1871 apply_nudged_OBC, needs_ext_seg_data)
1873 logical,
optional,
intent(in) :: apply_open_obc
1874 logical,
optional,
intent(in) :: apply_specified_obc
1875 logical,
optional,
intent(in) :: apply_flather_obc
1876 logical,
optional,
intent(in) :: apply_nudged_obc
1877 logical,
optional,
intent(in) :: needs_ext_seg_data
1878 open_boundary_query = .false.
1879 if (.not.
associated(obc))
return 1880 if (
present(apply_open_obc)) open_boundary_query = obc%open_u_BCs_exist_globally .or. &
1881 obc%open_v_BCs_exist_globally
1882 if (
present(apply_specified_obc)) open_boundary_query = obc%specified_u_BCs_exist_globally .or. &
1883 obc%specified_v_BCs_exist_globally
1884 if (
present(apply_flather_obc)) open_boundary_query = obc%Flather_u_BCs_exist_globally .or. &
1885 obc%Flather_v_BCs_exist_globally
1886 if (
present(apply_nudged_obc)) open_boundary_query = obc%nudged_u_BCs_exist_globally .or. &
1887 obc%nudged_v_BCs_exist_globally
1888 if (
present(needs_ext_seg_data)) open_boundary_query = obc%needs_IO_for_data
1890 end function open_boundary_query
1893 subroutine open_boundary_dealloc(OBC)
1898 if (.not.
associated(obc))
return 1900 do n=1, obc%number_of_segments
1901 segment => obc%segment(n)
1902 call deallocate_obc_segment_data(obc, segment)
1904 if (
associated(obc%segment))
deallocate(obc%segment)
1905 if (
associated(obc%segnum_u))
deallocate(obc%segnum_u)
1906 if (
associated(obc%segnum_v))
deallocate(obc%segnum_v)
1907 if (
associated(obc%rx_normal))
deallocate(obc%rx_normal)
1908 if (
associated(obc%ry_normal))
deallocate(obc%ry_normal)
1909 if (
associated(obc%rx_oblique))
deallocate(obc%rx_oblique)
1910 if (
associated(obc%ry_oblique))
deallocate(obc%ry_oblique)
1911 if (
associated(obc%cff_normal))
deallocate(obc%cff_normal)
1912 if (
associated(obc%tres_x))
deallocate(obc%tres_x)
1913 if (
associated(obc%tres_y))
deallocate(obc%tres_y)
1915 end subroutine open_boundary_dealloc
1918 subroutine open_boundary_end(OBC)
1920 call open_boundary_dealloc(obc)
1921 end subroutine open_boundary_end
1924 subroutine open_boundary_impose_normal_slope(OBC, G, depth)
1927 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: depth
1932 if (.not.
associated(obc))
return 1934 if (.not.(obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
1935 obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
1938 do n=1,obc%number_of_segments
1939 segment=>obc%segment(n)
1940 if (.not. segment%on_pe) cycle
1941 if (segment%direction == obc_direction_e)
then 1943 do j=segment%HI%jsd,segment%HI%jed
1944 depth(i+1,j) = depth(i,j)
1946 elseif (segment%direction == obc_direction_w)
then 1948 do j=segment%HI%jsd,segment%HI%jed
1949 depth(i,j) = depth(i+1,j)
1951 elseif (segment%direction == obc_direction_n)
then 1953 do i=segment%HI%isd,segment%HI%ied
1954 depth(i,j+1) = depth(i,j)
1956 elseif (segment%direction == obc_direction_s)
then 1958 do i=segment%HI%isd,segment%HI%ied
1959 depth(i,j) = depth(i,j+1)
1964 end subroutine open_boundary_impose_normal_slope
1969 subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US)
1973 real,
dimension(SZIB_(G),SZJ_(G)),
intent(inout) :: areacu
1974 real,
dimension(SZI_(G),SZJB_(G)),
intent(inout) :: areacv
1978 logical :: any_u, any_v
1980 if (.not.
associated(obc))
return 1982 do n=1,obc%number_of_segments
1983 segment=>obc%segment(n)
1984 if (.not. segment%on_pe) cycle
1985 if (segment%is_E_or_W)
then 1989 do j=segment%HI%jsd,segment%HI%jed
1990 if (g%mask2dCu(i,j) == 0) obc%segnum_u(i,j) = obc_none
1991 if (segment%direction == obc_direction_w)
then 1994 g%mask2dT(i+1,j) = 0
1997 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1998 if (segment%direction == obc_direction_w)
then 2001 g%mask2dCv(i+1,j) = 0
2007 do i=segment%HI%isd,segment%HI%ied
2008 if (g%mask2dCv(i,j) == 0) obc%segnum_v(i,j) = obc_none
2009 if (segment%direction == obc_direction_s)
then 2012 g%mask2dT(i,j+1) = 0
2015 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2016 if (segment%direction == obc_direction_s)
then 2019 g%mask2dCu(i,j+1) = 0
2025 do n=1,obc%number_of_segments
2026 segment=>obc%segment(n)
2027 if (.not. segment%on_pe .or. .not. segment%specified) cycle
2028 if (segment%is_E_or_W)
then 2031 do j=segment%HI%jsd,segment%HI%jed
2032 if (segment%direction == obc_direction_e)
then 2033 areacu(i,j) = g%areaT(i,j)
2035 areacu(i,j) = g%areaT(i+1,j)
2041 do i=segment%HI%isd,segment%HI%ied
2042 if (segment%direction == obc_direction_s)
then 2043 areacv(i,j) = g%areaT(i,j+1)
2045 areacu(i,j) = g%areaT(i,j)
2057 do n=1,obc%number_of_segments
2058 segment=>obc%segment(n)
2059 if (.not. segment%on_pe) cycle
2060 if (segment%is_E_or_W)
then 2062 do j=segment%HI%jsd,segment%HI%jed
2063 if (obc%segnum_u(i,j) /= obc_none) any_u = .true.
2067 do i=segment%HI%isd,segment%HI%ied
2068 if (obc%segnum_v(i,j) /= obc_none) any_v = .true.
2074 if (.not.(any_u .or. any_v)) obc%OBC_pe = .false.
2076 end subroutine open_boundary_impose_land_mask
2079 subroutine setup_obc_tracer_reservoirs(G, OBC)
2084 integer :: i, j, k, m, n
2086 do n=1,obc%number_of_segments
2087 segment=>obc%segment(n)
2088 if (
associated(segment%tr_Reg))
then 2089 if (segment%is_E_or_W)
then 2092 if (
associated(segment%tr_Reg%Tr(m)%tres))
then 2094 do j=segment%HI%jsd,segment%HI%jed
2095 obc%tres_x(i,j,k,m) = segment%tr_Reg%Tr(m)%t(i,j,k)
2103 if (
associated(segment%tr_Reg%Tr(m)%tres))
then 2105 do i=segment%HI%isd,segment%HI%ied
2106 obc%tres_y(i,j,k,m) = segment%tr_Reg%Tr(m)%t(i,j,k)
2115 end subroutine setup_obc_tracer_reservoirs
2118 subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, US, dt)
2121 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u_new
2124 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_old
2125 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v_new
2128 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_old
2130 real,
intent(in) :: dt
2132 real :: dhdt, dhdx, dhdy
2133 real :: gamma_u, gamma_2
2135 real :: rx_max, ry_max
2136 real :: rx_new, rx_avg
2137 real :: ry_new, ry_avg
2138 real :: cff_new, cff_avg
2139 real,
allocatable,
dimension(:,:,:) :: &
2149 integer :: i, j, k, is, ie, js, je, m, nz, n
2150 integer :: is_obc, ie_obc, js_obc, je_obc
2152 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2154 if (.not.
associated(obc))
return 2156 if (.not.(obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
2159 eps = 1.0e-20*us%m_s_to_L_T**2
2164 if (obc%gamma_uv < 1.0)
then 2165 do n=1,obc%number_of_segments
2166 segment=>obc%segment(n)
2167 if (.not. segment%on_pe) cycle
2168 if (segment%is_E_or_W .and. segment%radiation)
then 2171 do j=segment%HI%jsd,segment%HI%jed
2172 segment%rx_norm_rad(i,j,k) = obc%rx_normal(i,j,k)
2175 elseif (segment%is_N_or_S .and. segment%radiation)
then 2178 do i=segment%HI%isd,segment%HI%ied
2179 segment%ry_norm_rad(i,j,k) = obc%ry_normal(i,j,k)
2183 if (segment%is_E_or_W .and. segment%oblique)
then 2186 do j=segment%HI%jsd,segment%HI%jed
2187 segment%rx_norm_obl(i,j,k) = obc%rx_oblique(i,j,k)
2188 segment%ry_norm_obl(i,j,k) = obc%ry_oblique(i,j,k)
2189 segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
2192 elseif (segment%is_N_or_S .and. segment%oblique)
then 2195 do i=segment%HI%isd,segment%HI%ied
2196 segment%rx_norm_obl(i,j,k) = obc%rx_oblique(i,j,k)
2197 segment%ry_norm_obl(i,j,k) = obc%ry_oblique(i,j,k)
2198 segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
2206 do n=1,obc%number_of_segments
2207 segment=>obc%segment(n)
2208 if (
associated(segment%tr_Reg))
then 2209 if (segment%is_E_or_W)
then 2212 if (
associated(segment%tr_Reg%Tr(m)%tres))
then 2214 do j=segment%HI%jsd,segment%HI%jed
2215 segment%tr_Reg%Tr(m)%tres(i,j,k) = obc%tres_x(i,j,k,m)
2223 if (
associated(segment%tr_Reg%Tr(m)%tres))
then 2225 do i=segment%HI%isd,segment%HI%ied
2226 segment%tr_Reg%Tr(m)%tres(i,j,k) = obc%tres_y(i,j,k,m)
2235 gamma_u = obc%gamma_uv
2236 rx_max = obc%rx_max ; ry_max = obc%rx_max
2237 do n=1,obc%number_of_segments
2238 segment=>obc%segment(n)
2239 if (.not. segment%on_pe) cycle
2240 if (segment%oblique)
call gradient_at_q_points(g, segment, u_new(:,:,:), v_new(:,:,:))
2241 if (segment%direction == obc_direction_e)
then 2243 if (i<g%HI%IscB) cycle
2244 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
2245 if (segment%radiation)
then 2246 dhdt = (u_old(i-1,j,k) - u_new(i-1,j,k))
2247 dhdx = (u_new(i-1,j,k) - u_new(i-2,j,k))
2249 if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
2250 if (gamma_u < 1.0)
then 2251 rx_avg = (1.0-gamma_u)*segment%rx_norm_rad(i,j,k) + gamma_u*rx_new
2255 segment%rx_norm_rad(i,j,k) = rx_avg
2259 segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) / (1.0+rx_avg)
2262 if (gamma_u < 1.0)
then 2263 obc%rx_normal(i,j,k) = segment%rx_norm_rad(i,j,k)
2265 elseif (segment%oblique)
then 2266 dhdt = (u_old(i-1,j,k) - u_new(i-1,j,k))
2267 dhdx = (u_new(i-1,j,k) - u_new(i-2,j,k))
2268 if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0)
then 2269 dhdy = segment%grad_normal(j-1,1,k)
2270 elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0)
then 2273 dhdy = segment%grad_normal(j,1,k)
2275 if (dhdt*dhdx < 0.0) dhdt = 0.0
2276 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2277 rx_new = min(dhdt*dhdx, cff_new*rx_max)
2278 ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2279 if (gamma_u < 1.0)
then 2280 rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2281 ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2282 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2288 segment%rx_norm_obl(i,j,k) = rx_avg
2289 segment%ry_norm_obl(i,j,k) = ry_avg
2290 segment%cff_normal(i,j,k) = cff_avg
2291 segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) - &
2292 (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + &
2293 min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
2295 if (gamma_u < 1.0)
then 2298 obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2299 obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2300 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2302 elseif (segment%gradient)
then 2303 segment%normal_vel(i,j,k) = u_new(i-1,j,k)
2305 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then 2307 if (dhdt*dhdx <= 0.0)
then 2308 tau = segment%Velocity_nudging_timescale_in
2310 tau = segment%Velocity_nudging_timescale_out
2312 gamma_2 = dt / (tau + dt)
2313 segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2314 gamma_2 * segment%nudged_normal_vel(i,j,k)
2317 if (segment%radiation_tan .or. segment%radiation_grad)
then 2319 allocate(rx_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2321 if (gamma_u < 1.0)
then 2322 rx_tang_rad(i,segment%HI%JsdB,k) = segment%rx_norm_rad(i,segment%HI%jsd,k)
2323 rx_tang_rad(i,segment%HI%JedB,k) = segment%rx_norm_rad(i,segment%HI%jed,k)
2324 do j=segment%HI%JsdB+1,segment%HI%JedB-1
2325 rx_tang_rad(i,j,k) = 0.5*(segment%rx_norm_rad(i,j,k) + segment%rx_norm_rad(i,j+1,k))
2328 do j=segment%HI%JsdB,segment%HI%JedB
2329 dhdt = v_old(i,j,k)-v_new(i,j,k)
2330 dhdx = v_new(i,j,k)-v_new(i-1,j,k)
2331 rx_tang_rad(i,j,k) = 0.0
2332 if (dhdt*dhdx > 0.0) rx_tang_rad(i,j,k) = min( (dhdt/dhdx), rx_max)
2336 if (segment%radiation_tan)
then 2337 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2338 rx_avg = rx_tang_rad(i,j,k)
2339 segment%tangential_vel(i,j,k) = (v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) / (1.0+rx_avg)
2342 if (segment%nudged_tan)
then 2343 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2345 if (rx_tang_rad(i,j,k) <= 0.0)
then 2346 tau = segment%Velocity_nudging_timescale_in
2348 tau = segment%Velocity_nudging_timescale_out
2350 gamma_2 = dt / (tau + dt)
2351 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2352 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2355 if (segment%radiation_grad)
then 2356 js_obc = max(segment%HI%JsdB,g%jsd+1)
2357 je_obc = min(segment%HI%JedB,g%jed-1)
2358 do k=1,nz ;
do j=js_obc,je_obc
2359 rx_avg = rx_tang_rad(i,j,k)
2369 segment%tangential_grad(i,j,k) = ((v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
2370 rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) / (1.0+rx_avg)
2373 if (segment%nudged_grad)
then 2374 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2376 if (rx_tang_rad(i,j,k) <= 0.0)
then 2377 tau = segment%Velocity_nudging_timescale_in
2379 tau = segment%Velocity_nudging_timescale_out
2381 gamma_2 = dt / (tau + dt)
2382 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2383 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2386 deallocate(rx_tang_rad)
2388 if (segment%oblique_tan .or. segment%oblique_grad)
then 2390 allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2391 allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2392 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2394 if (gamma_u < 1.0)
then 2395 rx_tang_obl(i,segment%HI%JsdB,k) = segment%rx_norm_obl(i,segment%HI%jsd,k)
2396 rx_tang_obl(i,segment%HI%JedB,k) = segment%rx_norm_obl(i,segment%HI%jed,k)
2397 ry_tang_obl(i,segment%HI%JsdB,k) = segment%ry_norm_obl(i,segment%HI%jsd,k)
2398 ry_tang_obl(i,segment%HI%JedB,k) = segment%ry_norm_obl(i,segment%HI%jed,k)
2399 cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
2400 cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
2401 do j=segment%HI%JsdB+1,segment%HI%JedB-1
2402 rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i,j+1,k))
2403 ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i,j+1,k))
2404 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
2407 do j=segment%HI%JsdB,segment%HI%JedB
2408 dhdt = v_old(i,j,k)-v_new(i,j,k)
2409 dhdx = v_new(i,j,k)-v_new(i-1,j,k)
2410 if (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) > 0.0)
then 2411 dhdy = segment%grad_tan(j,1,k)
2412 elseif (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) == 0.0)
then 2415 dhdy = segment%grad_tan(j+1,1,k)
2417 if (dhdt*dhdx < 0.0) dhdt = 0.0
2418 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2419 rx_new = min(dhdt*dhdx, cff_new*rx_max)
2420 ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2421 rx_tang_obl(i,j,k) = rx_new
2422 ry_tang_obl(i,j,k) = ry_new
2423 cff_tangential(i,j,k) = cff_new
2427 if (segment%oblique_tan)
then 2428 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2429 rx_avg = rx_tang_obl(i,j,k)
2430 ry_avg = ry_tang_obl(i,j,k)
2431 cff_avg = cff_tangential(i,j,k)
2432 segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) - &
2433 (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + &
2434 min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
2438 if (segment%nudged_tan)
then 2439 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2441 if (rx_tang_obl(i,j,k) <= 0.0)
then 2442 tau = segment%Velocity_nudging_timescale_in
2444 tau = segment%Velocity_nudging_timescale_out
2446 gamma_2 = dt / (tau + dt)
2447 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2448 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2451 if (segment%oblique_grad)
then 2452 js_obc = max(segment%HI%JsdB,g%jsd+1)
2453 je_obc = min(segment%HI%JedB,g%jed-1)
2454 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
2455 rx_avg = rx_tang_obl(i,j,k)
2456 ry_avg = ry_tang_obl(i,j,k)
2457 cff_avg = cff_tangential(i,j,k)
2458 segment%tangential_grad(i,j,k) = &
2459 ((cff_avg*(v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
2460 rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) - &
2461 (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + &
2462 min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k)) ) / &
2466 if (segment%nudged_grad)
then 2467 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2469 if (rx_tang_obl(i,j,k) <= 0.0)
then 2470 tau = segment%Velocity_nudging_timescale_in
2472 tau = segment%Velocity_nudging_timescale_out
2474 gamma_2 = dt / (tau + dt)
2475 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2476 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2479 deallocate(rx_tang_obl)
2480 deallocate(ry_tang_obl)
2481 deallocate(cff_tangential)
2485 if (segment%direction == obc_direction_w)
then 2487 if (i>g%HI%IecB) cycle
2488 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
2489 if (segment%radiation)
then 2490 dhdt = (u_old(i+1,j,k) - u_new(i+1,j,k))
2491 dhdx = (u_new(i+1,j,k) - u_new(i+2,j,k))
2493 if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
2494 if (gamma_u < 1.0)
then 2495 rx_avg = (1.0-gamma_u)*segment%rx_norm_rad(i,j,k) + gamma_u*rx_new
2499 segment%rx_norm_rad(i,j,k) = rx_avg
2503 segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) / (1.0+rx_avg)
2504 if (gamma_u < 1.0)
then 2507 obc%rx_normal(i,j,k) = segment%rx_norm_rad(i,j,k)
2509 elseif (segment%oblique)
then 2510 dhdt = (u_old(i+1,j,k) - u_new(i+1,j,k))
2511 dhdx = (u_new(i+1,j,k) - u_new(i+2,j,k))
2512 if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0)
then 2513 dhdy = segment%grad_normal(j-1,1,k)
2514 elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0)
then 2517 dhdy = segment%grad_normal(j,1,k)
2519 if (dhdt*dhdx < 0.0) dhdt = 0.0
2521 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2522 rx_new = min(dhdt*dhdx, cff_new*rx_max)
2523 ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2524 if (gamma_u < 1.0)
then 2525 rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2526 ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2527 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2533 segment%rx_norm_obl(i,j,k) = rx_avg
2534 segment%ry_norm_obl(i,j,k) = ry_avg
2535 segment%cff_normal(i,j,k) = cff_avg
2536 segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) - &
2537 (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + &
2538 min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
2540 if (gamma_u < 1.0)
then 2543 obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2544 obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2545 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2547 elseif (segment%gradient)
then 2548 segment%normal_vel(i,j,k) = u_new(i+1,j,k)
2550 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then 2552 if (dhdt*dhdx <= 0.0)
then 2553 tau = segment%Velocity_nudging_timescale_in
2555 tau = segment%Velocity_nudging_timescale_out
2557 gamma_2 = dt / (tau + dt)
2558 segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2559 gamma_2 * segment%nudged_normal_vel(i,j,k)
2562 if (segment%radiation_tan .or. segment%radiation_grad)
then 2564 allocate(rx_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2566 if (gamma_u < 1.0)
then 2567 rx_tang_rad(i,segment%HI%JsdB,k) = segment%rx_norm_rad(i,segment%HI%jsd,k)
2568 rx_tang_rad(i,segment%HI%JedB,k) = segment%rx_norm_rad(i,segment%HI%jed,k)
2569 do j=segment%HI%JsdB+1,segment%HI%JedB-1
2570 rx_tang_rad(i,j,k) = 0.5*(segment%rx_norm_rad(i,j,k) + segment%rx_norm_rad(i,j+1,k))
2573 do j=segment%HI%JsdB,segment%HI%JedB
2574 dhdt = v_old(i+1,j,k)-v_new(i+1,j,k)
2575 dhdx = v_new(i+1,j,k)-v_new(i+2,j,k)
2576 rx_tang_rad(i,j,k) = 0.0
2577 if (dhdt*dhdx > 0.0) rx_tang_rad(i,j,k) = min( (dhdt/dhdx), rx_max)
2581 if (segment%radiation_tan)
then 2582 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2583 rx_avg = rx_tang_rad(i,j,k)
2584 segment%tangential_vel(i,j,k) = (v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) / (1.0+rx_avg)
2587 if (segment%nudged_tan)
then 2588 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2590 if (rx_tang_rad(i,j,k) <= 0.0)
then 2591 tau = segment%Velocity_nudging_timescale_in
2593 tau = segment%Velocity_nudging_timescale_out
2595 gamma_2 = dt / (tau + dt)
2596 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2597 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2600 if (segment%radiation_grad)
then 2601 js_obc = max(segment%HI%JsdB,g%jsd+1)
2602 je_obc = min(segment%HI%JedB,g%jed-1)
2603 do k=1,nz ;
do j=js_obc,je_obc
2604 rx_avg = rx_tang_rad(i,j,k)
2614 segment%tangential_grad(i,j,k) = ((v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
2615 rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) / (1.0+rx_avg)
2618 if (segment%nudged_grad)
then 2619 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2621 if (rx_tang_rad(i,j,k) <= 0.0)
then 2622 tau = segment%Velocity_nudging_timescale_in
2624 tau = segment%Velocity_nudging_timescale_out
2626 gamma_2 = dt / (tau + dt)
2627 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2628 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2631 deallocate(rx_tang_rad)
2633 if (segment%oblique_tan .or. segment%oblique_grad)
then 2635 allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2636 allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2637 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2639 if (gamma_u < 1.0)
then 2640 rx_tang_obl(i,segment%HI%JsdB,k) = segment%rx_norm_obl(i,segment%HI%jsd,k)
2641 rx_tang_obl(i,segment%HI%JedB,k) = segment%rx_norm_obl(i,segment%HI%jed,k)
2642 ry_tang_obl(i,segment%HI%JsdB,k) = segment%ry_norm_obl(i,segment%HI%jsd,k)
2643 ry_tang_obl(i,segment%HI%JedB,k) = segment%ry_norm_obl(i,segment%HI%jed,k)
2644 cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
2645 cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
2646 do j=segment%HI%JsdB+1,segment%HI%JedB-1
2647 rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i,j+1,k))
2648 ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i,j+1,k))
2649 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
2652 do j=segment%HI%JsdB,segment%HI%JedB
2653 dhdt = v_old(i+1,j,k)-v_new(i+1,j,k)
2654 dhdx = v_new(i+1,j,k)-v_new(i+2,j,k)
2655 if (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) > 0.0)
then 2656 dhdy = segment%grad_tan(j,1,k)
2657 elseif (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) == 0.0)
then 2660 dhdy = segment%grad_tan(j+1,1,k)
2662 if (dhdt*dhdx < 0.0) dhdt = 0.0
2663 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2664 rx_new = min(dhdt*dhdx, cff_new*rx_max)
2665 ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2666 rx_tang_obl(i,j,k) = rx_new
2667 ry_tang_obl(i,j,k) = ry_new
2668 cff_tangential(i,j,k) = cff_new
2672 if (segment%oblique_tan)
then 2673 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2674 rx_avg = rx_tang_obl(i,j,k)
2675 ry_avg = ry_tang_obl(i,j,k)
2676 cff_avg = cff_tangential(i,j,k)
2677 segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) - &
2678 (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + &
2679 min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
2683 if (segment%nudged_tan)
then 2684 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2686 if (rx_tang_obl(i,j,k) <= 0.0)
then 2687 tau = segment%Velocity_nudging_timescale_in
2689 tau = segment%Velocity_nudging_timescale_out
2691 gamma_2 = dt / (tau + dt)
2692 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2693 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2696 if (segment%oblique_grad)
then 2697 js_obc = max(segment%HI%JsdB,g%jsd+1)
2698 je_obc = min(segment%HI%JedB,g%jed-1)
2699 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
2700 rx_avg = rx_tang_obl(i,j,k)
2701 ry_avg = ry_tang_obl(i,j,k)
2702 cff_avg = cff_tangential(i,j,k)
2703 segment%tangential_grad(i,j,k) = &
2704 ((cff_avg*(v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
2705 rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) - &
2706 (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + &
2707 min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k))) / &
2711 if (segment%nudged_grad)
then 2712 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2714 if (rx_tang_obl(i,j,k) <= 0.0)
then 2715 tau = segment%Velocity_nudging_timescale_in
2717 tau = segment%Velocity_nudging_timescale_out
2719 gamma_2 = dt / (tau + dt)
2720 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2721 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2724 deallocate(rx_tang_obl)
2725 deallocate(ry_tang_obl)
2726 deallocate(cff_tangential)
2730 if (segment%direction == obc_direction_n)
then 2732 if (j<g%HI%JscB) cycle
2733 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
2734 if (segment%radiation)
then 2735 dhdt = (v_old(i,j-1,k) - v_new(i,j-1,k))
2736 dhdy = (v_new(i,j-1,k) - v_new(i,j-2,k))
2738 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2739 if (gamma_u < 1.0)
then 2740 ry_avg = (1.0-gamma_u)*segment%ry_norm_rad(i,j,k) + gamma_u*ry_new
2744 segment%ry_norm_rad(i,j,k) = ry_avg
2748 segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) / (1.0+ry_avg)
2749 if (gamma_u < 1.0)
then 2752 obc%ry_normal(i,j,k) = segment%ry_norm_rad(i,j,k)
2754 elseif (segment%oblique)
then 2755 dhdt = (v_old(i,j-1,k) - v_new(i,j-1,k))
2756 dhdy = (v_new(i,j-1,k) - v_new(i,j-2,k))
2757 if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0)
then 2758 dhdx = segment%grad_normal(i-1,1,k)
2759 elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0)
then 2762 dhdx = segment%grad_normal(i,1,k)
2764 if (dhdt*dhdy < 0.0) dhdt = 0.0
2765 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2766 ry_new = min(dhdt*dhdy, cff_new*ry_max)
2767 rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
2768 if (gamma_u < 1.0)
then 2769 rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2770 ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2771 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2777 segment%rx_norm_obl(i,j,k) = rx_avg
2778 segment%ry_norm_obl(i,j,k) = ry_avg
2779 segment%cff_normal(i,j,k) = cff_avg
2780 segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) - &
2781 (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) +&
2782 min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
2784 if (gamma_u < 1.0)
then 2787 obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2788 obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2789 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2791 elseif (segment%gradient)
then 2792 segment%normal_vel(i,j,k) = v_new(i,j-1,k)
2794 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then 2796 if (dhdt*dhdy <= 0.0)
then 2797 tau = segment%Velocity_nudging_timescale_in
2799 tau = segment%Velocity_nudging_timescale_out
2801 gamma_2 = dt / (tau + dt)
2802 segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2803 gamma_2 * segment%nudged_normal_vel(i,j,k)
2806 if (segment%radiation_tan .or. segment%radiation_grad)
then 2808 allocate(ry_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2810 if (gamma_u < 1.0)
then 2811 ry_tang_rad(segment%HI%IsdB,j,k) = segment%ry_norm_rad(segment%HI%isd,j,k)
2812 ry_tang_rad(segment%HI%IedB,j,k) = segment%ry_norm_rad(segment%HI%ied,j,k)
2813 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2814 ry_tang_rad(i,j,k) = 0.5*(segment%ry_norm_rad(i,j,k) + segment%ry_norm_rad(i+1,j,k))
2817 do i=segment%HI%IsdB,segment%HI%IedB
2818 dhdt = u_old(i,j-1,k)-u_new(i,j-1,k)
2819 dhdy = u_new(i,j-1,k)-u_new(i,j-2,k)
2820 ry_tang_rad(i,j,k) = 0.0
2821 if (dhdt*dhdy > 0.0) ry_tang_rad(i,j,k) = min( (dhdt/dhdy), rx_max)
2825 if (segment%radiation_tan)
then 2826 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2827 ry_avg = ry_tang_rad(i,j,k)
2828 segment%tangential_vel(i,j,k) = (u_new(i,j,k) + ry_avg*u_new(i,j-1,k)) / (1.0+ry_avg)
2831 if (segment%nudged_tan)
then 2832 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2834 if (ry_tang_rad(i,j,k) <= 0.0)
then 2835 tau = segment%Velocity_nudging_timescale_in
2837 tau = segment%Velocity_nudging_timescale_out
2839 gamma_2 = dt / (tau + dt)
2840 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2841 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2844 if (segment%radiation_grad)
then 2845 is_obc = max(segment%HI%IsdB,g%isd+1)
2846 ie_obc = min(segment%HI%IedB,g%ied-1)
2847 do k=1,nz ;
do i=is_obc,ie_obc
2848 ry_avg = ry_tang_rad(i,j,k)
2858 segment%tangential_grad(i,j,k) = ((u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
2859 ry_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) / (1.0+ry_avg)
2862 if (segment%nudged_grad)
then 2863 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2865 if (ry_tang_rad(i,j,k) <= 0.0)
then 2866 tau = segment%Velocity_nudging_timescale_in
2868 tau = segment%Velocity_nudging_timescale_out
2870 gamma_2 = dt / (tau + dt)
2871 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2872 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2875 deallocate(ry_tang_rad)
2877 if (segment%oblique_tan .or. segment%oblique_grad)
then 2879 allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2880 allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2881 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2883 if (gamma_u < 1.0)
then 2884 rx_tang_obl(segment%HI%IsdB,j,k) = segment%rx_norm_obl(segment%HI%isd,j,k)
2885 rx_tang_obl(segment%HI%IedB,j,k) = segment%rx_norm_obl(segment%HI%ied,j,k)
2886 ry_tang_obl(segment%HI%IsdB,j,k) = segment%ry_norm_obl(segment%HI%isd,j,k)
2887 ry_tang_obl(segment%HI%IedB,j,k) = segment%ry_norm_obl(segment%HI%ied,j,k)
2888 cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
2889 cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
2890 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2891 rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i+1,j,k))
2892 ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i+1,j,k))
2893 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
2896 do i=segment%HI%IsdB,segment%HI%IedB
2897 dhdt = u_old(i,j,k)-u_new(i,j,k)
2898 dhdy = u_new(i,j,k)-u_new(i,j-1,k)
2899 if (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) > 0.0)
then 2900 dhdx = segment%grad_tan(i,1,k)
2901 elseif (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) == 0.0)
then 2904 dhdx = segment%grad_tan(i+1,1,k)
2906 if (dhdt*dhdy < 0.0) dhdt = 0.0
2907 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2908 ry_new = min(dhdt*dhdy, cff_new*ry_max)
2909 rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
2910 rx_tang_obl(i,j,k) = rx_new
2911 ry_tang_obl(i,j,k) = ry_new
2912 cff_tangential(i,j,k) = cff_new
2916 if (segment%oblique_tan)
then 2917 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2918 rx_avg = rx_tang_obl(i,j,k)
2919 ry_avg = ry_tang_obl(i,j,k)
2920 cff_avg = cff_tangential(i,j,k)
2921 segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + ry_avg*u_new(i,j-1,k)) - &
2922 (max(rx_avg,0.0)*segment%grad_tan(i,2,k) + &
2923 min(rx_avg,0.0)*segment%grad_tan(i+1,2,k))) / &
2927 if (segment%nudged_tan)
then 2928 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2930 if (ry_tang_obl(i,j,k) <= 0.0)
then 2931 tau = segment%Velocity_nudging_timescale_in
2933 tau = segment%Velocity_nudging_timescale_out
2935 gamma_2 = dt / (tau + dt)
2936 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2937 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2940 if (segment%oblique_grad)
then 2941 is_obc = max(segment%HI%IsdB,g%isd+1)
2942 ie_obc = min(segment%HI%IedB,g%ied-1)
2943 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2944 rx_avg = rx_tang_obl(i,j,k)
2945 ry_avg = ry_tang_obl(i,j,k)
2946 cff_avg = cff_tangential(i,j,k)
2947 segment%tangential_grad(i,j,k) = &
2948 ((cff_avg*(u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
2949 ry_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) - &
2950 (max(rx_avg,0.0)*segment%grad_gradient(i,2,k) + &
2951 min(rx_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
2955 if (segment%nudged_grad)
then 2956 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2958 if (ry_tang_obl(i,j,k) <= 0.0)
then 2959 tau = segment%Velocity_nudging_timescale_in
2961 tau = segment%Velocity_nudging_timescale_out
2963 gamma_2 = dt / (tau + dt)
2964 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2965 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2968 deallocate(rx_tang_obl)
2969 deallocate(ry_tang_obl)
2970 deallocate(cff_tangential)
2974 if (segment%direction == obc_direction_s)
then 2976 if (j>g%HI%JecB) cycle
2977 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
2978 if (segment%radiation)
then 2979 dhdt = (v_old(i,j+1,k) - v_new(i,j+1,k))
2980 dhdy = (v_new(i,j+1,k) - v_new(i,j+2,k))
2982 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2983 if (gamma_u < 1.0)
then 2984 ry_avg = (1.0-gamma_u)*segment%ry_norm_rad(i,j,k) + gamma_u*ry_new
2988 segment%ry_norm_rad(i,j,k) = ry_avg
2992 segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) / (1.0+ry_avg)
2993 if (gamma_u < 1.0)
then 2996 obc%ry_normal(i,j,k) = segment%ry_norm_rad(i,j,k)
2998 elseif (segment%oblique)
then 2999 dhdt = (v_old(i,j+1,k) - v_new(i,j+1,k))
3000 dhdy = (v_new(i,j+1,k) - v_new(i,j+2,k))
3001 if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0)
then 3002 dhdx = segment%grad_normal(i-1,1,k)
3003 elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0)
then 3006 dhdx = segment%grad_normal(i,1,k)
3008 if (dhdt*dhdy < 0.0) dhdt = 0.0
3010 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
3011 ry_new = min(dhdt*dhdy, cff_new*ry_max)
3012 rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
3013 if (gamma_u < 1.0)
then 3014 rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
3015 ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
3016 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
3022 segment%rx_norm_obl(i,j,k) = rx_avg
3023 segment%ry_norm_obl(i,j,k) = ry_avg
3024 segment%cff_normal(i,j,k) = cff_avg
3025 segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) - &
3026 (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) + &
3027 min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
3029 if (gamma_u < 1.0)
then 3032 obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
3033 obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
3034 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
3036 elseif (segment%gradient)
then 3037 segment%normal_vel(i,j,k) = v_new(i,j+1,k)
3039 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then 3041 if (dhdt*dhdy <= 0.0)
then 3042 tau = segment%Velocity_nudging_timescale_in
3044 tau = segment%Velocity_nudging_timescale_out
3046 gamma_2 = dt / (tau + dt)
3047 segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
3048 gamma_2 * segment%nudged_normal_vel(i,j,k)
3051 if (segment%radiation_tan .or. segment%radiation_grad)
then 3053 allocate(ry_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
3055 if (gamma_u < 1.0)
then 3056 ry_tang_rad(segment%HI%IsdB,j,k) = segment%ry_norm_rad(segment%HI%isd,j,k)
3057 ry_tang_rad(segment%HI%IedB,j,k) = segment%ry_norm_rad(segment%HI%ied,j,k)
3058 do i=segment%HI%IsdB+1,segment%HI%IedB-1
3059 ry_tang_rad(i,j,k) = 0.5*(segment%ry_norm_rad(i,j,k) + segment%ry_norm_rad(i+1,j,k))
3062 do i=segment%HI%IsdB,segment%HI%IedB
3063 dhdt = u_old(i,j+1,k)-u_new(i,j+1,k)
3064 dhdy = u_new(i,j+1,k)-u_new(i,j+2,k)
3065 ry_tang_rad(i,j,k) = 0.0
3066 if (dhdt*dhdy > 0.0) ry_tang_rad(i,j,k) = min( (dhdt/dhdy), rx_max)
3070 if (segment%radiation_tan)
then 3071 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3072 ry_avg = ry_tang_rad(i,j,k)
3073 segment%tangential_vel(i,j,k) = (u_new(i,j+1,k) + ry_avg*u_new(i,j+2,k)) / (1.0+ry_avg)
3076 if (segment%nudged_tan)
then 3077 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3079 if (ry_tang_rad(i,j,k) <= 0.0)
then 3080 tau = segment%Velocity_nudging_timescale_in
3082 tau = segment%Velocity_nudging_timescale_out
3084 gamma_2 = dt / (tau + dt)
3085 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
3086 gamma_2 * segment%nudged_tangential_vel(i,j,k)
3089 if (segment%radiation_grad)
then 3090 is_obc = max(segment%HI%IsdB,g%isd+1)
3091 ie_obc = min(segment%HI%IedB,g%ied-1)
3092 do k=1,nz ;
do i=is_obc,ie_obc
3093 ry_avg = ry_tang_rad(i,j,k)
3103 segment%tangential_grad(i,j,k) = ((u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
3104 ry_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) / (1.0+ry_avg)
3107 if (segment%nudged_grad)
then 3108 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3110 if (ry_tang_rad(i,j,k) <= 0.0)
then 3111 tau = segment%Velocity_nudging_timescale_in
3113 tau = segment%Velocity_nudging_timescale_out
3115 gamma_2 = dt / (tau + dt)
3116 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
3117 gamma_2 * segment%nudged_tangential_grad(i,j,k)
3120 deallocate(ry_tang_rad)
3122 if (segment%oblique_tan .or. segment%oblique_grad)
then 3124 allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
3125 allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
3126 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
3128 if (gamma_u < 1.0)
then 3129 rx_tang_obl(segment%HI%IsdB,j,k) = segment%rx_norm_obl(segment%HI%isd,j,k)
3130 rx_tang_obl(segment%HI%IedB,j,k) = segment%rx_norm_obl(segment%HI%ied,j,k)
3131 ry_tang_obl(segment%HI%IsdB,j,k) = segment%ry_norm_obl(segment%HI%isd,j,k)
3132 ry_tang_obl(segment%HI%IedB,j,k) = segment%ry_norm_obl(segment%HI%ied,j,k)
3133 cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
3134 cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
3135 do i=segment%HI%IsdB+1,segment%HI%IedB-1
3136 rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i+1,j,k))
3137 ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i+1,j,k))
3138 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
3141 do i=segment%HI%IsdB,segment%HI%IedB
3142 dhdt = u_old(i,j+1,k)-u_new(i,j+1,k)
3143 dhdy = u_new(i,j+1,k)-u_new(i,j+2,k)
3144 if (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) > 0.0)
then 3145 dhdx = segment%grad_tan(i,1,k)
3146 elseif (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) == 0.0)
then 3149 dhdx = segment%grad_tan(i+1,1,k)
3151 if (dhdt*dhdy < 0.0) dhdt = 0.0
3152 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
3153 ry_new = min(dhdt*dhdy, cff_new*ry_max)
3154 rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
3155 rx_tang_obl(i,j,k) = rx_new
3156 ry_tang_obl(i,j,k) = ry_new
3157 cff_tangential(i,j,k) = cff_new
3161 if (segment%oblique_tan)
then 3162 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3163 rx_avg = rx_tang_obl(i,j,k)
3164 ry_avg = ry_tang_obl(i,j,k)
3165 cff_avg = cff_tangential(i,j,k)
3166 segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j+1,k) + ry_avg*u_new(i,j+2,k)) - &
3167 (max(rx_avg,0.0)*segment%grad_tan(i,2,k) + &
3168 min(rx_avg,0.0)*segment%grad_tan(i+1,2,k)) ) / &
3172 if (segment%nudged_tan)
then 3173 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3175 if (ry_tang_obl(i,j,k) <= 0.0)
then 3176 tau = segment%Velocity_nudging_timescale_in
3178 tau = segment%Velocity_nudging_timescale_out
3180 gamma_2 = dt / (tau + dt)
3181 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
3182 gamma_2 * segment%nudged_tangential_vel(i,j,k)
3185 if (segment%oblique_grad)
then 3186 is_obc = max(segment%HI%IsdB,g%isd+1)
3187 ie_obc = min(segment%HI%IedB,g%ied-1)
3188 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
3189 rx_avg = rx_tang_obl(i,j,k)
3190 ry_avg = ry_tang_obl(i,j,k)
3191 cff_avg = cff_tangential(i,j,k)
3192 segment%tangential_grad(i,j,k) = &
3193 ((cff_avg*(u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
3194 ry_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) - &
3195 (max(rx_avg,0.0)*segment%grad_gradient(i,2,k) + &
3196 min(rx_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
3200 if (segment%nudged_grad)
then 3201 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
3203 if (ry_tang_obl(i,j,k) <= 0.0)
then 3204 tau = segment%Velocity_nudging_timescale_in
3206 tau = segment%Velocity_nudging_timescale_out
3208 gamma_2 = dt / (tau + dt)
3209 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
3210 gamma_2 * segment%nudged_tangential_grad(i,j,k)
3213 deallocate(rx_tang_obl)
3214 deallocate(ry_tang_obl)
3215 deallocate(cff_tangential)
3221 call open_boundary_apply_normal_flow(obc, g, u_new, v_new)
3223 call pass_vector(u_new, v_new, g%Domain, clock=id_clock_pass)
3225 end subroutine radiation_open_bdry_conds
3228 subroutine open_boundary_apply_normal_flow(OBC, G, u, v)
3232 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
3234 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
3237 integer :: i, j, k, n
3240 if (.not.
associated(obc))
return 3242 do n=1,obc%number_of_segments
3243 segment => obc%segment(n)
3244 if (.not. segment%on_pe)
then 3246 elseif (segment%radiation .or. segment%oblique .or. segment%gradient)
then 3247 if (segment%is_E_or_W)
then 3249 do k=1,g%ke ;
do j=segment%HI%jsd,segment%HI%jed
3250 u(i,j,k) = segment%normal_vel(i,j,k)
3252 elseif (segment%is_N_or_S)
then 3254 do k=1,g%ke ;
do i=segment%HI%isd,segment%HI%ied
3255 v(i,j,k) = segment%normal_vel(i,j,k)
3261 end subroutine open_boundary_apply_normal_flow
3264 subroutine open_boundary_zero_normal_flow(OBC, G, u, v)
3268 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
3269 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
3271 integer :: i, j, k, n
3274 if (.not.
associated(obc))
return 3276 do n=1,obc%number_of_segments
3277 segment => obc%segment(n)
3278 if (.not. segment%on_pe)
then 3280 elseif (segment%is_E_or_W)
then 3282 do k=1,g%ke ;
do j=segment%HI%jsd,segment%HI%jed
3285 elseif (segment%is_N_or_S)
then 3287 do k=1,g%ke ;
do i=segment%HI%isd,segment%HI%ied
3293 end subroutine open_boundary_zero_normal_flow
3296 subroutine gradient_at_q_points(G, segment, uvel, vvel)
3299 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: uvel
3300 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: vvel
3303 if (.not. segment%on_pe)
return 3305 if (segment%is_E_or_W)
then 3306 if (segment%direction == obc_direction_e)
then 3309 do j=max(segment%HI%JsdB, g%HI%JsdB+1),min(segment%HI%JedB, g%HI%JedB-1)
3310 segment%grad_normal(j,1,k) = (uvel(i-1,j+1,k)-uvel(i-1,j,k)) * g%mask2dBu(i-1,j)
3311 segment%grad_normal(j,2,k) = (uvel(i,j+1,k)-uvel(i,j,k)) * g%mask2dBu(i,j)
3314 if (segment%oblique_tan)
then 3316 do j=max(segment%HI%jsd-1, g%HI%jsd),min(segment%HI%jed+1, g%HI%jed)
3317 segment%grad_tan(j,1,k) = (vvel(i-1,j,k)-vvel(i-1,j-1,k)) * g%mask2dT(i-1,j)
3318 segment%grad_tan(j,2,k) = (vvel(i,j,k)-vvel(i,j-1,k)) * g%mask2dT(i,j)
3322 if (segment%oblique_grad)
then 3324 do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
3325 segment%grad_gradient(j,1,k) = (((vvel(i-1,j,k) - vvel(i-2,j,k))*g%IdxBu(i-2,j)) - &
3326 (vvel(i-1,j-1,k) - vvel(i-2,j-1,k))*g%IdxBu(i-2,j-1)) * g%mask2dCu(i-2,j)
3327 segment%grad_gradient(j,2,k) = (((vvel(i,j,k) - vvel(i-1,j,k))*g%IdxBu(i-1,j)) - &
3328 (vvel(i,j-1,k) - vvel(i-1,j-1,k))*g%IdxBu(i-1,j-1)) * g%mask2dCu(i-1,j)
3335 do j=max(segment%HI%JsdB, g%HI%JsdB+1),min(segment%HI%JedB, g%HI%JedB-1)
3336 segment%grad_normal(j,1,k) = (uvel(i+1,j+1,k)-uvel(i+1,j,k)) * g%mask2dBu(i+1,j)
3337 segment%grad_normal(j,2,k) = (uvel(i,j+1,k)-uvel(i,j,k)) * g%mask2dBu(i,j)
3340 if (segment%oblique_tan)
then 3342 do j=max(segment%HI%jsd-1, g%HI%jsd),min(segment%HI%jed+1, g%HI%jed)
3343 segment%grad_tan(j,1,k) = (vvel(i+2,j,k)-vvel(i+2,j-1,k)) * g%mask2dT(i+2,j)
3344 segment%grad_tan(j,2,k) = (vvel(i+1,j,k)-vvel(i+1,j-1,k)) * g%mask2dT(i+1,j)
3348 if (segment%oblique_grad)
then 3350 do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
3351 segment%grad_gradient(j,1,k) = (((vvel(i+3,j,k) - vvel(i+2,j,k))*g%IdxBu(i+2,j)) - &
3352 (vvel(i+3,j-1,k) - vvel(i+2,j-1,k))*g%IdxBu(i+2,j-1)) * g%mask2dCu(i+2,j)
3353 segment%grad_gradient(j,2,k) = (((vvel(i+2,j,k) - vvel(i+1,j,k))*g%IdxBu(i+1,j)) - &
3354 (vvel(i+2,j-1,k) - vvel(i+1,j-1,k))*g%IdxBu(i+1,j-1)) * g%mask2dCu(i+1,j)
3359 elseif (segment%is_N_or_S)
then 3360 if (segment%direction == obc_direction_n)
then 3363 do i=max(segment%HI%IsdB, g%HI%IsdB+1),min(segment%HI%IedB, g%HI%IedB-1)
3364 segment%grad_normal(i,1,k) = (vvel(i+1,j-1,k)-vvel(i,j-1,k)) * g%mask2dBu(i,j-1)
3365 segment%grad_normal(i,2,k) = (vvel(i+1,j,k)-vvel(i,j,k)) * g%mask2dBu(i,j)
3368 if (segment%oblique_tan)
then 3370 do i=max(segment%HI%isd-1, g%HI%isd),min(segment%HI%ied+1, g%HI%ied)
3371 segment%grad_tan(i,1,k) = (uvel(i,j-1,k)-uvel(i-1,j-1,k)) * g%mask2dT(i,j-1)
3372 segment%grad_tan(i,2,k) = (uvel(i,j,k)-uvel(i-1,j,k)) * g%mask2dT(i,j)
3376 if (segment%oblique_grad)
then 3378 do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
3379 segment%grad_gradient(i,1,k) = (((uvel(i,j-1,k) - uvel(i,j-2,k))*g%IdyBu(i,j-2)) - &
3380 (uvel(i-1,j-1,k) - uvel(i-1,j-2,k))*g%IdyBu(i-1,j-2)) * g%mask2dCv(i,j-2)
3381 segment%grad_gradient(i,2,k) = (((uvel(i,j,k) - uvel(i,j-1,k))*g%IdyBu(i,j-1)) - &
3382 (uvel(i-1,j,k) - uvel(i-1,j-1,k))*g%IdyBu(i-1,j-1)) * g%mask2dCv(i,j-1)
3389 do i=max(segment%HI%IsdB, g%HI%IsdB+1),min(segment%HI%IedB, g%HI%IedB-1)
3390 segment%grad_normal(i,1,k) = (vvel(i+1,j+1,k)-vvel(i,j+1,k)) * g%mask2dBu(i,j+1)
3391 segment%grad_normal(i,2,k) = (vvel(i+1,j,k)-vvel(i,j,k)) * g%mask2dBu(i,j)
3394 if (segment%oblique_tan)
then 3396 do i=max(segment%HI%isd-1, g%HI%isd),min(segment%HI%ied+1, g%HI%ied)
3397 segment%grad_tan(i,1,k) = (uvel(i,j+2,k)-uvel(i-1,j+2,k)) * g%mask2dT(i,j+2)
3398 segment%grad_tan(i,2,k) = (uvel(i,j+1,k)-uvel(i-1,j+1,k)) * g%mask2dT(i,j+1)
3402 if (segment%oblique_grad)
then 3404 do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
3405 segment%grad_gradient(i,1,k) = (((uvel(i,j+3,k) - uvel(i,j+2,k))*g%IdyBu(i,j+2)) - &
3406 (uvel(i-1,j+3,k) - uvel(i-1,j+2,k))*g%IdyBu(i-1,j+2)) * g%mask2dCv(i,j+2)
3407 segment%grad_gradient(i,2,k) = (((uvel(i,j+2,k) - uvel(i,j+1,k))*g%IdyBu(i,j+1)) - &
3408 (uvel(i-1,j+2,k) - uvel(i-1,j+1,k))*g%IdyBu(i-1,j+1)) * g%mask2dCv(i,j+1)
3415 end subroutine gradient_at_q_points
3420 subroutine set_tracer_data(OBC, tv, h, G, PF, tracer_Reg)
3424 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(inout) :: h
3428 integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, n
3429 integer :: isd_off, jsd_off
3430 integer :: isdb, iedb, jsdb, jedb
3432 character(len=40) :: mdl =
"set_tracer_data" 3433 character(len=200) :: filename, obc_file, inputdir
3435 real :: temp_u(g%domain%niglobal+1,g%domain%njglobal)
3436 real :: temp_v(g%domain%niglobal,g%domain%njglobal+1)
3438 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3439 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3440 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3446 if (
associated(tv%T))
then 3451 do n=1,obc%number_of_segments
3452 segment => obc%segment(n)
3453 if (.not. segment%on_pe) cycle
3455 if (segment%direction == obc_direction_e)
then 3457 do k=1,g%ke ;
do j=segment%HI%jsd,segment%HI%jed
3458 tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k)
3460 elseif (segment%direction == obc_direction_w)
then 3462 do k=1,g%ke ;
do j=segment%HI%jsd,segment%HI%jed
3463 tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k)
3465 elseif (segment%direction == obc_direction_n)
then 3467 do k=1,g%ke ;
do i=segment%HI%isd,segment%HI%ied
3468 tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k)
3470 elseif (segment%direction == obc_direction_s)
then 3472 do k=1,g%ke ;
do i=segment%HI%isd,segment%HI%ied
3473 tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k)
3479 end subroutine set_tracer_data
3482 function lookup_seg_field(OBC_seg,field)
3484 character(len=32),
intent(in) :: field
3485 integer :: lookup_seg_field
3490 do n=1,obc_seg%num_fields
3491 if (trim(field) == obc_seg%field(n)%name)
then 3497 end function lookup_seg_field
3501 subroutine allocate_obc_segment_data(OBC, segment)
3505 integer :: isd, ied, jsd, jed
3506 integer :: IsdB, IedB, JsdB, JedB
3507 integer :: IscB, IecB, JscB, JecB
3508 character(len=40) :: mdl =
"allocate_OBC_segment_data" 3510 isd = segment%HI%isd ; ied = segment%HI%ied
3511 jsd = segment%HI%jsd ; jed = segment%HI%jed
3512 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
3513 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
3514 iscb = segment%HI%IscB ; iecb = segment%HI%IecB
3515 jscb = segment%HI%JscB ; jecb = segment%HI%JecB
3517 if (.not. segment%on_pe)
return 3519 if (segment%is_E_or_W)
then 3521 allocate(segment%Cg(isdb:iedb,jsd:jed)); segment%Cg(:,:)=0.
3522 allocate(segment%Htot(isdb:iedb,jsd:jed)); segment%Htot(:,:)=0.0
3523 allocate(segment%h(isdb:iedb,jsd:jed,obc%ke)); segment%h(:,:,:)=0.0
3524 allocate(segment%eta(isdb:iedb,jsd:jed)); segment%eta(:,:)=0.0
3525 if (segment%radiation)
then 3526 allocate(segment%rx_norm_rad(isdb:iedb,jsd:jed,obc%ke)); segment%rx_norm_rad(:,:,:)=0.0
3528 allocate(segment%normal_vel(isdb:iedb,jsd:jed,obc%ke)); segment%normal_vel(:,:,:)=0.0
3529 allocate(segment%normal_vel_bt(isdb:iedb,jsd:jed)); segment%normal_vel_bt(:,:)=0.0
3530 allocate(segment%normal_trans(isdb:iedb,jsd:jed,obc%ke)); segment%normal_trans(:,:,:)=0.0
3531 if (segment%nudged)
then 3532 allocate(segment%nudged_normal_vel(isdb:iedb,jsd:jed,obc%ke)); segment%nudged_normal_vel(:,:,:)=0.0
3534 if (segment%radiation_tan .or. segment%nudged_tan .or. segment%specified_tan .or. &
3535 segment%oblique_tan .or. obc%computed_vorticity .or. obc%computed_strain)
then 3536 allocate(segment%tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_vel(:,:,:)=0.0
3538 if (segment%nudged_tan)
then 3539 allocate(segment%nudged_tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_vel(:,:,:)=0.0
3541 if (segment%nudged_grad)
then 3542 allocate(segment%nudged_tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
3544 if (obc%specified_vorticity .or. obc%specified_strain .or. segment%radiation_grad .or. &
3545 segment%oblique_grad .or. segment%specified_grad)
then 3546 allocate(segment%tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_grad(:,:,:)=0.0
3548 if (segment%oblique)
then 3549 allocate(segment%grad_normal(jsdb:jedb,2,obc%ke)); segment%grad_normal(:,:,:) = 0.0
3550 allocate(segment%rx_norm_obl(isdb:iedb,jsd:jed,obc%ke)); segment%rx_norm_obl(:,:,:)=0.0
3551 allocate(segment%ry_norm_obl(isdb:iedb,jsd:jed,obc%ke)); segment%ry_norm_obl(:,:,:)=0.0
3552 allocate(segment%cff_normal(isdb:iedb,jsd:jed,obc%ke)); segment%cff_normal(:,:,:)=0.0
3554 if (segment%oblique_tan)
then 3555 allocate(segment%grad_tan(jsd-1:jed+1,2,obc%ke)); segment%grad_tan(:,:,:) = 0.0
3557 if (segment%oblique_grad)
then 3558 allocate(segment%grad_gradient(jsd:jed,2,obc%ke)); segment%grad_gradient(:,:,:) = 0.0
3562 if (segment%is_N_or_S)
then 3564 allocate(segment%Cg(isd:ied,jsdb:jedb)); segment%Cg(:,:)=0.
3565 allocate(segment%Htot(isd:ied,jsdb:jedb)); segment%Htot(:,:)=0.0
3566 allocate(segment%h(isd:ied,jsdb:jedb,obc%ke)); segment%h(:,:,:)=0.0
3567 allocate(segment%eta(isd:ied,jsdb:jedb)); segment%eta(:,:)=0.0
3568 if (segment%radiation)
then 3569 allocate(segment%ry_norm_rad(isd:ied,jsdb:jedb,obc%ke)); segment%ry_norm_rad(:,:,:)=0.0
3571 allocate(segment%normal_vel(isd:ied,jsdb:jedb,obc%ke)); segment%normal_vel(:,:,:)=0.0
3572 allocate(segment%normal_vel_bt(isd:ied,jsdb:jedb)); segment%normal_vel_bt(:,:)=0.0
3573 allocate(segment%normal_trans(isd:ied,jsdb:jedb,obc%ke)); segment%normal_trans(:,:,:)=0.0
3574 if (segment%nudged)
then 3575 allocate(segment%nudged_normal_vel(isd:ied,jsdb:jedb,obc%ke)); segment%nudged_normal_vel(:,:,:)=0.0
3577 if (segment%radiation_tan .or. segment%nudged_tan .or. segment%specified_tan .or. &
3578 segment%oblique_tan .or. obc%computed_vorticity .or. obc%computed_strain)
then 3579 allocate(segment%tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_vel(:,:,:)=0.0
3581 if (segment%nudged_tan)
then 3582 allocate(segment%nudged_tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_vel(:,:,:)=0.0
3584 if (segment%nudged_grad)
then 3585 allocate(segment%nudged_tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
3587 if (obc%specified_vorticity .or. obc%specified_strain .or. segment%radiation_grad .or. &
3588 segment%oblique_grad .or. segment%specified_grad)
then 3589 allocate(segment%tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_grad(:,:,:)=0.0
3591 if (segment%oblique)
then 3592 allocate(segment%grad_normal(isdb:iedb,2,obc%ke)); segment%grad_normal(:,:,:) = 0.0
3593 allocate(segment%rx_norm_obl(isd:ied,jsdb:jedb,obc%ke)); segment%rx_norm_obl(:,:,:)=0.0
3594 allocate(segment%ry_norm_obl(isd:ied,jsdb:jedb,obc%ke)); segment%ry_norm_obl(:,:,:)=0.0
3595 allocate(segment%cff_normal(isd:ied,jsdb:jedb,obc%ke)); segment%cff_normal(:,:,:)=0.0
3597 if (segment%oblique_tan)
then 3598 allocate(segment%grad_tan(isd-1:ied+1,2,obc%ke)); segment%grad_tan(:,:,:) = 0.0
3600 if (segment%oblique_grad)
then 3601 allocate(segment%grad_gradient(isd:ied,2,obc%ke)); segment%grad_gradient(:,:,:) = 0.0
3605 end subroutine allocate_obc_segment_data
3608 subroutine deallocate_obc_segment_data(OBC, segment)
3612 character(len=40) :: mdl =
"deallocate_OBC_segment_data" 3614 if (.not. segment%on_pe)
return 3616 if (
associated (segment%Cg))
deallocate(segment%Cg)
3617 if (
associated (segment%Htot))
deallocate(segment%Htot)
3618 if (
associated (segment%h))
deallocate(segment%h)
3619 if (
associated (segment%eta))
deallocate(segment%eta)
3620 if (
associated (segment%rx_norm_rad))
deallocate(segment%rx_norm_rad)
3621 if (
associated (segment%ry_norm_rad))
deallocate(segment%ry_norm_rad)
3622 if (
associated (segment%rx_norm_obl))
deallocate(segment%rx_norm_obl)
3623 if (
associated (segment%ry_norm_obl))
deallocate(segment%ry_norm_obl)
3624 if (
associated (segment%cff_normal))
deallocate(segment%cff_normal)
3625 if (
associated (segment%grad_normal))
deallocate(segment%grad_normal)
3626 if (
associated (segment%grad_tan))
deallocate(segment%grad_tan)
3627 if (
associated (segment%grad_gradient))
deallocate(segment%grad_gradient)
3628 if (
associated (segment%normal_vel))
deallocate(segment%normal_vel)
3629 if (
associated (segment%normal_vel_bt))
deallocate(segment%normal_vel_bt)
3630 if (
associated (segment%normal_trans))
deallocate(segment%normal_trans)
3631 if (
associated (segment%nudged_normal_vel))
deallocate(segment%nudged_normal_vel)
3632 if (
associated (segment%tangential_vel))
deallocate(segment%tangential_vel)
3633 if (
associated (segment%nudged_tangential_vel))
deallocate(segment%nudged_tangential_vel)
3634 if (
associated (segment%nudged_tangential_grad))
deallocate(segment%nudged_tangential_grad)
3635 if (
associated (segment%tangential_grad))
deallocate(segment%tangential_grad)
3636 if (
associated (segment%tr_Reg))
call segment_tracer_registry_end(segment%tr_Reg)
3639 end subroutine deallocate_obc_segment_data
3644 subroutine open_boundary_test_extern_uv(G, OBC, u, v)
3647 real,
dimension(SZIB_(G),SZJ_(G), SZK_(G)),
intent(inout) :: u
3648 real,
dimension(SZI_(G),SZJB_(G), SZK_(G)),
intent(inout) :: v
3650 integer :: i, j, k, n
3652 if (.not.
associated(obc))
return 3654 do n = 1, obc%number_of_segments
3656 if (obc%segment(n)%is_N_or_S)
then 3657 j = obc%segment(n)%HI%JsdB
3658 if (obc%segment(n)%direction == obc_direction_n)
then 3659 do i = obc%segment(n)%HI%IsdB, obc%segment(n)%HI%IedB
3660 u(i,j+1,k) = obc%silly_u
3663 do i = obc%segment(n)%HI%IsdB, obc%segment(n)%HI%IedB
3664 u(i,j,k) = obc%silly_u
3667 elseif (obc%segment(n)%is_E_or_W)
then 3668 i = obc%segment(n)%HI%IsdB
3669 if (obc%segment(n)%direction == obc_direction_e)
then 3670 do j = obc%segment(n)%HI%JsdB, obc%segment(n)%HI%JedB
3671 v(i+1,j,k) = obc%silly_u
3674 do j = obc%segment(n)%HI%JsdB, obc%segment(n)%HI%JedB
3675 v(i,j,k) = obc%silly_u
3682 end subroutine open_boundary_test_extern_uv
3687 subroutine open_boundary_test_extern_h(G, GV, OBC, h)
3691 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(inout) :: h
3694 integer :: i, j, k, n
3696 if (.not.
associated(obc))
return 3698 silly_h = gv%Z_to_H*obc%silly_h
3700 do n = 1, obc%number_of_segments
3702 if (obc%segment(n)%is_N_or_S)
then 3703 j = obc%segment(n)%HI%JsdB
3704 if (obc%segment(n)%direction == obc_direction_n)
then 3705 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
3706 h(i,j+1,k) = silly_h
3709 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
3713 elseif (obc%segment(n)%is_E_or_W)
then 3714 i = obc%segment(n)%HI%IsdB
3715 if (obc%segment(n)%direction == obc_direction_e)
then 3716 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
3717 h(i+1,j,k) = silly_h
3720 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
3728 end subroutine open_boundary_test_extern_h
3731 subroutine update_obc_segment_data(G, GV, US, OBC, tv, h, Time)
3737 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(inout) :: h
3738 type(time_type),
intent(in) :: time
3740 integer :: c, i, j, k, is, ie, js, je, isd, ied, jsd, jed
3741 integer :: isdb, iedb, jsdb, jedb, n, m, nz
3742 character(len=40) :: mdl =
"update_OBC_segment_data" 3743 character(len=200) :: filename, obc_file, inputdir
3745 integer,
dimension(4) :: siz
3746 real,
dimension(:,:,:),
pointer :: tmp_buffer_in => null()
3747 integer :: ni_seg, nj_seg
3748 integer :: ni_buf, nj_buf
3750 integer :: is_obc, ie_obc, js_obc, je_obc
3751 integer :: ishift, jshift
3752 real,
dimension(:,:),
pointer :: seg_vel => null()
3753 real,
dimension(:,:),
pointer :: seg_trans => null()
3754 real,
dimension(:,:,:),
allocatable,
target :: tmp_buffer
3755 real,
dimension(:),
allocatable :: h_stack
3756 integer :: is_obc2, js_obc2
3757 real :: net_h_src, net_h_int, scl_fac
3758 real :: tidal_vel, tidal_elev
3759 real,
pointer,
dimension(:,:) :: normal_trans_bt=>null()
3763 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3764 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3765 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3770 if (.not.
associated(obc))
return 3772 if (obc%add_tide_constituents) time_delta = time_type_to_real(time - obc%time_ref)
3774 do n = 1, obc%number_of_segments
3775 segment => obc%segment(n)
3777 if (.not. segment%on_pe) cycle
3780 ni_seg = segment%ie_obc-segment%is_obc+1
3781 nj_seg = segment%je_obc-segment%js_obc+1
3782 is_obc = max(segment%is_obc,isd-1)
3783 ie_obc = min(segment%ie_obc,ied)
3784 js_obc = max(segment%js_obc,jsd-1)
3785 je_obc = min(segment%je_obc,jed)
3798 if (segment%is_E_or_W)
then 3799 allocate(normal_trans_bt(segment%HI%IsdB:segment%HI%IedB,segment%HI%jsd:segment%HI%jed))
3800 normal_trans_bt(:,:)=0.0
3801 if (segment%direction == obc_direction_w) ishift=1
3803 do j=segment%HI%jsd,segment%HI%jed
3804 segment%Cg(i,j) = sqrt(gv%g_prime(1)*g%bathyT(i+ishift,j))
3805 segment%Htot(i,j)=0.0
3807 segment%h(i,j,k) = h(i+ishift,j,k)
3808 segment%Htot(i,j)=segment%Htot(i,j)+segment%h(i,j,k)
3812 allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB))
3813 normal_trans_bt(:,:)=0.0
3814 if (segment%direction == obc_direction_s) jshift=1
3816 do i=segment%HI%isd,segment%HI%ied
3817 segment%Cg(i,j) = sqrt(gv%g_prime(1)*g%bathyT(i,j+jshift))
3818 segment%Htot(i,j)=0.0
3820 segment%h(i,j,k) = h(i,j+jshift,k)
3821 segment%Htot(i,j)=segment%Htot(i,j)+segment%h(i,j,k)
3826 allocate(h_stack(g%ke))
3828 do m = 1,segment%num_fields
3829 if (segment%field(m)%fid > 0)
then 3830 siz(1)=
size(segment%field(m)%buffer_src,1)
3831 siz(2)=
size(segment%field(m)%buffer_src,2)
3832 siz(3)=
size(segment%field(m)%buffer_src,3)
3833 if (.not.
associated(segment%field(m)%buffer_dst))
then 3834 if (siz(3) /= segment%field(m)%nk_src)
call mom_error(fatal,
'nk_src inconsistency')
3835 if (segment%field(m)%nk_src > 1)
then 3836 if (segment%is_E_or_W)
then 3837 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then 3838 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3839 elseif (segment%field(m)%name ==
'Vamp' .or. segment%field(m)%name ==
'Vphase')
then 3840 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,siz(3)))
3841 elseif (segment%field(m)%name ==
'Uamp' .or. segment%field(m)%name ==
'Uphase' .or. &
3842 segment%field(m)%name ==
'SSHamp' .or. segment%field(m)%name ==
'SSHphase')
then 3843 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,siz(3)))
3845 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
3848 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then 3849 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3850 elseif (segment%field(m)%name ==
'Uamp' .or. segment%field(m)%name ==
'Uphase')
then 3851 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,siz(3)))
3852 elseif (segment%field(m)%name ==
'Vamp' .or. segment%field(m)%name ==
'Vphase' .or. &
3853 segment%field(m)%name ==
'SSHamp' .or. segment%field(m)%name ==
'SSHphase')
then 3854 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,siz(3)))
3856 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
3860 if (segment%is_E_or_W)
then 3861 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX' .or. &
3862 segment%field(m)%name ==
'Vamp' .or. segment%field(m)%name ==
'Vphase')
then 3863 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
3865 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,1))
3868 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY' .or. &
3869 segment%field(m)%name ==
'Uamp' .or. segment%field(m)%name ==
'Uphase')
then 3870 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
3872 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,1))
3876 segment%field(m)%buffer_dst(:,:,:)=0.0
3881 if (obc%brushcutter_mode)
then 3882 allocate(tmp_buffer(1,nj_seg*2-1,segment%field(m)%nk_src))
3884 allocate(tmp_buffer(1,nj_seg,segment%field(m)%nk_src))
3887 if (obc%brushcutter_mode)
then 3888 allocate(tmp_buffer(ni_seg*2-1,1,segment%field(m)%nk_src))
3890 allocate(tmp_buffer(ni_seg,1,segment%field(m)%nk_src))
3898 if (turns /= 0)
then 3899 if (modulo(turns, 2) /= 0)
then 3900 allocate(tmp_buffer_in(
size(tmp_buffer, 2),
size(tmp_buffer, 1),
size(tmp_buffer, 3)))
3902 allocate(tmp_buffer_in(
size(tmp_buffer, 1),
size(tmp_buffer, 2),
size(tmp_buffer, 3)))
3905 tmp_buffer_in => tmp_buffer
3908 call time_interp_external(segment%field(m)%fid,time, tmp_buffer_in)
3910 if (turns /= 0)
then 3912 if (segment%is_E_or_W &
3913 .and. .not. (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'Vamp' &
3914 .or. segment%field(m)%name ==
'Vphase' .or. segment%field(m)%name ==
'DVDX'))
then 3915 nj_buf =
size(tmp_buffer, 2) - 1
3916 call rotate_array(tmp_buffer_in(:nj_buf,:,:), turns, tmp_buffer(:,:nj_buf,:))
3917 elseif (segment%is_N_or_S &
3918 .and. .not. (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'Uamp' &
3919 .or. segment%field(m)%name ==
'Uphase' .or. segment%field(m)%name ==
'DUDY'))
then 3920 ni_buf =
size(tmp_buffer, 1) - 1
3921 call rotate_array(tmp_buffer_in(:,:ni_buf,:), turns, tmp_buffer(:ni_buf,:,:))
3927 if (segment%field(m)%name ==
'U' &
3928 .or. segment%field(m)%name ==
'DVDX' &
3929 .or. segment%field(m)%name ==
'DUDY' &
3930 .or. segment%field(m)%name ==
'Uamp')
then 3931 tmp_buffer(:,:,:) = -tmp_buffer(:,:,:)
3935 if (obc%brushcutter_mode)
then 3936 if (segment%is_E_or_W)
then 3937 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX' .or. &
3938 segment%field(m)%name ==
'Vamp' .or. segment%field(m)%name ==
'Vphase')
then 3939 segment%field(m)%buffer_src(is_obc,:,:) = &
3940 tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset)+1:2,:)
3942 segment%field(m)%buffer_src(is_obc,:,:) = &
3943 tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset):2,:)
3946 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY' .or. &
3947 segment%field(m)%name ==
'Uamp' .or. segment%field(m)%name ==
'Uphase')
then 3948 segment%field(m)%buffer_src(:,js_obc,:) = &
3949 tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset)+1:2,1,:)
3951 segment%field(m)%buffer_src(:,js_obc,:) = &
3952 tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset):2,1,:)
3956 if (segment%is_E_or_W)
then 3957 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX' .or. &
3958 segment%field(m)%name ==
'Vamp' .or. segment%field(m)%name ==
'Vphase')
then 3959 segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset+1,:)
3961 segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset,:)
3964 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY' .or. &
3965 segment%field(m)%name ==
'Uamp' .or. segment%field(m)%name ==
'Uphase')
then 3966 segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset+1,1,:)
3968 segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset,1,:)
3973 if (segment%field(m)%nk_src > 1 .and.&
3974 (index(segment%field(m)%name,
'phase') .le. 0 .and. index(segment%field(m)%name,
'amp') .le. 0))
then 3975 call time_interp_external(segment%field(m)%fid_dz,time, tmp_buffer_in)
3976 if (turns /= 0)
then 3978 if (segment%is_E_or_W &
3979 .and. .not. (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX'))
then 3980 nj_buf =
size(tmp_buffer, 2) - 1
3981 call rotate_array(tmp_buffer_in(:nj_buf,:,:), turns, tmp_buffer(:,:nj_buf,:))
3982 elseif (segment%is_N_or_S &
3983 .and. .not. (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY'))
then 3984 ni_buf =
size(tmp_buffer, 1) - 1
3985 call rotate_array(tmp_buffer_in(:,:ni_buf,:), turns, tmp_buffer(:ni_buf,:,:))
3990 if (obc%brushcutter_mode)
then 3991 if (segment%is_E_or_W)
then 3992 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then 3993 segment%field(m)%dz_src(is_obc,:,:) = &
3994 tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset)+1:2,:)
3996 segment%field(m)%dz_src(is_obc,:,:) = &
3997 tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset):2,:)
4000 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then 4001 segment%field(m)%dz_src(:,js_obc,:) = &
4002 tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset)+1:2,1,:)
4004 segment%field(m)%dz_src(:,js_obc,:) = &
4005 tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset):2,1,:)
4009 if (segment%is_E_or_W)
then 4010 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then 4011 segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset+1,:)
4013 segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset,:)
4016 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then 4017 segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset+1,1,:)
4019 segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset,1,:)
4024 call adjustsegmentetatofitbathymetry(g,gv,us,segment,m)
4026 if (segment%is_E_or_W)
then 4028 if (segment%direction == obc_direction_e) ishift=0
4030 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then 4032 do j=max(js_obc,jsd),min(je_obc,jed-1)
4035 segment%field(m)%buffer_dst(i,j,:)=0.0
4036 if (g%mask2dCu(i,j)>0. .and. g%mask2dCu(i,j+1)>0.)
then 4037 h_stack(:) = 0.5*(h(i+ishift,j,:) + h(i+ishift,j+1,:))
4038 call remapping_core_h(obc%remap_CS, &
4039 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
4040 segment%field(m)%buffer_src(i,j,:), &
4041 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
4042 elseif (g%mask2dCu(i,j)>0.)
then 4043 h_stack(:) = h(i+ishift,j,:)
4044 call remapping_core_h(obc%remap_CS, &
4045 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
4046 segment%field(m)%buffer_src(i,j,:), &
4047 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
4048 elseif (g%mask2dCu(i,j+1)>0.)
then 4049 h_stack(:) = h(i+ishift,j+1,:)
4050 call remapping_core_h(obc%remap_CS, &
4051 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
4052 segment%field(m)%buffer_src(i,j,:), &
4053 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
4057 do j=js_obc+1,je_obc
4060 segment%field(m)%buffer_dst(i,j,:)=0.0
4061 if (g%mask2dCu(i,j)>0.)
then 4062 net_h_src = sum( segment%field(m)%dz_src(i,j,:) )
4063 net_h_int = sum( h(i+ishift,j,:) )
4064 scl_fac = net_h_int / net_h_src
4065 call remapping_core_h(obc%remap_CS, &
4066 segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,j,:), &
4067 segment%field(m)%buffer_src(i,j,:), &
4068 g%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(i,j,:))
4074 if (segment%direction == obc_direction_n) jshift=0
4076 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then 4078 do i=max(is_obc,isd),min(ie_obc,ied-1)
4079 segment%field(m)%buffer_dst(i,j,:)=0.0
4080 if (g%mask2dCv(i,j)>0. .and. g%mask2dCv(i+1,j)>0.)
then 4083 h_stack(:) = 0.5*(h(i,j+jshift,:) + h(i+1,j+jshift,:))
4084 call remapping_core_h(obc%remap_CS, &
4085 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
4086 segment%field(m)%buffer_src(i,j,:), &
4087 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
4088 elseif (g%mask2dCv(i,j)>0.)
then 4089 h_stack(:) = h(i,j+jshift,:)
4090 call remapping_core_h(obc%remap_CS, &
4091 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
4092 segment%field(m)%buffer_src(i,j,:), &
4093 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
4094 elseif (g%mask2dCv(i+1,j)>0.)
then 4095 h_stack(:) = h(i+1,j+jshift,:)
4096 call remapping_core_h(obc%remap_CS, &
4097 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
4098 segment%field(m)%buffer_src(i,j,:), &
4099 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
4103 do i=is_obc+1,ie_obc
4106 segment%field(m)%buffer_dst(i,j,:)=0.0
4107 if (g%mask2dCv(i,j)>0.)
then 4108 net_h_src = sum( segment%field(m)%dz_src(i,j,:) )
4109 net_h_int = sum( h(i,j+jshift,:) )
4110 scl_fac = net_h_int / net_h_src
4111 call remapping_core_h(obc%remap_CS, &
4112 segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,j,:), &
4113 segment%field(m)%buffer_src(i,j,:), &
4114 g%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,j,:))
4119 elseif (segment%field(m)%nk_src > 1 .and. &
4120 (index(segment%field(m)%name,
'phase') > 0 .or. index(segment%field(m)%name,
'amp') > 0))
then 4122 segment%field(m)%buffer_dst(:,:,:) = segment%field(m)%buffer_src(:,:,:)
4124 segment%field(m)%buffer_dst(:,:,1) = segment%field(m)%buffer_src(:,:,1)
4126 deallocate(tmp_buffer)
4128 deallocate(tmp_buffer_in)
4130 if (.not.
associated(segment%field(m)%buffer_dst))
then 4131 if (segment%is_E_or_W)
then 4132 if (segment%field(m)%name ==
'V')
then 4133 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
4134 else if (segment%field(m)%name ==
'Vamp' .or. segment%field(m)%name ==
'Vphase')
then 4135 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
4136 elseif (segment%field(m)%name ==
'U')
then 4137 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
4138 elseif (segment%field(m)%name ==
'Uamp' .or. segment%field(m)%name ==
'Uphase')
then 4139 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,1))
4140 elseif (segment%field(m)%name ==
'DVDX')
then 4141 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
4142 elseif (segment%field(m)%name ==
'SSH' .or. segment%field(m)%name ==
'SSHamp' &
4143 .or. segment%field(m)%name ==
'SSHphase')
then 4144 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
4146 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
4149 if (segment%field(m)%name ==
'U')
then 4150 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
4151 elseif (segment%field(m)%name ==
'Uamp' .or. segment%field(m)%name ==
'Uphase')
then 4152 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
4153 elseif (segment%field(m)%name ==
'V')
then 4154 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
4155 elseif (segment%field(m)%name ==
'Vamp' .or. segment%field(m)%name ==
'Vphase')
then 4156 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,1))
4157 elseif (segment%field(m)%name ==
'DUDY')
then 4158 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
4159 elseif (segment%field(m)%name ==
'SSH' .or. segment%field(m)%name ==
'SSHamp' &
4160 .or. segment%field(m)%name ==
'SSHphase')
then 4161 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
4163 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
4166 segment%field(m)%buffer_dst(:,:,:) = segment%field(m)%value
4172 do m = 1,segment%num_fields
4175 if (trim(segment%field(m)%name) ==
'U' .or. trim(segment%field(m)%name) ==
'V')
then 4176 if (trim(segment%field(m)%name) ==
'U' .and. segment%is_E_or_W)
then 4178 do j=js_obc+1,je_obc
4179 normal_trans_bt(i,j) = 0.0
4181 if (obc%add_tide_constituents)
then 4182 do c=1,obc%n_tide_constituents
4183 tidal_vel = tidal_vel + (obc%tide_fn(c) * segment%field(segment%uamp_index)%buffer_dst(i,j,c)) * &
4184 cos((time_delta*obc%tide_frequencies(c) - segment%field(segment%uphase_index)%buffer_dst(i,j,c)) &
4185 + (obc%tide_eq_phases(c) + obc%tide_un(c)))
4189 segment%normal_vel(i,j,k) = us%m_s_to_L_T*(segment%field(m)%buffer_dst(i,j,k) + tidal_vel)
4190 segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k)*segment%h(i,j,k) * g%dyCu(i,j)
4191 normal_trans_bt(i,j) = normal_trans_bt(i,j) + segment%normal_trans(i,j,k)
4193 segment%normal_vel_bt(i,j) = normal_trans_bt(i,j) &
4194 / (max(segment%Htot(i,j), 1.e-12 * gv%m_to_H) * g%dyCu(i,j))
4195 if (
associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,j,:) = segment%normal_vel(i,j,:)
4197 elseif (trim(segment%field(m)%name) ==
'V' .and. segment%is_N_or_S)
then 4199 do i=is_obc+1,ie_obc
4200 normal_trans_bt(i,j) = 0.0
4202 if (obc%add_tide_constituents)
then 4203 do c=1,obc%n_tide_constituents
4204 tidal_vel = tidal_vel + (obc%tide_fn(c) * segment%field(segment%vamp_index)%buffer_dst(i,j,c)) * &
4205 cos((time_delta*obc%tide_frequencies(c) - segment%field(segment%vphase_index)%buffer_dst(i,j,c)) &
4206 + (obc%tide_eq_phases(c) + obc%tide_un(c)))
4210 segment%normal_vel(i,j,k) = us%m_s_to_L_T*(segment%field(m)%buffer_dst(i,j,k) + tidal_vel)
4211 segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k)*segment%h(i,j,k) * &
4213 normal_trans_bt(i,j) = normal_trans_bt(i,j) + segment%normal_trans(i,j,k)
4215 segment%normal_vel_bt(i,j) = normal_trans_bt(i,j) &
4216 / (max(segment%Htot(i,j), 1.e-12 * gv%m_to_H) * g%dxCv(i,j))
4217 if (
associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,j,:) = segment%normal_vel(i,j,:)
4219 elseif (trim(segment%field(m)%name) ==
'V' .and. segment%is_E_or_W .and. &
4220 associated(segment%tangential_vel))
then 4224 if (obc%add_tide_constituents)
then 4225 do c=1,obc%n_tide_constituents
4226 tidal_vel = tidal_vel + (obc%tide_fn(c) * segment%field(segment%vamp_index)%buffer_dst(i,j,c)) * &
4227 cos((time_delta*obc%tide_frequencies(c) - segment%field(segment%vphase_index)%buffer_dst(i,j,c)) &
4228 + (obc%tide_eq_phases(c) + obc%tide_un(c)))
4232 segment%tangential_vel(i,j,k) = us%m_s_to_L_T*(segment%field(m)%buffer_dst(i,j,k) + tidal_vel)
4234 if (
associated(segment%nudged_tangential_vel)) &
4235 segment%nudged_tangential_vel(i,j,:) = segment%tangential_vel(i,j,:)
4237 elseif (trim(segment%field(m)%name) ==
'U' .and. segment%is_N_or_S .and. &
4238 associated(segment%tangential_vel))
then 4242 if (obc%add_tide_constituents)
then 4243 do c=1,obc%n_tide_constituents
4244 tidal_vel = tidal_vel + (obc%tide_fn(c) * segment%field(segment%uamp_index)%buffer_dst(i,j,c)) * &
4245 cos((time_delta*obc%tide_frequencies(c) - segment%field(segment%uphase_index)%buffer_dst(i,j,c)) &
4246 + (obc%tide_eq_phases(c) + obc%tide_un(c)))
4250 segment%tangential_vel(i,j,k) = us%m_s_to_L_T*(segment%field(m)%buffer_dst(i,j,k) + tidal_vel)
4252 if (
associated(segment%nudged_tangential_vel)) &
4253 segment%nudged_tangential_vel(i,j,:) = segment%tangential_vel(i,j,:)
4256 elseif (trim(segment%field(m)%name) ==
'DVDX' .and. segment%is_E_or_W .and. &
4257 associated(segment%tangential_grad))
then 4261 segment%tangential_grad(i,j,k) = us%T_to_s*segment%field(m)%buffer_dst(i,j,k)
4262 if (
associated(segment%nudged_tangential_grad)) &
4263 segment%nudged_tangential_grad(i,j,:) = segment%tangential_grad(i,j,:)
4266 elseif (trim(segment%field(m)%name) ==
'DUDY' .and. segment%is_N_or_S .and. &
4267 associated(segment%tangential_grad))
then 4271 segment%tangential_grad(i,j,k) = us%T_to_s*segment%field(m)%buffer_dst(i,j,k)
4272 if (
associated(segment%nudged_tangential_grad)) &
4273 segment%nudged_tangential_grad(i,j,:) = segment%tangential_grad(i,j,:)
4282 if (segment%is_E_or_W)
then 4289 if (segment%is_N_or_S)
then 4297 if (trim(segment%field(m)%name) ==
'SSH')
then 4302 if (obc%add_tide_constituents)
then 4303 do c=1,obc%n_tide_constituents
4304 tidal_elev = tidal_elev + (obc%tide_fn(c) * segment%field(segment%zamp_index)%buffer_dst(i,j,c)) * &
4305 cos((time_delta*obc%tide_frequencies(c) - segment%field(segment%zphase_index)%buffer_dst(i,j,c)) &
4306 + (obc%tide_eq_phases(c) + obc%tide_un(c)))
4309 segment%eta(i,j) = gv%m_to_H * obc%ramp_value &
4310 * (segment%field(m)%buffer_dst(i,j,1) + tidal_elev)
4317 if (obc%add_tide_constituents)
then 4318 do c=1,obc%n_tide_constituents
4319 tidal_elev = tidal_elev + (obc%tide_fn(c) * segment%field(segment%zamp_index)%buffer_dst(i,j,c)) * &
4320 cos((time_delta*obc%tide_frequencies(c) - segment%field(segment%zphase_index)%buffer_dst(i,j,c)) &
4321 + (obc%tide_eq_phases(c) + obc%tide_un(c)))
4324 segment%eta(i,j) = gv%m_to_H * (segment%field(m)%buffer_dst(i,j,1) + tidal_elev)
4330 if (trim(segment%field(m)%name) ==
'TEMP')
then 4331 if (
associated(segment%field(m)%buffer_dst))
then 4332 do k=1,nz;
do j=js_obc2, je_obc;
do i=is_obc2,ie_obc
4333 segment%tr_Reg%Tr(1)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
4334 enddo ;
enddo ;
enddo 4335 if (.not. segment%tr_Reg%Tr(1)%is_initialized)
then 4337 do k=1,nz;
do j=js_obc2, je_obc;
do i=is_obc2,ie_obc
4338 segment%tr_Reg%Tr(1)%tres(i,j,k) = segment%tr_Reg%Tr(1)%t(i,j,k)
4339 enddo ;
enddo ;
enddo 4340 segment%tr_Reg%Tr(1)%is_initialized=.true.
4343 segment%tr_Reg%Tr(1)%OBC_inflow_conc = segment%field(m)%value
4345 elseif (trim(segment%field(m)%name) ==
'SALT')
then 4346 if (
associated(segment%field(m)%buffer_dst))
then 4347 do k=1,nz;
do j=js_obc2, je_obc;
do i=is_obc2,ie_obc
4348 segment%tr_Reg%Tr(2)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
4349 enddo ;
enddo ;
enddo 4350 if (.not. segment%tr_Reg%Tr(2)%is_initialized)
then 4352 do k=1,nz;
do j=js_obc2, je_obc;
do i=is_obc2,ie_obc
4353 segment%tr_Reg%Tr(2)%tres(i,j,k) = segment%tr_Reg%Tr(2)%t(i,j,k)
4354 enddo ;
enddo ;
enddo 4355 segment%tr_Reg%Tr(2)%is_initialized=.true.
4358 segment%tr_Reg%Tr(2)%OBC_inflow_conc = segment%field(m)%value
4364 deallocate(normal_trans_bt)
4368 end subroutine update_obc_segment_data
4373 subroutine update_obc_ramp(Time, OBC, activate)
4374 type(time_type),
target,
intent(in) :: time
4376 logical,
optional,
intent(in) :: activate
4380 real :: deltatime, wghta
4381 character(len=12) :: msg
4383 if (.not. obc%ramp)
return 4387 if (
present(activate))
then 4389 obc%ramp_start_time = time
4390 obc%ramping_is_activated = .true.
4391 obc%trunc_ramp_time = obc%ramp_timescale
4394 if (.not.obc%ramping_is_activated)
return 4395 deltatime = max( 0., time_type_to_real( time - obc%ramp_start_time ) )
4396 if (deltatime >= obc%trunc_ramp_time)
then 4397 obc%ramp_value = 1.0
4400 wghta = min( 1., deltatime / obc%ramp_timescale )
4405 obc%ramp_value = wghta
4407 write(msg(1:12),
'(es12.3)') obc%ramp_value
4408 call mom_error(note,
"MOM_open_boundary: update_OBC_ramp set OBC"// &
4409 " ramp to "//trim(msg))
4410 end subroutine update_obc_ramp
4413 subroutine register_obc(name, param_file, Reg)
4414 character(len=32),
intent(in) :: name
4418 character(len=256) :: mesg
4420 if (.not.
associated(reg))
call obc_registry_init(param_file, reg)
4422 if (reg%nobc>=max_fields_)
then 4423 write(mesg,
'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for & 4424 &all the open boundaries being registered via register_OBC.")') reg%nobc+1
4425 call mom_error(fatal,
"MOM register_tracer: "//mesg)
4427 reg%nobc = reg%nobc + 1
4430 reg%OB(nobc)%name = name
4432 if (reg%locked)
call mom_error(fatal, &
4433 "MOM register_tracer was called for variable "//trim(reg%OB(nobc)%name)//&
4434 " with a locked tracer registry.")
4436 end subroutine register_obc
4439 subroutine obc_registry_init(param_file, Reg)
4443 integer,
save :: init_calls = 0
4445 #include "version_variable.h" 4446 character(len=40) :: mdl =
"MOM_open_boundary" 4447 character(len=256) :: mesg
4449 if (.not.
associated(reg))
then ;
allocate(reg)
4450 else ;
return ;
endif 4455 init_calls = init_calls + 1
4456 if (init_calls > 1)
then 4457 write(mesg,
'("OBC_registry_init called ",I3, & 4458 &" times with different registry pointers.")') init_calls
4459 if (is_root_pe())
call mom_error(warning,
"MOM_open_boundary"//mesg)
4462 end subroutine obc_registry_init
4465 function register_file_obc(param_file, CS, OBC_Reg)
4469 logical :: register_file_obc
4470 character(len=32) :: casename =
"OBC file" 4472 if (
associated(cs))
then 4473 call mom_error(warning,
"register_file_OBC called with an "// &
4474 "associated control structure.")
4480 call register_obc(casename, param_file, obc_reg)
4481 register_file_obc = .true.
4483 end function register_file_obc
4486 subroutine file_obc_end(CS)
4489 if (
associated(cs))
then 4492 end subroutine file_obc_end
4495 subroutine segment_tracer_registry_init(param_file, segment)
4499 integer,
save :: init_calls = 0
4502 #include "version_variable.h" 4503 character(len=40) :: mdl =
"segment_tracer_registry_init" 4504 character(len=256) :: mesg
4506 if (.not.
associated(segment%tr_Reg))
then 4507 allocate(segment%tr_Reg)
4512 init_calls = init_calls + 1
4515 if (init_calls == 1)
call log_version(param_file, mdl, version,
"")
4524 end subroutine segment_tracer_registry_init
4526 subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, &
4527 OBC_scalar, OBC_array)
4538 real,
optional,
intent(in) :: obc_scalar
4540 logical,
optional,
intent(in) :: obc_array
4546 integer :: isd, ied, jsd, jed
4547 integer :: isdb, iedb, jsdb, jedb
4548 character(len=256) :: mesg
4550 call segment_tracer_registry_init(param_file, segment)
4552 if (segment%tr_Reg%ntseg>=max_fields_)
then 4553 write(mesg,
'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for & 4554 &all the tracers being registered via register_tracer.")') segment%tr_Reg%ntseg+1
4555 call mom_error(fatal,
"MOM register_tracer: "//mesg)
4557 segment%tr_Reg%ntseg = segment%tr_Reg%ntseg + 1
4558 ntseg = segment%tr_Reg%ntseg
4560 isd = segment%HI%isd ; ied = segment%HI%ied
4561 jsd = segment%HI%jsd ; jed = segment%HI%jed
4562 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
4563 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
4565 segment%tr_Reg%Tr(ntseg)%Tr => tr_ptr
4566 segment%tr_Reg%Tr(ntseg)%name = tr_ptr%name
4568 if (segment%tr_Reg%locked)
call mom_error(fatal, &
4569 "MOM register_tracer was called for variable "//trim(segment%tr_Reg%Tr(ntseg)%name)//&
4570 " with a locked tracer registry.")
4572 if (
present(obc_scalar)) segment%tr_Reg%Tr(ntseg)%OBC_inflow_conc = obc_scalar
4573 if (
present(obc_array))
then 4574 if (segment%is_E_or_W)
then 4575 allocate(segment%tr_Reg%Tr(ntseg)%t(isdb:iedb,jsd:jed,1:gv%ke));segment%tr_Reg%Tr(ntseg)%t(:,:,:)=0.0
4576 allocate(segment%tr_Reg%Tr(ntseg)%tres(isdb:iedb,jsd:jed,1:gv%ke));segment%tr_Reg%Tr(ntseg)%tres(:,:,:)=0.0
4577 segment%tr_Reg%Tr(ntseg)%is_initialized=.false.
4578 elseif (segment%is_N_or_S)
then 4579 allocate(segment%tr_Reg%Tr(ntseg)%t(isd:ied,jsdb:jedb,1:gv%ke));segment%tr_Reg%Tr(ntseg)%t(:,:,:)=0.0
4580 allocate(segment%tr_Reg%Tr(ntseg)%tres(isd:ied,jsdb:jedb,1:gv%ke));segment%tr_Reg%Tr(ntseg)%tres(:,:,:)=0.0
4581 segment%tr_Reg%Tr(ntseg)%is_initialized=.false.
4585 end subroutine register_segment_tracer
4588 subroutine segment_tracer_registry_end(Reg)
4594 if (
associated(reg))
then 4596 if (
associated(reg%Tr(n)%t))
deallocate(reg%Tr(n)%t)
4600 end subroutine segment_tracer_registry_end
4602 subroutine register_temp_salt_segments(GV, OBC, tr_Reg, param_file)
4609 integer :: isd, ied, isdb, iedb, jsd, jed, jsdb, jedb, nz, nf
4610 integer :: i, j, k, n
4611 character(len=32) :: name
4615 if (.not.
associated(obc))
return 4617 do n=1, obc%number_of_segments
4618 segment=>obc%segment(n)
4619 if (.not. segment%on_pe) cycle
4621 if (
associated(segment%tr_Reg)) &
4622 call mom_error(fatal,
"register_temp_salt_segments: tracer array was previously allocated")
4625 call tracer_name_lookup(tr_reg, tr_ptr, name)
4626 call register_segment_tracer(tr_ptr, param_file, gv, segment, &
4627 obc_array=segment%temp_segment_data_exists)
4629 call tracer_name_lookup(tr_reg, tr_ptr, name)
4630 call register_segment_tracer(tr_ptr, param_file, gv, segment, &
4631 obc_array=segment%salt_segment_data_exists)
4634 end subroutine register_temp_salt_segments
4636 subroutine fill_temp_salt_segments(G, OBC, tv)
4642 integer :: isd, ied, isdb, iedb, jsd, jed, jsdb, jedb, n, nz
4646 if (.not.
associated(obc))
return 4647 if (.not.
associated(tv%T) .and.
associated(tv%S))
return 4655 do n=1, obc%number_of_segments
4656 segment => obc%segment(n)
4657 if (.not. segment%on_pe) cycle
4659 isd = segment%HI%isd ; ied = segment%HI%ied
4660 jsd = segment%HI%jsd ; jed = segment%HI%jed
4661 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
4662 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
4665 if (segment%is_E_or_W)
then 4667 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
4668 if (segment%direction == obc_direction_w)
then 4669 segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i+1,j,k)
4670 segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i+1,j,k)
4672 segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j,k)
4673 segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j,k)
4678 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
4679 if (segment%direction == obc_direction_s)
then 4680 segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j+1,k)
4681 segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j+1,k)
4683 segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j,k)
4684 segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j,k)
4688 segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:)
4689 segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:)
4692 call setup_obc_tracer_reservoirs(g, obc)
4693 end subroutine fill_temp_salt_segments
4698 subroutine mask_outside_obcs(G, US, param_file, OBC)
4705 integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n
4708 logical :: fatal_error = .false.
4710 integer,
parameter :: cin = 3, cout = 4, cland = -1, cedge = -2
4711 character(len=256) :: mesg
4713 real,
allocatable,
dimension(:,:) :: color, color2
4716 if (.not.
associated(obc))
return 4718 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
4719 units=
"m", default=0.0, scale=us%m_to_Z, do_not_log=.true.)
4721 allocate(color(g%isd:g%ied, g%jsd:g%jed)) ; color = 0
4722 allocate(color2(g%isd:g%ied, g%jsd:g%jed)) ; color2 = 0
4727 color(g%isd,j) = cedge
4728 color(g%ied,j) = cedge
4729 color2(g%isd,j) = cedge
4730 color2(g%ied,j) = cedge
4733 color(i,g%jsd) = cedge
4734 color(i,g%jed) = cedge
4735 color2(i,g%jsd) = cedge
4736 color2(i,g%jed) = cedge
4743 if (g%bathyT(i,j) <= min_depth) color(i,j) = cland
4744 if (g%bathyT(i,j) <= min_depth) color2(i,j) = cland
4748 do j=g%jsd,g%jed ;
do i=g%IsdB+1,g%IedB-1
4749 l_seg = obc%segnum_u(i,j)
4750 if (l_seg == obc_none) cycle
4752 if (obc%segment(l_seg)%direction == obc_direction_w)
then 4753 if (color(i,j) == 0.0) color(i,j) = cout
4754 if (color(i+1,j) == 0.0) color(i+1,j) = cin
4755 elseif (obc%segment(l_seg)%direction == obc_direction_e)
then 4756 if (color(i,j) == 0.0) color(i,j) = cin
4757 if (color(i+1,j) == 0.0) color(i+1,j) = cout
4760 do j=g%JsdB+1,g%JedB-1 ;
do i=g%isd,g%ied
4761 l_seg = obc%segnum_v(i,j)
4762 if (l_seg == obc_none) cycle
4764 if (obc%segment(l_seg)%direction == obc_direction_s)
then 4765 if (color(i,j) == 0.0) color(i,j) = cout
4766 if (color(i,j+1) == 0.0) color(i,j+1) = cin
4767 elseif (obc%segment(l_seg)%direction == obc_direction_n)
then 4768 if (color(i,j) == 0.0) color(i,j) = cin
4769 if (color(i,j+1) == 0.0) color(i,j+1) = cout
4773 do j=g%JsdB+1,g%JedB-1 ;
do i=g%isd,g%ied
4774 l_seg = obc%segnum_v(i,j)
4775 if (l_seg == obc_none) cycle
4777 if (obc%segment(l_seg)%direction == obc_direction_s)
then 4778 if (color2(i,j) == 0.0) color2(i,j) = cout
4779 if (color2(i,j+1) == 0.0) color2(i,j+1) = cin
4780 elseif (obc%segment(l_seg)%direction == obc_direction_n)
then 4781 if (color2(i,j) == 0.0) color2(i,j) = cin
4782 if (color2(i,j+1) == 0.0) color2(i,j+1) = cout
4785 do j=g%jsd,g%jed ;
do i=g%IsdB+1,g%IedB-1
4786 l_seg = obc%segnum_u(i,j)
4787 if (l_seg == obc_none) cycle
4789 if (obc%segment(l_seg)%direction == obc_direction_w)
then 4790 if (color2(i,j) == 0.0) color2(i,j) = cout
4791 if (color2(i+1,j) == 0.0) color2(i+1,j) = cin
4792 elseif (obc%segment(l_seg)%direction == obc_direction_e)
then 4793 if (color2(i,j) == 0.0) color2(i,j) = cin
4794 if (color2(i+1,j) == 0.0) color2(i+1,j) = cout
4799 call flood_fill(g, color, cin, cout, cland)
4800 call flood_fill2(g, color2, cin, cout, cland)
4803 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
4804 if (color(i,j) /= color2(i,j))
then 4805 fatal_error = .true.
4806 write(mesg,
'("MOM_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", & 4807 "the masking of the outside grid points.")') i, j
4808 call mom_error(warning,
"MOM register_tracer: "//mesg, all_print=.true.)
4810 if (color(i,j) == cout) g%bathyT(i,j) = min_depth
4812 if (fatal_error)
call mom_error(fatal, &
4813 "MOM_open_boundary: inconsistent OBC segments.")
4817 end subroutine mask_outside_obcs
4820 subroutine flood_fill(G, color, cin, cout, cland)
4822 real,
dimension(:,:),
intent(inout) :: color
4823 integer,
intent(in) :: cin
4824 integer,
intent(in) :: cout
4825 integer,
intent(in) :: cland
4828 integer :: i, j, ncount
4831 do while (ncount > 0)
4833 do j=g%jsd+1,g%jed-1
4834 do i=g%isd+1,g%ied-1
4835 if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0)
then 4836 color(i,j) = color(i-1,j)
4839 if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0)
then 4840 color(i,j) = color(i+1,j)
4843 if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0)
then 4844 color(i,j) = color(i,j-1)
4847 if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0)
then 4848 color(i,j) = color(i,j+1)
4853 do j=g%jed-1,g%jsd+1,-1
4854 do i=g%ied-1,g%isd+1,-1
4855 if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0)
then 4856 color(i,j) = color(i-1,j)
4859 if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0)
then 4860 color(i,j) = color(i+1,j)
4863 if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0)
then 4864 color(i,j) = color(i,j-1)
4867 if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0)
then 4868 color(i,j) = color(i,j+1)
4874 call sum_across_pes(ncount)
4877 end subroutine flood_fill
4880 subroutine flood_fill2(G, color, cin, cout, cland)
4882 real,
dimension(:,:),
intent(inout) :: color
4883 integer,
intent(in) :: cin
4884 integer,
intent(in) :: cout
4885 integer,
intent(in) :: cland
4888 integer :: i, j, ncount
4891 do while (ncount > 0)
4893 do i=g%isd+1,g%ied-1
4894 do j=g%jsd+1,g%jed-1
4895 if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0)
then 4896 color(i,j) = color(i-1,j)
4899 if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0)
then 4900 color(i,j) = color(i+1,j)
4903 if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0)
then 4904 color(i,j) = color(i,j-1)
4907 if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0)
then 4908 color(i,j) = color(i,j+1)
4913 do i=g%ied-1,g%isd+1,-1
4914 do j=g%jed-1,g%jsd+1,-1
4915 if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0)
then 4916 color(i,j) = color(i-1,j)
4919 if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0)
then 4920 color(i,j) = color(i+1,j)
4923 if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0)
then 4924 color(i,j) = color(i,j-1)
4927 if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0)
then 4928 color(i,j) = color(i,j+1)
4934 call sum_across_pes(ncount)
4937 end subroutine flood_fill2
4940 subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart_CSp, &
4942 type(hor_index_type),
intent(in) :: hi
4948 logical,
intent(in) :: use_temperature
4952 character(len=100) :: mesg
4955 if (.not.
associated(obc)) &
4956 call mom_error(fatal,
"open_boundary_register_restarts: Called with "//&
4957 "uninitialized OBC control structure")
4959 if (
associated(obc%rx_normal) .or.
associated(obc%ry_normal) .or. &
4960 associated(obc%rx_oblique) .or.
associated(obc%ry_oblique) .or.
associated(obc%cff_normal)) &
4961 call mom_error(fatal,
"open_boundary_register_restarts: Restart "//&
4962 "arrays were previously allocated")
4964 if (
associated(obc%tres_x) .or.
associated(obc%tres_y)) &
4965 call mom_error(fatal,
"open_boundary_register_restarts: Restart "//&
4966 "arrays were previously allocated")
4972 if (obc%radiation_BCs_exist_globally)
then 4973 allocate(obc%rx_normal(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke))
4974 allocate(obc%ry_normal(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke))
4975 obc%rx_normal(:,:,:) = 0.0
4976 obc%ry_normal(:,:,:) = 0.0
4978 vd(1) = var_desc(
"rx_normal",
"m s-1",
"Normal Phase Speed for EW radiation OBCs",
'u',
'L')
4979 vd(2) = var_desc(
"ry_normal",
"m s-1",
"Normal Phase Speed for NS radiation OBCs",
'v',
'L')
4981 .false., restart_csp)
4984 if (obc%oblique_BCs_exist_globally)
then 4985 allocate(obc%rx_oblique(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke))
4986 allocate(obc%ry_oblique(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke))
4987 obc%rx_oblique(:,:,:) = 0.0
4988 obc%ry_oblique(:,:,:) = 0.0
4990 vd(1) = var_desc(
"rx_oblique",
"m2 s-2",
"Radiation Speed Squared for EW oblique OBCs",
'u',
'L')
4991 vd(2) = var_desc(
"ry_oblique",
"m2 s-2",
"Radiation Speed Squared for NS oblique OBCs",
'v',
'L')
4993 .false., restart_csp)
4995 allocate(obc%cff_normal(hi%IsdB:hi%IedB,hi%jsdB:hi%jedB,gv%ke))
4996 obc%cff_normal(:,:,:) = 0.0
4997 vd(1) = var_desc(
"cff_normal",
"m2 s-2",
"denominator for oblique OBCs",
'q',
'L')
5001 if (reg%ntr == 0)
return 5002 if (.not.
associated(obc%tracer_x_reservoirs_used))
then 5004 allocate(obc%tracer_x_reservoirs_used(reg%ntr))
5005 allocate(obc%tracer_y_reservoirs_used(reg%ntr))
5006 obc%tracer_x_reservoirs_used(:) = .false.
5007 obc%tracer_y_reservoirs_used(:) = .false.
5008 call parse_for_tracer_reservoirs(obc, param_file, use_temperature)
5011 if (obc%ntr /= reg%ntr)
then 5013 write(mesg,
'("Inconsistent values for ntr ", I8," and ",I8,".")') obc%ntr, reg%ntr
5014 call mom_error(warning,
'open_boundary_register_restarts: '//mesg)
5019 if (any(obc%tracer_x_reservoirs_used))
then 5020 allocate(obc%tres_x(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke,obc%ntr))
5021 obc%tres_x(:,:,:,:) = 0.0
5023 if (obc%tracer_x_reservoirs_used(m))
then 5024 if (modulo(hi%turns, 2) /= 0)
then 5025 write(mesg,
'("tres_y_",I3.3)') m
5026 vd(1) = var_desc(mesg,
"Conc",
"Tracer concentration for NS OBCs",
'v',
'L')
5029 write(mesg,
'("tres_x_",I3.3)') m
5030 vd(1) = var_desc(mesg,
"Conc",
"Tracer concentration for EW OBCs",
'u',
'L')
5036 if (any(obc%tracer_y_reservoirs_used))
then 5037 allocate(obc%tres_y(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke,obc%ntr))
5038 obc%tres_y(:,:,:,:) = 0.0
5040 if (obc%tracer_y_reservoirs_used(m))
then 5041 if (modulo(hi%turns, 2) /= 0)
then 5042 write(mesg,
'("tres_x_",I3.3)') m
5043 vd(1) = var_desc(mesg,
"Conc",
"Tracer concentration for EW OBCs",
'u',
'L')
5046 write(mesg,
'("tres_y_",I3.3)') m
5047 vd(1) = var_desc(mesg,
"Conc",
"Tracer concentration for NS OBCs",
'v',
'L')
5053 end subroutine open_boundary_register_restarts
5056 subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
5059 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: uhr
5061 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: vhr
5063 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
5066 real,
intent(in) :: dt
5070 real :: u_l_in, u_l_out
5071 real :: v_l_in, v_l_out
5073 integer :: i, j, k, m, n, ntr, nz
5074 integer :: ishift, idir, jshift, jdir
5078 if (
associated(obc))
then ;
if (obc%OBC_pe)
then ;
do n=1,obc%number_of_segments
5079 segment=>obc%segment(n)
5080 if (.not.
associated(segment%tr_Reg)) cycle
5081 if (segment%is_E_or_W)
then 5083 do j=segment%HI%jsd,segment%HI%jed
5086 if (segment%direction == obc_direction_w)
then 5087 ishift = 1 ; idir = -1
5089 ishift = 0 ; idir = 1
5092 if (g%mask2dT(i+ishift,j) == 0.0) cycle
5094 do m=1,ntr ;
if (
associated(segment%tr_Reg%Tr(m)%tres))
then ;
do k=1,nz
5095 u_l_out = max(0.0, (idir*uhr(i,j,k))*segment%Tr_InvLscale_out / &
5096 ((h(i+ishift,j,k) + gv%H_subroundoff)*g%dyCu(i,j)))
5097 u_l_in = min(0.0, (idir*uhr(i,j,k))*segment%Tr_InvLscale_in / &
5098 ((h(i+ishift,j,k) + gv%H_subroundoff)*g%dyCu(i,j)))
5099 fac1 = 1.0 + (u_l_out-u_l_in)
5100 segment%tr_Reg%Tr(m)%tres(i,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,j,k) + &
5101 (u_l_out*reg%Tr(m)%t(i+ishift,j,k) - &
5102 u_l_in*segment%tr_Reg%Tr(m)%t(i,j,k)))
5103 if (
associated(obc%tres_x)) obc%tres_x(i,j,k,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
5104 enddo ;
endif ;
enddo 5106 elseif (segment%is_N_or_S)
then 5108 do i=segment%HI%isd,segment%HI%ied
5111 if (segment%direction == obc_direction_s)
then 5112 jshift = 1 ; jdir = -1
5114 jshift = 0 ; jdir = 1
5117 if (g%mask2dT(i,j+jshift) == 0.0) cycle
5119 do m=1,ntr ;
if (
associated(segment%tr_Reg%Tr(m)%tres))
then ;
do k=1,nz
5120 v_l_out = max(0.0, (jdir*vhr(i,j,k))*segment%Tr_InvLscale_out / &
5121 ((h(i,j+jshift,k) + gv%H_subroundoff)*g%dxCv(i,j)))
5122 v_l_in = min(0.0, (jdir*vhr(i,j,k))*segment%Tr_InvLscale_in / &
5123 ((h(i,j+jshift,k) + gv%H_subroundoff)*g%dxCv(i,j)))
5124 fac1 = 1.0 + (v_l_out-v_l_in)
5125 segment%tr_Reg%Tr(m)%tres(i,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,j,k) + &
5126 (v_l_out*reg%Tr(m)%t(i,j+jshift,k) - &
5127 v_l_in*segment%tr_Reg%Tr(m)%t(i,j,k)))
5128 if (
associated(obc%tres_y)) obc%tres_y(i,j,k,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
5129 enddo ;
endif ;
enddo 5132 enddo ;
endif ;
endif 5134 end subroutine update_segment_tracer_reservoirs
5144 subroutine adjustsegmentetatofitbathymetry(G, GV, US, segment,fld)
5149 integer,
intent(in) :: fld
5151 integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
5153 real,
allocatable,
dimension(:,:,:) :: eta
5154 real :: hTolerance = 0.1
5155 real :: hTmp, eTmp, dilate
5156 character(len=100) :: mesg
5158 htolerance = 0.1*us%m_to_Z
5160 nz =
size(segment%field(fld)%dz_src,3)
5162 if (segment%is_E_or_W)
then 5164 is = segment%HI%isdB ; ie = segment%HI%iedB
5165 js = segment%HI%jsd ; je = segment%HI%jed
5167 is = segment%HI%isd ; ie = segment%HI%ied
5168 js = segment%HI%jsdB ; je = segment%HI%jedB
5170 allocate(eta(is:ie,js:je,nz+1))
5171 contractions=0; dilations=0
5172 do j=js,je ;
do i=is,ie
5180 eta(i,j,k)=eta(i,j,k-1)-segment%field(fld)%dz_src(i,j,k-1)
5185 if (-eta(i,j,k) > segment%Htot(i,j)*gv%H_to_Z + htolerance)
then 5186 eta(i,j,k) = -segment%Htot(i,j)*gv%H_to_Z
5187 contractions = contractions + 1
5194 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z))
then 5195 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
5196 segment%field(fld)%dz_src(i,j,k) = gv%Angstrom_Z
5198 segment%field(fld)%dz_src(i,j,k) = (eta(i,j,k) - eta(i,j,k+1))
5204 if (-eta(i,j,nz+1) < (segment%Htot(i,j) * gv%H_to_Z) - htolerance)
then 5205 dilations = dilations + 1
5207 eta(i,j,nz+1) = -(segment%Htot(i,j) * gv%H_to_Z)
5208 segment%field(fld)%dz_src(i,j,nz)= eta(i,j,nz)-eta(i,j,nz+1)
5219 segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k)*gv%Z_to_H
5241 end subroutine adjustsegmentetatofitbathymetry
5244 subroutine rotate_obc_config(OBC_in, G_in, OBC, G, turns)
5246 type(dyn_horgrid_type),
intent(in) :: g_in
5248 type(dyn_horgrid_type),
intent(in) :: g
5249 integer,
intent(in) :: turns
5253 if (obc_in%number_of_segments==0)
return 5256 obc%number_of_segments = obc_in%number_of_segments
5258 obc%user_BCs_set_globally = obc_in%user_BCs_set_globally
5261 obc%zero_vorticity = obc_in%zero_vorticity
5262 obc%freeslip_vorticity = obc_in%freeslip_vorticity
5263 obc%computed_vorticity = obc_in%computed_vorticity
5264 obc%specified_vorticity = obc_in%specified_vorticity
5265 obc%zero_strain = obc_in%zero_strain
5266 obc%freeslip_strain = obc_in%freeslip_strain
5267 obc%computed_strain = obc_in%computed_strain
5268 obc%specified_strain = obc_in%specified_strain
5269 obc%zero_biharmonic = obc_in%zero_biharmonic
5270 obc%silly_h = obc_in%silly_h
5271 obc%silly_u = obc_in%silly_u
5274 allocate(obc%segment(0:obc%number_of_segments))
5275 do l = 1, obc%number_of_segments
5276 call rotate_obc_segment_config(obc_in%segment(l), g_in, obc%segment(l), g, turns)
5278 call allocate_obc_segment_data(obc, obc%segment(l))
5279 call rotate_obc_segment_data(obc_in%segment(l), obc%segment(l), turns)
5283 allocate(obc%segnum_u(g%IsdB:g%IedB,g%jsd:g%jed))
5284 allocate(obc%segnum_v(g%isd:g%ied,g%JsdB:g%JedB))
5285 call rotate_array_pair(obc_in%segnum_u, obc_in%segnum_v, turns, &
5286 obc%segnum_u, obc%segnum_v)
5289 obc%open_u_BCs_exist_globally = obc_in%open_v_BCs_exist_globally
5290 obc%open_v_BCs_exist_globally = obc_in%open_u_BCs_exist_globally
5291 obc%Flather_u_BCs_exist_globally = obc_in%Flather_v_BCs_exist_globally
5292 obc%Flather_v_BCs_exist_globally = obc_in%Flather_u_BCs_exist_globally
5293 obc%oblique_BCs_exist_globally = obc_in%oblique_BCs_exist_globally
5294 obc%nudged_u_BCs_exist_globally = obc_in%nudged_v_BCs_exist_globally
5295 obc%nudged_v_BCs_exist_globally = obc_in%nudged_u_BCs_exist_globally
5296 obc%specified_u_BCs_exist_globally= obc_in%specified_v_BCs_exist_globally
5297 obc%specified_v_BCs_exist_globally= obc_in%specified_u_BCs_exist_globally
5298 obc%radiation_BCs_exist_globally = obc_in%radiation_BCs_exist_globally
5301 obc%brushcutter_mode = obc_in%brushcutter_mode
5302 obc%update_OBC = obc_in%update_OBC
5303 obc%needs_IO_for_data = obc_in%needs_IO_for_data
5305 obc%ntr = obc_in%ntr
5307 obc%gamma_uv = obc_in%gamma_uv
5308 obc%rx_max = obc_in%rx_max
5309 obc%OBC_pe = obc_in%OBC_pe
5312 if (
ASSOCIATED(obc_in%remap_CS))
then 5313 allocate(obc%remap_CS)
5314 obc%remap_CS = obc_in%remap_CS
5320 end subroutine rotate_obc_config
5323 subroutine rotate_obc_segment_config(segment_in, G_in, segment, G, turns)
5325 type(dyn_horgrid_type),
intent(in) :: G_in
5327 type(dyn_horgrid_type),
intent(in) :: G
5328 integer,
intent(in) :: turns
5331 integer :: Is_obc_in, Ie_obc_in, Js_obc_in, Je_obc_in
5332 integer :: Is_obc, Ie_obc, Js_obc, Je_obc
5339 segment%on_pe = segment_in%on_pe
5342 segment%Flather = segment_in%Flather
5343 segment%radiation = segment_in%radiation
5344 segment%radiation_tan = segment_in%radiation_tan
5345 segment%radiation_grad = segment_in%radiation_grad
5346 segment%oblique = segment_in%oblique
5347 segment%oblique_tan = segment_in%oblique_tan
5348 segment%oblique_grad = segment_in%oblique_grad
5349 segment%nudged = segment_in%nudged
5350 segment%nudged_tan = segment_in%nudged_tan
5351 segment%nudged_grad = segment_in%nudged_grad
5352 segment%specified = segment_in%specified
5353 segment%specified_tan = segment_in%specified_tan
5354 segment%specified_grad = segment_in%specified_grad
5355 segment%open = segment_in%open
5356 segment%gradient = segment_in%gradient
5359 segment%u_values_needed = segment_in%v_values_needed
5360 segment%v_values_needed = segment_in%u_values_needed
5361 segment%z_values_needed = segment_in%z_values_needed
5362 segment%g_values_needed = segment_in%g_values_needed
5363 segment%t_values_needed = segment_in%t_values_needed
5364 segment%s_values_needed = segment_in%s_values_needed
5366 segment%values_needed = segment_in%values_needed
5369 segment%Velocity_nudging_timescale_in = segment_in%Velocity_nudging_timescale_in
5370 segment%Velocity_nudging_timescale_out= segment_in%Velocity_nudging_timescale_out
5379 if (segment_in%direction == obc_direction_n)
then 5380 is_obc_in = segment_in%Ie_obc + g_in%idg_offset
5381 ie_obc_in = segment_in%Is_obc + g_in%idg_offset
5383 is_obc_in = segment_in%Is_obc + g_in%idg_offset
5384 ie_obc_in = segment_in%Ie_obc + g_in%idg_offset
5387 if (segment_in%direction == obc_direction_w)
then 5388 js_obc_in = segment_in%Je_obc + g_in%jdg_offset
5389 je_obc_in = segment_in%Js_obc + g_in%jdg_offset
5391 js_obc_in = segment_in%Js_obc + g_in%jdg_offset
5392 je_obc_in = segment_in%Je_obc + g_in%jdg_offset
5396 is_obc = g_in%jegB - js_obc_in
5397 ie_obc = g_in%JegB - je_obc_in
5404 call setup_segment_indices(g, segment, is_obc, ie_obc, js_obc, je_obc)
5407 if (is_obc > ie_obc)
then 5408 segment%Is_obc = ie_obc - g%idg_offset
5409 segment%Ie_obc = is_obc - g%idg_offset
5411 segment%Is_obc = is_obc - g%idg_offset
5412 segment%Ie_obc = ie_obc - g%idg_offset
5415 if (js_obc > je_obc)
then 5416 segment%Js_obc = je_obc - g%jdg_offset
5417 segment%Je_obc = js_obc - g%jdg_offset
5419 segment%Js_obc = js_obc - g%jdg_offset
5420 segment%Je_obc = je_obc - g%jdg_offset
5425 select case (segment_in%direction)
5426 case (obc_direction_n)
5427 segment%direction = obc_direction_w
5428 segment%is_E_or_W_2 = segment_in%is_N_or_S
5429 segment%is_E_or_W = segment_in%is_N_or_S .and. segment_in%on_pe
5430 segment%is_N_or_S = .false.
5431 case (obc_direction_w)
5432 segment%direction = obc_direction_s
5433 segment%is_N_or_S = segment_in%is_E_or_W
5434 segment%is_E_or_W = .false.
5435 segment%is_E_or_W_2 = .false.
5436 case (obc_direction_s)
5437 segment%direction = obc_direction_e
5438 segment%is_E_or_W_2 = segment_in%is_N_or_S
5439 segment%is_E_or_W = segment_in%is_N_or_S .and. segment_in%on_pe
5440 segment%is_N_or_S = .false.
5441 case (obc_direction_e)
5442 segment%direction = obc_direction_n
5443 segment%is_N_or_S = segment_in%is_E_or_W
5444 segment%is_E_or_W = .false.
5445 segment%is_E_or_W_2 = .false.
5447 segment%direction = obc_none
5451 segment%Tr_InvLscale_in = segment_in%Tr_InvLscale_in
5452 segment%Tr_InvLscale_out = segment_in%Tr_InvLscale_out
5453 end subroutine rotate_obc_segment_config
5457 subroutine rotate_obc_init(OBC_in, G, GV, US, param_file, tv, restart_CSp, OBC)
5459 type(ocean_grid_type),
intent(in) :: g
5460 type(verticalgrid_type),
intent(in) :: gv
5461 type(unit_scale_type),
intent(in) :: us
5462 type(param_file_type),
intent(in) :: param_file
5463 type(thermo_var_ptrs),
intent(inout) :: tv
5464 type(mom_restart_cs),
pointer,
intent(in) :: restart_csp
5467 logical :: use_temperature
5470 call get_param(param_file,
"MOM",
"ENABLE_THERMODYNAMICS", use_temperature, &
5471 "If true, Temperature and salinity are used as state "//&
5472 "variables.", default=.true., do_not_log=.true.)
5474 do l = 1, obc%number_of_segments
5475 call rotate_obc_segment_data(obc_in%segment(l), obc%segment(l), g%HI%turns)
5478 if (use_temperature) &
5479 call fill_temp_salt_segments(g, obc, tv)
5481 call open_boundary_init(g, gv, us, param_file, obc, restart_csp)
5482 end subroutine rotate_obc_init
5486 subroutine rotate_obc_segment_data(segment_in, segment, turns)
5489 integer,
intent(in) :: turns
5492 integer :: is, ie, js, je, nk
5493 integer :: num_fields
5496 num_fields = segment_in%num_fields
5497 allocate(segment%field(num_fields))
5499 segment%num_fields = segment_in%num_fields
5500 do n = 1, num_fields
5501 segment%field(n)%fid = segment_in%field(n)%fid
5502 segment%field(n)%fid_dz = segment_in%field(n)%fid_dz
5504 if (modulo(turns, 2) /= 0)
then 5505 select case (segment_in%field(n)%name)
5507 segment%field(n)%name =
'V' 5509 segment%field(n)%name =
'Vamp' 5511 segment%field(n)%name =
'Vphase' 5513 segment%field(n)%name =
'U' 5515 segment%field(n)%name =
'Uamp' 5517 segment%field(n)%name =
'Uphase' 5519 segment%field(n)%name =
'DUDY' 5521 segment%field(n)%name =
'DVDX' 5523 segment%field(n)%name = segment_in%field(n)%name
5526 segment%field(n)%name = segment_in%field(n)%name
5529 if (
allocated(segment_in%field(n)%buffer_src))
then 5530 call allocate_rotated_array(segment_in%field(n)%buffer_src, &
5531 lbound(segment_in%field(n)%buffer_src), turns, &
5532 segment%field(n)%buffer_src)
5533 call rotate_array(segment_in%field(n)%buffer_src, turns, &
5534 segment%field(n)%buffer_src)
5537 segment%field(n)%nk_src = segment_in%field(n)%nk_src
5539 if (
allocated(segment_in%field(n)%dz_src))
then 5540 call allocate_rotated_array(segment_in%field(n)%dz_src, &
5541 lbound(segment_in%field(n)%dz_src), turns, &
5542 segment%field(n)%dz_src)
5543 call rotate_array(segment_in%field(n)%dz_src, turns, &
5544 segment%field(n)%dz_src)
5547 segment%field(n)%buffer_dst => null()
5548 segment%field(n)%value = segment_in%field(n)%value
5551 segment%temp_segment_data_exists = segment_in%temp_segment_data_exists
5552 segment%salt_segment_data_exists = segment_in%salt_segment_data_exists
5553 end subroutine rotate_obc_segment_data
Open boundary segment data from files (mostly).
Wraps the FMS time manager functions.
Methods for testing for, and list of, obsolete run-time parameters.
Ocean grid type. See mom_grid for details.
A structure that can be parsed to read and document run-time parameters.
Provides the ocean grid type.
Open boundary segment data structure.
Registry type for tracers on segments.
Generates vertical grids as part of the ALE algorithm.
Provides column-wise vertical remapping functions.
Tracer on OBC segment data structure, for putting into a segment tracer registry.
Wraps the MPP cpu clock functions.
Register fields for restarts.
Describes the horizontal ocean grid with only dynamic memory arrays.
This module contains I/O framework code.
The MOM6 facility to parse input files for runtime parameters.
Register a pair of restart fieilds whose rotations map onto each other.
An overloaded interface to log the values of various types of parameters.
Do a halo update on a pair of arrays representing the two components of a vector.
Container for remapping parameters.
Interfaces to non-domain-oriented communication subroutines, including the MOM6 reproducing sums faci...
A restart registry and the control structure for restarts.
Describes various unit conversion factors.
This module contains the tracer_registry_type and the subroutines that handle registration of tracers...
Regridding control structure.
Tidal contributions to geopotential.
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Type to carry basic OBC information needed for updating values.
Type to carry something (what] for the OBC registry.
Describes the decomposed MOM domain and has routines for communications across PEs.
Type to carry basic tracer information.
Routines for error handling and I/O management.
The MOM6 facility for reading and writing restart files, and querying what has been read.
Type for describing a variable, typically a tracer.
An overloaded interface to log version information about modules.
Describes the vertical ocean grid, including unit conversion factors.
Control structure for open boundaries that read from files. Probably lots to update here.
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Simple type to store astronomical longitudes used to calculate tidal phases.
Indicate whether a field has been read from a restart file.
Contains a shareable dynamic type for describing horizontal grids and metric data and utilty routines...
Controls where open boundary conditions are applied.
Handy functions for manipulating strings.
Provides a transparent vertical ocean grid type and supporting routines.
Provides transparent structures with groups of MOM6 variables and supporting routines.
Do a halo update on an array.
An overloaded interface to read and log the values of various types of parameters.