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)
1644 type(param_file_type),
intent(in) :: PF
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
1664 call get_param(pf, mdl, segnam, segstr)
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)
1792 type(ocean_grid_type),
intent(in) :: g
1793 type(verticalgrid_type),
intent(in) :: gv
1794 type(unit_scale_type),
intent(in) :: us
1795 type(param_file_type),
intent(in) :: param_file
1797 type(mom_restart_cs),
pointer :: 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
1851 if (query_initialized(obc%rx_oblique,
"rx_oblique", restart_csp))
then
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
1856 if (query_initialized(obc%ry_oblique,
"ry_oblique", restart_csp))
then
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
1861 if (query_initialized(obc%cff_normal,
"cff_normal", restart_csp))
then
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)
1926 type(dyn_horgrid_type),
intent(in) :: g
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)
1971 type(dyn_horgrid_type),
intent(inout) :: g
1972 type(unit_scale_type),
intent(in) :: 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)
2080 type(ocean_grid_type),
intent(in) :: G
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)
2119 type(ocean_grid_type),
intent(inout) :: g
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
2129 type(unit_scale_type),
intent(in) :: us
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)
3231 type(ocean_grid_type),
intent(inout) :: g
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)
3267 type(ocean_grid_type),
intent(inout) :: g
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)
3297 type(ocean_grid_type),
intent(in) :: G
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)
3421 type(ocean_grid_type),
intent(inout) :: g
3423 type(thermo_var_ptrs),
intent(inout) :: tv
3424 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(inout) :: h
3425 type(param_file_type),
intent(in) :: pf
3426 type(tracer_registry_type),
pointer :: tracer_reg
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
3448 call pass_var(tv%T, g%Domain)
3449 call pass_var(tv%S, g%Domain)
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)
3645 type(ocean_grid_type),
intent(in) :: g
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)
3688 type(ocean_grid_type),
intent(in) :: g
3689 type(verticalgrid_type),
intent(in) :: gv
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)
3732 type(ocean_grid_type),
intent(in) :: g
3733 type(verticalgrid_type),
intent(in) :: gv
3734 type(unit_scale_type),
intent(in) :: us
3736 type(thermo_var_ptrs),
intent(in) :: tv
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,:,:))
3923 call rotate_array(tmp_buffer_in, turns, tmp_buffer)
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,:,:))
3987 call rotate_array(tmp_buffer_in, turns, tmp_buffer)
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
4415 type(param_file_type),
intent(in) :: param_file
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)
4440 type(param_file_type),
intent(in) :: param_file
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)
4466 type(param_file_type),
intent(in) :: param_file
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)
4496 type(param_file_type),
intent(in) :: param_file
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)
4528 type(verticalgrid_type),
intent(in) :: gv
4529 type(tracer_type),
target :: tr_ptr
4536 type(param_file_type),
intent(in) :: param_file
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)
4603 type(verticalgrid_type),
intent(in) :: gv
4605 type(tracer_registry_type),
pointer :: tr_reg
4606 type(param_file_type),
intent(in) :: 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
4613 type(tracer_type),
pointer :: tr_ptr => null()
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)
4637 type(ocean_grid_type),
intent(in) :: g
4639 type(thermo_var_ptrs),
intent(inout) :: 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
4650 call pass_var(tv%T, g%Domain)
4651 call pass_var(tv%S, g%Domain)
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)
4699 type(dyn_horgrid_type),
intent(inout) :: G
4700 type(param_file_type),
intent(in) :: param_file
4702 type(unit_scale_type),
intent(in) :: US
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)
4821 type(dyn_horgrid_type),
intent(inout) :: G
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)
4873 call pass_var(color, g%Domain)
4874 call sum_across_pes(ncount)
4877 end subroutine flood_fill
4880 subroutine flood_fill2(G, color, cin, cout, cland)
4881 type(dyn_horgrid_type),
intent(inout) :: G
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)
4933 call pass_var(color, g%Domain)
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
4943 type(verticalgrid_type),
pointer :: gv
4945 type(tracer_registry_type),
pointer :: reg
4946 type(param_file_type),
intent(in) :: param_file
4947 type(mom_restart_cs),
pointer :: restart_csp
4948 logical,
intent(in) :: use_temperature
4950 type(vardesc) :: vd(2)
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')
4980 call register_restart_pair(obc%rx_normal, obc%ry_normal, vd(1), vd(2), &
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')
4992 call register_restart_pair(obc%rx_oblique, obc%ry_oblique, vd(1), vd(2), &
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')
4998 call register_restart_field(obc%cff_normal, vd(1), .false., restart_csp)
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')
5027 call register_restart_field(obc%tres_x(:,:,:,m), vd(1), .false., restart_csp)
5029 write(mesg,
'("tres_x_",I3.3)') m
5030 vd(1) = var_desc(mesg,
"Conc",
"Tracer concentration for EW OBCs",
'u',
'L')
5031 call register_restart_field(obc%tres_x(:,:,:,m), vd(1), .false., restart_csp)
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')
5044 call register_restart_field(obc%tres_y(:,:,:,m), vd(1), .false., restart_csp)
5046 write(mesg,
'("tres_y_",I3.3)') m
5047 vd(1) = var_desc(mesg,
"Conc",
"Tracer concentration for NS OBCs",
'v',
'L')
5048 call register_restart_field(obc%tres_y(:,:,:,m), vd(1), .false., restart_csp)
5053 end subroutine open_boundary_register_restarts
5056 subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
5057 type(ocean_grid_type),
intent(in) :: g
5058 type(verticalgrid_type),
intent(in) :: gv
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
5067 type(tracer_registry_type),
pointer :: reg
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)
5145 type(ocean_grid_type),
intent(in) :: G
5146 type(verticalgrid_type),
intent(in) :: GV
5147 type(unit_scale_type),
intent(in) :: US
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