7 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
10 use mom_domains,
only : to_all, scalar_pair, cgrid_ne
15 use mom_io,
only : east_face, north_face
16 use mom_io,
only : slasher, read_data, field_size, single_file
19 use mom_obsolete_params,
only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char
22 use time_interp_external_mod,
only : init_external_field, time_interp_external
23 use time_interp_external_mod,
only : time_interp_external_init
25 use mom_remapping,
only : initialize_remapping, remapping_core_h, end_remapping
31 implicit none ;
private
33 #include <MOM_memory.h>
35 public open_boundary_apply_normal_flow
36 public open_boundary_config
37 public open_boundary_init
38 public open_boundary_query
39 public open_boundary_end
40 public open_boundary_impose_normal_slope
41 public open_boundary_impose_land_mask
42 public radiation_open_bdry_conds
43 public set_tracer_data
44 public update_obc_segment_data
45 public open_boundary_test_extern_uv
46 public open_boundary_test_extern_h
47 public open_boundary_zero_normal_flow
48 public register_obc, obc_registry_init
49 public register_file_obc, file_obc_end
50 public segment_tracer_registry_init
51 public segment_tracer_registry_end
52 public register_segment_tracer
53 public register_temp_salt_segments
54 public fill_temp_salt_segments
55 public open_boundary_register_restarts
57 integer,
parameter,
public :: obc_none = 0
58 integer,
parameter,
public :: obc_simple = 1
59 integer,
parameter,
public :: obc_wall = 2
60 integer,
parameter,
public :: obc_flather = 3
61 integer,
parameter,
public :: obc_radiation = 4
62 integer,
parameter,
public :: obc_direction_n = 100
63 integer,
parameter,
public :: obc_direction_s = 200
64 integer,
parameter,
public :: obc_direction_e = 300
65 integer,
parameter,
public :: obc_direction_w = 400
66 integer,
parameter :: max_obc_fields = 100
72 character(len=8) :: name
73 real,
pointer,
dimension(:,:,:) :: buffer_src=>null()
76 real,
dimension(:,:,:),
pointer :: dz_src=>null()
77 real,
dimension(:,:,:),
pointer :: buffer_dst=>null()
78 real,
dimension(:,:),
pointer :: bt_vel=>null()
89 real,
dimension(:,:,:),
pointer :: t => null()
90 real :: obc_inflow_conc = 0.0
91 character(len=32) :: name
93 real,
dimension(:,:,:),
pointer :: tres => null()
94 logical :: is_initialized
101 logical :: locked = .false.
111 logical :: radiation_tan
113 logical :: radiation_grad
116 logical :: oblique_tan
118 logical :: oblique_grad
121 logical :: nudged_tan
122 logical :: nudged_grad
124 logical :: specified_tan
127 logical :: values_needed
132 integer :: num_fields
133 character(len=32),
pointer,
dimension(:) :: field_names=>null()
138 real :: velocity_nudging_timescale_in
139 real :: velocity_nudging_timescale_out
141 logical :: temp_segment_data_exists
142 logical :: salt_segment_data_exists
143 real,
pointer,
dimension(:,:) :: cg=>null()
145 real,
pointer,
dimension(:,:) :: htot=>null()
146 real,
pointer,
dimension(:,:,:) :: h=>null()
147 real,
pointer,
dimension(:,:,:) :: normal_vel=>null()
149 real,
pointer,
dimension(:,:,:) :: tangential_vel=>null()
151 real,
pointer,
dimension(:,:,:) :: tangential_grad=>null()
153 real,
pointer,
dimension(:,:,:) :: normal_trans=>null()
155 real,
pointer,
dimension(:,:) :: normal_vel_bt=>null()
157 real,
pointer,
dimension(:,:) :: eta=>null()
158 real,
pointer,
dimension(:,:,:) :: grad_normal=>null()
160 real,
pointer,
dimension(:,:,:) :: grad_tan=>null()
162 real,
pointer,
dimension(:,:,:) :: grad_gradient=>null()
164 real,
pointer,
dimension(:,:,:) :: rx_normal=>null()
166 real,
pointer,
dimension(:,:,:) :: ry_normal=>null()
168 real,
pointer,
dimension(:,:,:) :: cff_normal=>null()
170 real,
pointer,
dimension(:,:,:) :: nudged_normal_vel=>null()
172 real,
pointer,
dimension(:,:,:) :: nudged_tangential_vel=>null()
174 real,
pointer,
dimension(:,:,:) :: nudged_tangential_grad=>null()
177 type(hor_index_type) :: hi
178 real :: tr_invlscale3_out
179 real :: tr_invlscale3_in
187 integer :: number_of_segments = 0
189 logical :: open_u_bcs_exist_globally = .false.
191 logical :: open_v_bcs_exist_globally = .false.
193 logical :: flather_u_bcs_exist_globally = .false.
195 logical :: flather_v_bcs_exist_globally = .false.
197 logical :: oblique_bcs_exist_globally = .false.
199 logical :: nudged_u_bcs_exist_globally = .false.
201 logical :: nudged_v_bcs_exist_globally = .false.
203 logical :: specified_u_bcs_exist_globally = .false.
205 logical :: specified_v_bcs_exist_globally = .false.
207 logical :: radiation_bcs_exist_globally = .false.
208 logical :: user_bcs_set_globally = .false.
210 logical :: update_obc = .false.
211 logical :: needs_io_for_data = .false.
212 logical :: zero_vorticity = .false.
213 logical :: freeslip_vorticity = .false.
215 logical :: computed_vorticity = .false.
217 logical :: specified_vorticity = .false.
219 logical :: zero_strain = .false.
220 logical :: freeslip_strain = .false.
222 logical :: computed_strain = .false.
224 logical :: specified_strain = .false.
226 logical :: zero_biharmonic = .false.
228 logical :: brushcutter_mode = .false.
234 integer,
pointer,
dimension(:,:) :: &
235 segnum_u => null(), &
249 real,
pointer,
dimension(:,:,:) :: rx_normal => null()
250 real,
pointer,
dimension(:,:,:) :: ry_normal => null()
251 real,
pointer,
dimension(:,:,:) :: cff_normal => null()
263 real :: tide_flow = 3.0e6
268 character(len=32) :: name
275 logical :: locked = .false.
279 integer :: id_clock_pass
281 character(len=40) :: mdl =
"MOM_open_boundary"
283 #include "version_variable.h"
294 subroutine open_boundary_config(G, US, param_file, OBC)
301 logical :: debug_obc, debug, mask_outside, reentrant_x, reentrant_y
302 character(len=15) :: segment_param_str
303 character(len=100) :: segment_str
304 character(len=200) :: config1
305 real :: lscale_in, lscale_out
309 "Controls where open boundaries are located, what kind of boundary condition "//&
310 "to impose, and what data to apply, if any.")
311 call get_param(param_file, mdl,
"OBC_NUMBER_OF_SEGMENTS", obc%number_of_segments, &
312 "The number of open boundary segments.", &
314 call get_param(param_file, mdl,
"G_EARTH", obc%g_Earth, &
315 "The gravitational acceleration of the Earth.", &
316 units=
"m s-2", default = 9.80)
317 call get_param(param_file, mdl,
"OBC_USER_CONFIG", config1, &
318 "A string that sets how the open boundary conditions are "//&
319 " configured: \n", default=
"none", do_not_log=.true.)
320 call get_param(param_file, mdl,
"NK", obc%ke, &
321 "The number of model layers", default=0, do_not_log=.true.)
323 if (config1 /=
"none") obc%user_BCs_set_globally = .true.
325 if (obc%number_of_segments > 0)
then
326 call get_param(param_file, mdl,
"OBC_ZERO_VORTICITY", obc%zero_vorticity, &
327 "If true, sets relative vorticity to zero on open boundaries.", &
329 call get_param(param_file, mdl,
"OBC_FREESLIP_VORTICITY", obc%freeslip_vorticity, &
330 "If true, sets the normal gradient of tangential velocity to "//&
331 "zero in the relative vorticity on open boundaries. This cannot "//&
332 "be true if another OBC_XXX_VORTICITY option is True.", default=.true.)
333 call get_param(param_file, mdl,
"OBC_COMPUTED_VORTICITY", obc%computed_vorticity, &
334 "If true, uses the external values of tangential velocity "//&
335 "in the relative vorticity on open boundaries. This cannot "//&
336 "be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
337 call get_param(param_file, mdl,
"OBC_SPECIFIED_VORTICITY", obc%specified_vorticity, &
338 "If true, uses the external values of tangential velocity "//&
339 "in the relative vorticity on open boundaries. This cannot "//&
340 "be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
341 if ((obc%zero_vorticity .and. obc%freeslip_vorticity) .or. &
342 (obc%zero_vorticity .and. obc%computed_vorticity) .or. &
343 (obc%zero_vorticity .and. obc%specified_vorticity) .or. &
344 (obc%freeslip_vorticity .and. obc%computed_vorticity) .or. &
345 (obc%freeslip_vorticity .and. obc%specified_vorticity) .or. &
346 (obc%computed_vorticity .and. obc%specified_vorticity)) &
347 call mom_error(fatal,
"MOM_open_boundary.F90, open_boundary_config:\n"//&
348 "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"//&
349 "and OBC_IMPORTED_VORTICITY can be True at once.")
350 call get_param(param_file, mdl,
"OBC_ZERO_STRAIN", obc%zero_strain, &
351 "If true, sets the strain used in the stress tensor to zero on open boundaries.", &
353 call get_param(param_file, mdl,
"OBC_FREESLIP_STRAIN", obc%freeslip_strain, &
354 "If true, sets the normal gradient of tangential velocity to "//&
355 "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
356 "be true if another OBC_XXX_STRAIN option is True.", default=.true.)
357 call get_param(param_file, mdl,
"OBC_COMPUTED_STRAIN", obc%computed_strain, &
358 "If true, sets the normal gradient of tangential velocity to "//&
359 "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
360 "be true if another OBC_XXX_STRAIN option is True.", default=.false.)
361 call get_param(param_file, mdl,
"OBC_SPECIFIED_STRAIN", obc%specified_strain, &
362 "If true, sets the normal gradient of tangential velocity to "//&
363 "zero in the strain use in the stress tensor on open boundaries. This cannot "//&
364 "be true if another OBC_XXX_STRAIN option is True.", default=.false.)
365 if ((obc%zero_strain .and. obc%freeslip_strain) .or. &
366 (obc%zero_strain .and. obc%computed_strain) .or. &
367 (obc%zero_strain .and. obc%specified_strain) .or. &
368 (obc%freeslip_strain .and. obc%computed_strain) .or. &
369 (obc%freeslip_strain .and. obc%specified_strain) .or. &
370 (obc%computed_strain .and. obc%specified_strain)) &
371 call mom_error(fatal,
"MOM_open_boundary.F90, open_boundary_config: \n"//&
372 "Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN, OBC_COMPUTED_STRAIN \n"//&
373 "and OBC_IMPORTED_STRAIN can be True at once.")
374 call get_param(param_file, mdl,
"OBC_ZERO_BIHARMONIC", obc%zero_biharmonic, &
375 "If true, zeros the Laplacian of flow on open boundaries in the biharmonic "//&
376 "viscosity term.", default=.false.)
377 call get_param(param_file, mdl,
"MASK_OUTSIDE_OBCS", mask_outside, &
378 "If true, set the areas outside open boundaries to be land.", &
381 call get_param(param_file, mdl,
"DEBUG", debug, default=.false.)
382 call get_param(param_file, mdl,
"DEBUG_OBC", debug_obc, default=.false.)
383 if (debug_obc .or. debug) &
384 call log_param(param_file, mdl,
"DEBUG_OBC", debug_obc, &
385 "If true, do additional calls to help debug the performance "//&
386 "of the open boundary condition code.", default=.false., &
387 debuggingparam=.true.)
389 call get_param(param_file, mdl,
"OBC_SILLY_THICK", obc%silly_h, &
390 "A silly value of thicknesses used outside of open boundary "//&
391 "conditions for debugging.", units=
"m", default=0.0, &
392 do_not_log=.not.debug_obc, debuggingparam=.true.)
393 call get_param(param_file, mdl,
"OBC_SILLY_VEL", obc%silly_u, &
394 "A silly value of velocities used outside of open boundary "//&
395 "conditions for debugging.", units=
"m/s", default=0.0, &
396 do_not_log=.not.debug_obc, debuggingparam=.true.)
397 reentrant_x = .false.
398 call get_param(param_file, mdl,
"REENTRANT_X", reentrant_x, default=.true.)
399 reentrant_y = .false.
400 call get_param(param_file, mdl,
"REENTRANT_Y", reentrant_y, default=.false.)
404 allocate(obc%segment(0:obc%number_of_segments))
405 do l=0,obc%number_of_segments
406 obc%segment(l)%Flather = .false.
407 obc%segment(l)%radiation = .false.
408 obc%segment(l)%radiation_tan = .false.
409 obc%segment(l)%radiation_grad = .false.
410 obc%segment(l)%oblique = .false.
411 obc%segment(l)%oblique_tan = .false.
412 obc%segment(l)%oblique_grad = .false.
413 obc%segment(l)%nudged = .false.
414 obc%segment(l)%nudged_tan = .false.
415 obc%segment(l)%nudged_grad = .false.
416 obc%segment(l)%specified = .false.
417 obc%segment(l)%specified_tan = .false.
418 obc%segment(l)%open = .false.
419 obc%segment(l)%gradient = .false.
420 obc%segment(l)%values_needed = .false.
421 obc%segment(l)%direction = obc_none
422 obc%segment(l)%is_N_or_S = .false.
423 obc%segment(l)%is_E_or_W = .false.
424 obc%segment(l)%Velocity_nudging_timescale_in = 0.0
425 obc%segment(l)%Velocity_nudging_timescale_out = 0.0
426 obc%segment(l)%num_fields = 0.0
428 allocate(obc%segnum_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; obc%segnum_u(:,:) = obc_none
429 allocate(obc%segnum_v(g%isd:g%ied,g%JsdB:g%JedB)) ; obc%segnum_v(:,:) = obc_none
431 do l = 1, obc%number_of_segments
432 write(segment_param_str(1:15),
"('OBC_SEGMENT_',i3.3)") l
433 call get_param(param_file, mdl, segment_param_str, segment_str, &
434 "Documentation needs to be dynamic?????", &
435 fail_if_missing=.true.)
436 segment_str = remove_spaces(segment_str)
437 if (segment_str(1:2) ==
'I=')
then
438 call setup_u_point_obc(obc, g, segment_str, l, param_file, reentrant_y)
439 elseif (segment_str(1:2) ==
'J=')
then
440 call setup_v_point_obc(obc, g, segment_str, l, param_file, reentrant_x)
442 call mom_error(fatal,
"MOM_open_boundary.F90, open_boundary_config: "//&
443 "Unable to interpret "//segment_param_str//
" = "//trim(segment_str))
448 call initialize_segment_data(g, obc, param_file)
450 if (open_boundary_query(obc, apply_open_obc=.true.))
then
451 call get_param(param_file, mdl,
"OBC_RADIATION_MAX", obc%rx_max, &
452 "The maximum magnitude of the baroclinic radiation "//&
453 "velocity (or speed of characteristics). This is only "//&
454 "used if one of the open boundary segments is using Orlanski.", &
455 units=
"m s-1", default=10.0)
456 call get_param(param_file, mdl,
"OBC_RAD_VEL_WT", obc%gamma_uv, &
457 "The relative weighting for the baroclinic radiation "//&
458 "velocities (or speed of characteristics) at the new "//&
459 "time level (1) or the running mean (0) for velocities. "//&
460 "Valid values range from 0 to 1. This is only used if "//&
461 "one of the open boundary segments is using Orlanski.", &
462 units=
"nondim", default=0.3)
467 if (open_boundary_query(obc, apply_open_obc=.true.))
then
468 call get_param(param_file, mdl,
"OBC_TRACER_RESERVOIR_LENGTH_SCALE_OUT ", lscale_out, &
469 "An effective length scale for restoring the tracer concentration "//&
470 "at the boundaries to externally imposed values when the flow "//&
471 "is exiting the domain.", units=
"m", default=0.0)
473 call get_param(param_file, mdl,
"OBC_TRACER_RESERVOIR_LENGTH_SCALE_IN ", lscale_in, &
474 "An effective length scale for restoring the tracer concentration "//&
475 "at the boundaries to values from the interior when the flow "//&
476 "is entering the domain.", units=
"m", default=0.0)
479 if (mask_outside)
call mask_outside_obcs(g, us, param_file, obc)
484 do l = 1, obc%number_of_segments
485 obc%segment(l)%Tr_InvLscale3_in=0.0
486 if (lscale_in>0.) obc%segment(l)%Tr_InvLscale3_in = 1.0/(lscale_in*lscale_in*lscale_in)
487 obc%segment(l)%Tr_InvLscale3_out=0.0
488 if (lscale_out>0.) obc%segment(l)%Tr_InvLscale3_out = 1.0/(lscale_out*lscale_out*lscale_out)
494 if ((obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally) .and. &
495 .not.g%symmetric )
call mom_error(fatal, &
496 "MOM_open_boundary, open_boundary_config: "//&
497 "Symmetric memory must be used when using Flather OBCs.")
499 if (.not.(obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
500 obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally))
then
502 call open_boundary_dealloc(obc)
505 call time_interp_external_init()
508 end subroutine open_boundary_config
512 subroutine initialize_segment_data(G, OBC, PF)
513 use mpp_mod,
only : mpp_pe, mpp_set_current_pelist, mpp_get_current_pelist,mpp_npes
519 integer :: n,m,num_fields
520 character(len=256) :: segstr, filename
521 character(len=20) :: segnam, suffix
522 character(len=32) :: varnam, fieldname
525 character(len=32),
dimension(MAX_OBC_FIELDS) :: fields
526 character(len=128) :: inputdir
528 character(len=32) :: remappingScheme
529 logical :: check_reconstruction, check_remapping, force_bounds_in_subcell
530 integer,
dimension(4) :: siz,siz2
531 integer :: is, ie, js, je
532 integer :: isd, ied, jsd, jed
533 integer :: IsdB, IedB, JsdB, JedB
534 integer,
dimension(:),
allocatable :: saved_pelist
535 integer :: current_pe
536 integer,
dimension(1) :: single_pelist
539 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
546 call get_param(pf, mdl,
"INPUTDIR", inputdir, default=
".")
547 inputdir = slasher(inputdir)
549 call get_param(pf, mdl,
"REMAPPING_SCHEME", remappingscheme, &
550 "This sets the reconstruction scheme used "//&
551 "for vertical remapping for all variables. "//&
552 "It can be one of the following schemes: \n"//&
553 trim(remappingschemesdoc), default=remappingdefaultscheme,do_not_log=.true.)
554 call get_param(pf, mdl,
"FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
555 "If true, cell-by-cell reconstructions are checked for "//&
556 "consistency and if non-monotonicity or an inconsistency is "//&
557 "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
558 call get_param(pf, mdl,
"FATAL_CHECK_REMAPPING", check_remapping, &
559 "If true, the results of remapping are checked for "//&
560 "conservation and new extrema and if an inconsistency is "//&
561 "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
562 call get_param(pf, mdl,
"REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
563 "If true, the values on the intermediate grid used for remapping "//&
564 "are forced to be bounded, which might not be the case due to "//&
565 "round off.", default=.false.,do_not_log=.true.)
566 call get_param(pf, mdl,
"BRUSHCUTTER_MODE", obc%brushcutter_mode, &
567 "If true, read external OBC data on the supergrid.", &
570 allocate(obc%remap_CS)
571 call initialize_remapping(obc%remap_CS, remappingscheme, boundary_extrapolation = .false., &
572 check_reconstruction=check_reconstruction, &
573 check_remapping=check_remapping, force_bounds_in_subcell=force_bounds_in_subcell)
575 if (obc%user_BCs_set_globally)
return
578 do n=1, obc%number_of_segments
579 segment => obc%segment(n)
580 write(segnam,
"('OBC_SEGMENT_',i3.3,'_DATA')") n
581 call get_param(pf, mdl, segnam, segstr,
'OBC segment docs')
586 allocate(saved_pelist(0:mpp_npes()-1))
587 call mpp_get_current_pelist(saved_pelist)
588 current_pe = mpp_pe()
589 single_pelist(1) = current_pe
590 call mpp_set_current_pelist(single_pelist)
592 do n=1, obc%number_of_segments
593 segment => obc%segment(n)
595 write(segnam,
"('OBC_SEGMENT_',i3.3,'_DATA')") n
596 write(suffix,
"('_segment_',i3.3)") n
602 call parse_segment_data_str(trim(segstr), fields=fields, num_fields=num_fields)
603 if (num_fields == 0)
then
604 call mom_mesg(
'initialize_segment_data: num_fields = 0')
608 allocate(segment%field(num_fields))
610 if (segment%Flather)
then
611 if (num_fields < 3)
call mom_error(fatal, &
612 "MOM_open_boundary, initialize_segment_data: "//&
613 "Need at least three inputs for Flather")
615 segment%num_fields = num_fields
617 segment%temp_segment_data_exists=.false.
618 segment%salt_segment_data_exists=.false.
623 isd = segment%HI%isd ; ied = segment%HI%ied
624 jsd = segment%HI%jsd ; jed = segment%HI%jed
625 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
626 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
629 call parse_segment_data_str(trim(segstr), var=trim(fields(m)),
value=
value, filenam=filename, fieldnam=fieldname)
630 if (trim(filename) /=
'none')
then
631 obc%update_OBC = .true.
632 obc%needs_IO_for_data = .true.
633 segment%values_needed = .true.
634 segment%field(m)%name = trim(fields(m))
635 if (segment%field(m)%name ==
'TEMP') &
636 segment%temp_segment_data_exists=.true.
637 if (segment%field(m)%name ==
'SALT') &
638 segment%salt_segment_data_exists=.true.
639 filename = trim(inputdir)//trim(filename)
640 fieldname = trim(fieldname)//trim(suffix)
641 call field_size(filename,fieldname,siz,no_domain=.true.)
642 if (siz(4) == 1) segment%values_needed = .false.
643 if (segment%on_pe)
then
644 if (obc%brushcutter_mode .and. (modulo(siz(1),2) == 0 .or. modulo(siz(2),2) == 0))
then
645 call mom_error(fatal,
'segment data are not on the supergrid')
650 if (obc%brushcutter_mode)
then
658 if (obc%brushcutter_mode)
then
666 if (segment%is_E_or_W)
then
667 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
668 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
670 allocate(segment%field(m)%buffer_src(isdb:iedb,jsd:jed,siz2(3)))
673 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
674 allocate(segment%field(m)%buffer_src(isdb:iedb,jsdb:jedb,siz2(3)))
676 allocate(segment%field(m)%buffer_src(isd:ied,jsdb:jedb,siz2(3)))
679 segment%field(m)%buffer_src(:,:,:)=0.0
680 segment%field(m)%fid = init_external_field(trim(filename),&
681 trim(fieldname),ignore_axis_atts=.true.,threading=single_file)
683 fieldname =
'dz_'//trim(fieldname)
684 call field_size(filename,fieldname,siz,no_domain=.true.)
685 if (segment%is_E_or_W)
then
686 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
687 allocate(segment%field(m)%dz_src(isdb:iedb,jsdb:jedb,siz(3)))
689 allocate(segment%field(m)%dz_src(isdb:iedb,jsd:jed,siz(3)))
692 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
693 allocate(segment%field(m)%dz_src(isdb:iedb,jsdb:jedb,siz(3)))
695 allocate(segment%field(m)%dz_src(isd:ied,jsdb:jedb,siz(3)))
698 segment%field(m)%dz_src(:,:,:)=0.0
699 segment%field(m)%nk_src=siz(3)
700 segment%field(m)%fid_dz = init_external_field(trim(filename),trim(fieldname),&
701 ignore_axis_atts=.true.,threading=single_file)
703 segment%field(m)%nk_src=1
707 segment%field(m)%fid = -1
708 segment%field(m)%value =
value
713 call mpp_set_current_pelist(saved_pelist)
715 end subroutine initialize_segment_data
719 subroutine setup_segment_indices(G, seg, Is_obc, Ie_obc, Js_obc, Je_obc)
722 integer,
intent(in) :: Is_obc
723 integer,
intent(in) :: Ie_obc
724 integer,
intent(in) :: Js_obc
725 integer,
intent(in) :: Je_obc
727 integer :: Isg,Ieg,Jsg,Jeg
730 if (ie_obc<is_obc)
then
731 isg=ie_obc;ieg=is_obc
733 isg=is_obc;ieg=ie_obc
735 if (je_obc<js_obc)
then
736 jsg=je_obc;jeg=js_obc
738 jsg=js_obc;jeg=je_obc
742 seg%HI%IsgB = isg ; seg%HI%IegB = ieg
743 seg%HI%isg = isg+1 ; seg%HI%ieg = ieg
744 seg%HI%JsgB = jsg ; seg%HI%JegB = jeg
745 seg%HI%jsg = jsg+1 ; seg%HI%Jeg = jeg
748 isg = isg - g%idg_offset
749 jsg = jsg - g%jdg_offset
750 ieg = ieg - g%idg_offset
751 jeg = jeg - g%jdg_offset
755 seg%HI%IsdB = min( max(isg, g%HI%IsdB), g%HI%IedB)
756 seg%HI%IedB = min( max(ieg, g%HI%IsdB), g%HI%IedB)
757 seg%HI%isd = min( max(isg+1, g%HI%isd), g%HI%ied)
758 seg%HI%ied = min( max(ieg, g%HI%isd), g%HI%ied)
759 seg%HI%IscB = min( max(isg, g%HI%IscB), g%HI%IecB)
760 seg%HI%IecB = min( max(ieg, g%HI%IscB), g%HI%IecB)
761 seg%HI%isc = min( max(isg+1, g%HI%isc), g%HI%iec)
762 seg%HI%iec = min( max(ieg, g%HI%isc), g%HI%iec)
766 seg%HI%JsdB = min( max(jsg, g%HI%JsdB), g%HI%JedB)
767 seg%HI%JedB = min( max(jeg, g%HI%JsdB), g%HI%JedB)
768 seg%HI%jsd = min( max(jsg+1, g%HI%jsd), g%HI%jed)
769 seg%HI%jed = min( max(jeg, g%HI%jsd), g%HI%jed)
770 seg%HI%JscB = min( max(jsg, g%HI%JscB), g%HI%JecB)
771 seg%HI%JecB = min( max(jeg, g%HI%JscB), g%HI%JecB)
772 seg%HI%jsc = min( max(jsg+1, g%HI%jsc), g%HI%jec)
773 seg%HI%jec = min( max(jeg, g%HI%jsc), g%HI%jec)
775 end subroutine setup_segment_indices
778 subroutine setup_u_point_obc(OBC, G, segment_str, l_seg, PF, reentrant_y)
781 character(len=*),
intent(in) :: segment_str
782 integer,
intent(in) :: l_seg
784 logical,
intent(in) :: reentrant_y
786 integer :: I_obc, Js_obc, Je_obc
788 character(len=32) :: action_str(8)
789 character(len=128) :: segment_param_str
790 real,
allocatable,
dimension(:) :: tnudge
792 call parse_segment_str(g%ieg, g%jeg, segment_str, i_obc, js_obc, je_obc, action_str, reentrant_y)
794 call setup_segment_indices(g, obc%segment(l_seg),i_obc,i_obc,js_obc,je_obc)
796 i_obc = i_obc - g%idg_offset
797 js_obc = js_obc - g%jdg_offset
798 je_obc = je_obc - g%jdg_offset
800 if (je_obc>js_obc)
then
801 obc%segment(l_seg)%direction = obc_direction_e
802 elseif (je_obc<js_obc)
then
803 obc%segment(l_seg)%direction = obc_direction_w
804 j=js_obc;js_obc=je_obc;je_obc=j
807 obc%segment(l_seg)%on_pe = .false.
810 if (len_trim(action_str(a_loop)) == 0)
then
812 elseif (trim(action_str(a_loop)) ==
'FLATHER')
then
813 obc%segment(l_seg)%Flather = .true.
814 obc%segment(l_seg)%open = .true.
815 obc%Flather_u_BCs_exist_globally = .true.
816 obc%open_u_BCs_exist_globally = .true.
817 elseif (trim(action_str(a_loop)) ==
'ORLANSKI')
then
818 obc%segment(l_seg)%radiation = .true.
819 obc%segment(l_seg)%open = .true.
820 obc%open_u_BCs_exist_globally = .true.
821 obc%radiation_BCs_exist_globally = .true.
822 elseif (trim(action_str(a_loop)) ==
'ORLANSKI_TAN')
then
823 obc%segment(l_seg)%radiation = .true.
824 obc%segment(l_seg)%radiation_tan = .true.
825 obc%radiation_BCs_exist_globally = .true.
826 elseif (trim(action_str(a_loop)) ==
'ORLANSKI_GRAD')
then
827 obc%segment(l_seg)%radiation = .true.
828 obc%segment(l_seg)%radiation_grad = .true.
829 elseif (trim(action_str(a_loop)) ==
'OBLIQUE')
then
830 obc%segment(l_seg)%oblique = .true.
831 obc%segment(l_seg)%open = .true.
832 obc%oblique_BCs_exist_globally = .true.
833 obc%open_u_BCs_exist_globally = .true.
834 elseif (trim(action_str(a_loop)) ==
'OBLIQUE_TAN')
then
835 obc%segment(l_seg)%oblique = .true.
836 obc%segment(l_seg)%oblique_tan = .true.
837 obc%oblique_BCs_exist_globally = .true.
838 elseif (trim(action_str(a_loop)) ==
'OBLIQUE_GRAD')
then
839 obc%segment(l_seg)%oblique = .true.
840 obc%segment(l_seg)%oblique_grad = .true.
841 elseif (trim(action_str(a_loop)) ==
'NUDGED')
then
842 obc%segment(l_seg)%nudged = .true.
843 obc%nudged_u_BCs_exist_globally = .true.
844 elseif (trim(action_str(a_loop)) ==
'NUDGED_TAN')
then
845 obc%segment(l_seg)%nudged_tan = .true.
846 obc%nudged_u_BCs_exist_globally = .true.
847 elseif (trim(action_str(a_loop)) ==
'NUDGED_GRAD')
then
848 obc%segment(l_seg)%nudged_grad = .true.
849 elseif (trim(action_str(a_loop)) ==
'GRADIENT')
then
850 obc%segment(l_seg)%gradient = .true.
851 obc%segment(l_seg)%open = .true.
852 obc%open_u_BCs_exist_globally = .true.
853 elseif (trim(action_str(a_loop)) ==
'SIMPLE')
then
854 obc%segment(l_seg)%specified = .true.
855 obc%specified_u_BCs_exist_globally = .true.
856 elseif (trim(action_str(a_loop)) ==
'SIMPLE_TAN')
then
857 obc%segment(l_seg)%specified_tan = .true.
859 call mom_error(fatal,
"MOM_open_boundary.F90, setup_u_point_obc: "//&
860 "String '"//trim(action_str(a_loop))//
"' not understood.")
862 if (obc%segment(l_seg)%nudged .or. obc%segment(l_seg)%nudged_tan)
then
863 write(segment_param_str(1:43),
"('OBC_SEGMENT_',i3.3,'_VELOCITY_NUDGING_TIMESCALES')") l_seg
865 call get_param(pf, mdl, segment_param_str(1:43), tnudge, &
866 "Timescales in days for nudging along a segment, "//&
867 "for inflow, then outflow. Setting both to zero should "//&
868 "behave like SIMPLE obcs for the baroclinic velocities.", &
869 fail_if_missing=.true.,default=0.,units=
"days")
870 obc%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1)*86400.
871 obc%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2)*86400.
877 if (i_obc<=g%HI%IsdB+1 .or. i_obc>=g%HI%IedB-1)
return
878 if (je_obc<=g%HI%JsdB .or. js_obc>=g%HI%JedB)
return
880 obc%segment(l_seg)%on_pe = .true.
881 obc%segment(l_seg)%is_E_or_W = .true.
883 do j=g%HI%jsd, g%HI%jed
884 if (j>js_obc .and. j<=je_obc)
then
885 obc%segnum_u(i_obc,j) = l_seg
888 obc%segment(l_seg)%Is_obc = i_obc
889 obc%segment(l_seg)%Ie_obc = i_obc
890 obc%segment(l_seg)%Js_obc = js_obc
891 obc%segment(l_seg)%Je_obc = je_obc
892 call allocate_obc_segment_data(obc, obc%segment(l_seg))
894 if (obc%segment(l_seg)%oblique .and. obc%segment(l_seg)%radiation) &
895 call mom_error(fatal,
"MOM_open_boundary.F90, setup_u_point_obc: \n"//&
896 "Orlanski and Oblique OBC options cannot be used together on one segment.")
898 end subroutine setup_u_point_obc
901 subroutine setup_v_point_obc(OBC, G, segment_str, l_seg, PF, reentrant_x)
904 character(len=*),
intent(in) :: segment_str
905 integer,
intent(in) :: l_seg
907 logical,
intent(in) :: reentrant_x
909 integer :: J_obc, Is_obc, Ie_obc
911 character(len=32) :: action_str(8)
912 character(len=128) :: segment_param_str
913 real,
allocatable,
dimension(:) :: tnudge
916 call parse_segment_str(g%ieg, g%jeg, segment_str, j_obc, is_obc, ie_obc, action_str, reentrant_x)
918 call setup_segment_indices(g, obc%segment(l_seg),is_obc,ie_obc,j_obc,j_obc)
920 j_obc = j_obc - g%jdg_offset
921 is_obc = is_obc - g%idg_offset
922 ie_obc = ie_obc - g%idg_offset
924 if (ie_obc>is_obc)
then
925 obc%segment(l_seg)%direction = obc_direction_s
926 elseif (ie_obc<is_obc)
then
927 obc%segment(l_seg)%direction = obc_direction_n
928 i=is_obc;is_obc=ie_obc;ie_obc=i
931 obc%segment(l_seg)%on_pe = .false.
934 if (len_trim(action_str(a_loop)) == 0)
then
936 elseif (trim(action_str(a_loop)) ==
'FLATHER')
then
937 obc%segment(l_seg)%Flather = .true.
938 obc%segment(l_seg)%open = .true.
939 obc%Flather_v_BCs_exist_globally = .true.
940 obc%open_v_BCs_exist_globally = .true.
941 elseif (trim(action_str(a_loop)) ==
'ORLANSKI')
then
942 obc%segment(l_seg)%radiation = .true.
943 obc%segment(l_seg)%open = .true.
944 obc%open_v_BCs_exist_globally = .true.
945 obc%radiation_BCs_exist_globally = .true.
946 elseif (trim(action_str(a_loop)) ==
'ORLANSKI_TAN')
then
947 obc%segment(l_seg)%radiation = .true.
948 obc%segment(l_seg)%radiation_tan = .true.
949 obc%radiation_BCs_exist_globally = .true.
950 elseif (trim(action_str(a_loop)) ==
'ORLANSKI_GRAD')
then
951 obc%segment(l_seg)%radiation = .true.
952 obc%segment(l_seg)%radiation_grad = .true.
953 elseif (trim(action_str(a_loop)) ==
'OBLIQUE')
then
954 obc%segment(l_seg)%oblique = .true.
955 obc%segment(l_seg)%open = .true.
956 obc%oblique_BCs_exist_globally = .true.
957 obc%open_v_BCs_exist_globally = .true.
958 elseif (trim(action_str(a_loop)) ==
'OBLIQUE_TAN')
then
959 obc%segment(l_seg)%oblique = .true.
960 obc%segment(l_seg)%oblique_tan = .true.
961 obc%oblique_BCs_exist_globally = .true.
962 elseif (trim(action_str(a_loop)) ==
'OBLIQUE_GRAD')
then
963 obc%segment(l_seg)%oblique = .true.
964 obc%segment(l_seg)%oblique_grad = .true.
965 elseif (trim(action_str(a_loop)) ==
'NUDGED')
then
966 obc%segment(l_seg)%nudged = .true.
967 obc%nudged_v_BCs_exist_globally = .true.
968 elseif (trim(action_str(a_loop)) ==
'NUDGED_TAN')
then
969 obc%segment(l_seg)%nudged_tan = .true.
970 obc%nudged_v_BCs_exist_globally = .true.
971 elseif (trim(action_str(a_loop)) ==
'NUDGED_GRAD')
then
972 obc%segment(l_seg)%nudged_grad = .true.
973 elseif (trim(action_str(a_loop)) ==
'GRADIENT')
then
974 obc%segment(l_seg)%gradient = .true.
975 obc%segment(l_seg)%open = .true.
976 obc%open_v_BCs_exist_globally = .true.
977 elseif (trim(action_str(a_loop)) ==
'SIMPLE')
then
978 obc%segment(l_seg)%specified = .true.
979 obc%specified_v_BCs_exist_globally = .true.
980 elseif (trim(action_str(a_loop)) ==
'SIMPLE_TAN')
then
981 obc%segment(l_seg)%specified_tan = .true.
983 call mom_error(fatal,
"MOM_open_boundary.F90, setup_v_point_obc: "//&
984 "String '"//trim(action_str(a_loop))//
"' not understood.")
986 if (obc%segment(l_seg)%nudged .or. obc%segment(l_seg)%nudged_tan)
then
987 write(segment_param_str(1:43),
"('OBC_SEGMENT_',i3.3,'_VELOCITY_NUDGING_TIMESCALES')") l_seg
989 call get_param(pf, mdl, segment_param_str(1:43), tnudge, &
990 "Timescales in days for nudging along a segment, "//&
991 "for inflow, then outflow. Setting both to zero should "//&
992 "behave like SIMPLE obcs for the baroclinic velocities.", &
993 fail_if_missing=.true.,default=0.,units=
"days")
994 obc%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1)*86400.
995 obc%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2)*86400.
1001 if (j_obc<=g%HI%JsdB+1 .or. j_obc>=g%HI%JedB-1)
return
1002 if (ie_obc<=g%HI%IsdB .or. is_obc>=g%HI%IedB)
return
1004 obc%segment(l_seg)%on_pe = .true.
1005 obc%segment(l_seg)%is_N_or_S = .true.
1007 do i=g%HI%isd, g%HI%ied
1008 if (i>is_obc .and. i<=ie_obc)
then
1009 obc%segnum_v(i,j_obc) = l_seg
1012 obc%segment(l_seg)%Is_obc = is_obc
1013 obc%segment(l_seg)%Ie_obc = ie_obc
1014 obc%segment(l_seg)%Js_obc = j_obc
1015 obc%segment(l_seg)%Je_obc = j_obc
1016 call allocate_obc_segment_data(obc, obc%segment(l_seg))
1018 if (obc%segment(l_seg)%oblique .and. obc%segment(l_seg)%radiation) &
1019 call mom_error(fatal,
"MOM_open_boundary.F90, setup_v_point_obc: \n"//&
1020 "Orlanski and Oblique OBC options cannot be used together on one segment.")
1022 end subroutine setup_v_point_obc
1025 subroutine parse_segment_str(ni_global, nj_global, segment_str, l, m, n, action_str, reentrant)
1026 integer,
intent(in) :: ni_global
1027 integer,
intent(in) :: nj_global
1028 character(len=*),
intent(in) :: segment_str
1029 integer,
intent(out) :: l
1030 integer,
intent(out) :: m
1031 integer,
intent(out) :: n
1032 character(len=*),
intent(out) :: action_str(:)
1033 logical,
intent(in) :: reentrant
1035 character(len=24) :: word1, word2, m_word, n_word
1040 integer,
parameter :: halo = 10
1043 word1 = extract_word(segment_str,
',',1)
1044 word2 = extract_word(segment_str,
',',2)
1045 if (word1(1:2)==
'I=')
then
1048 if (.not. (word2(1:2)==
'J='))
call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1049 "Second word of string '"//trim(segment_str)//
"' must start with 'J='.")
1050 elseif (word1(1:2)==
'J=')
then
1053 if (.not. (word2(1:2)==
'I='))
call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1054 "Second word of string '"//trim(segment_str)//
"' must start with 'I='.")
1056 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str"//&
1057 "String '"//segment_str//
"' must start with 'I=' or 'J='.")
1061 l = interpret_int_expr( word1(3:24), l_max )
1062 if (l<0 .or. l>l_max)
then
1063 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1064 "First value from string '"//trim(segment_str)//
"' is outside of the physical domain.")
1068 m_word = extract_word(word2(3:24),
':',1)
1069 m = interpret_int_expr( m_word, mn_max )
1071 if (m<-halo .or. m>mn_max+halo)
then
1072 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1073 "Beginning of range in string '"//trim(segment_str)//
"' is outside of the physical domain.")
1076 if (m<-1 .or. m>mn_max+1)
then
1077 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1078 "Beginning of range in string '"//trim(segment_str)//
"' is outside of the physical domain.")
1083 n_word = extract_word(word2(3:24),
':',2)
1084 n = interpret_int_expr( n_word, mn_max )
1086 if (n<-halo .or. n>mn_max+halo)
then
1087 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1088 "End of range in string '"//trim(segment_str)//
"' is outside of the physical domain.")
1091 if (n<-1 .or. n>mn_max+1)
then
1092 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1093 "End of range in string '"//trim(segment_str)//
"' is outside of the physical domain.")
1097 if (abs(n-m)==0)
then
1098 call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str: "//&
1099 "Range in string '"//trim(segment_str)//
"' must span one cell.")
1103 do j = 1,
size(action_str)
1104 action_str(j) = extract_word(segment_str,
',',2+j)
1110 integer function interpret_int_expr(string, imax)
1111 character(len=*),
intent(in) :: string
1112 integer,
intent(in) :: imax
1116 slen = len_trim(string)
1117 if (slen==0)
call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str"//&
1118 "Parsed string was empty!")
1119 if (len_trim(string)==1 .and. string(1:1)==
'N')
then
1120 interpret_int_expr = imax
1121 elseif (string(1:1)==
'N')
then
1122 if (string(2:2)==
'+')
then
1123 read(string(3:slen),*,err=911) interpret_int_expr
1124 interpret_int_expr = imax + interpret_int_expr
1125 elseif (string(2:2)==
'-')
then
1126 read(string(3:slen),*,err=911) interpret_int_expr
1127 interpret_int_expr = imax - interpret_int_expr
1130 read(string(1:slen),*,err=911) interpret_int_expr
1133 911
call mom_error(fatal,
"MOM_open_boundary.F90, parse_segment_str"//&
1134 "Problem reading value from string '"//trim(string)//
"'.")
1135 end function interpret_int_expr
1136 end subroutine parse_segment_str
1139 subroutine parse_segment_data_str(segment_str, var, value, filenam, fieldnam, fields, num_fields, debug )
1140 character(len=*),
intent(in) :: segment_str
1142 character(len=*),
optional,
intent(in) :: var
1143 character(len=*),
optional,
intent(out) :: filenam
1144 character(len=*),
optional,
intent(out) :: fieldnam
1146 real,
optional,
intent(out) ::
value
1147 character(len=*),
dimension(MAX_OBC_FIELDS), &
1148 optional,
intent(out) :: fields
1149 integer,
optional,
intent(out) :: num_fields
1150 logical,
optional,
intent(in) :: debug
1152 character(len=128) :: word1, word2, word3, method
1153 integer :: lword, nfields, n, m, orient
1154 logical :: continue,dbg
1155 character(len=32),
dimension(MAX_OBC_FIELDS) :: flds
1160 if (
PRESENT(debug)) dbg=debug
1163 word1 = extract_word(segment_str,
',',nfields+1)
1164 if (trim(word1) ==
'')
exit
1166 word2 = extract_word(word1,
'=',1)
1167 flds(nfields) = trim(word2)
1170 if (
PRESENT(fields))
then
1176 if (
PRESENT(num_fields))
then
1182 if (
PRESENT(var))
then
1184 if (trim(var)==trim(flds(n)))
then
1194 word3 = extract_word(segment_str,
',',m)
1195 word1 = extract_word(word3,
':',1)
1197 word2 = extract_word(word1,
'=',1)
1198 if (trim(word2) == trim(var))
then
1199 method=trim(extract_word(word1,
'=',2))
1200 lword=len_trim(method)
1201 if (method(lword-3:lword) ==
'file')
then
1203 word1 = extract_word(word3,
':',2)
1204 filenam = extract_word(word1,
'(',1)
1205 fieldnam = extract_word(word1,
'(',2)
1206 lword=len_trim(fieldnam)
1207 fieldnam = fieldnam(1:lword-1)
1209 elseif (method(lword-4:lword) ==
'value')
then
1212 word1 = extract_word(word3,
':',2)
1213 lword=len_trim(word1)
1214 read(word1(1:lword),*,
end=986,err=987) value
1220 986
call mom_error(fatal,
'End of record while parsing segment data specification! '//trim(segment_str))
1221 987
call mom_error(fatal,
'Error while parsing segment data specification! '//trim(segment_str))
1223 end subroutine parse_segment_data_str
1227 subroutine parse_segment_param_real(segment_str, var, param_value, debug )
1228 character(len=*),
intent(in) :: segment_str
1230 character(len=*),
intent(in) :: var
1231 real,
intent(out) :: param_value
1232 logical,
optional,
intent(in) :: debug
1234 character(len=128) :: word1, word2, word3, method
1235 integer :: lword, nfields, n, m, orient
1236 logical :: continue,dbg
1237 character(len=32),
dimension(MAX_OBC_FIELDS) :: flds
1242 if (
PRESENT(debug)) dbg=debug
1245 word1 = extract_word(segment_str,
',',nfields+1)
1246 if (trim(word1) ==
'')
exit
1248 word2 = extract_word(word1,
'=',1)
1249 flds(nfields) = trim(word2)
1266 if (trim(var)==trim(flds(n)))
then
1276 word3 = extract_word(segment_str,
',',m)
1279 word2 = extract_word(word1,
'=',1)
1280 if (trim(word2) == trim(var))
then
1281 method=trim(extract_word(word1,
'=',2))
1282 lword=len_trim(method)
1283 read(method(1:lword),*,err=987) param_value
1303 986
call mom_error(fatal,
'End of record while parsing segment data specification! '//trim(segment_str))
1304 987
call mom_error(fatal,
'Error while parsing segment parameter specification! '//trim(segment_str))
1306 end subroutine parse_segment_param_real
1309 subroutine open_boundary_init(G, param_file, OBC)
1310 type(ocean_grid_type),
intent(in) :: g
1311 type(param_file_type),
intent(in) :: param_file
1315 if (.not.
associated(obc))
return
1317 id_clock_pass = cpu_clock_id(
'(Ocean OBC halo updates)', grain=clock_routine)
1319 end subroutine open_boundary_init
1321 logical function open_boundary_query(OBC, apply_open_OBC, apply_specified_OBC, apply_Flather_OBC, &
1322 apply_nudged_OBC, needs_ext_seg_data)
1324 logical,
optional,
intent(in) :: apply_open_obc
1325 logical,
optional,
intent(in) :: apply_specified_obc
1326 logical,
optional,
intent(in) :: apply_flather_obc
1327 logical,
optional,
intent(in) :: apply_nudged_obc
1328 logical,
optional,
intent(in) :: needs_ext_seg_data
1329 open_boundary_query = .false.
1330 if (.not.
associated(obc))
return
1331 if (
present(apply_open_obc)) open_boundary_query = obc%open_u_BCs_exist_globally .or. &
1332 obc%open_v_BCs_exist_globally
1333 if (
present(apply_specified_obc)) open_boundary_query = obc%specified_u_BCs_exist_globally .or. &
1334 obc%specified_v_BCs_exist_globally
1335 if (
present(apply_flather_obc)) open_boundary_query = obc%Flather_u_BCs_exist_globally .or. &
1336 obc%Flather_v_BCs_exist_globally
1337 if (
present(apply_nudged_obc)) open_boundary_query = obc%nudged_u_BCs_exist_globally .or. &
1338 obc%nudged_v_BCs_exist_globally
1339 if (
present(needs_ext_seg_data)) open_boundary_query = obc%needs_IO_for_data
1341 end function open_boundary_query
1344 subroutine open_boundary_dealloc(OBC)
1349 if (.not.
associated(obc))
return
1351 do n=1, obc%number_of_segments
1352 segment => obc%segment(n)
1353 call deallocate_obc_segment_data(obc, segment)
1355 if (
associated(obc%segment))
deallocate(obc%segment)
1356 if (
associated(obc%segnum_u))
deallocate(obc%segnum_u)
1357 if (
associated(obc%segnum_v))
deallocate(obc%segnum_v)
1359 end subroutine open_boundary_dealloc
1362 subroutine open_boundary_end(OBC)
1364 call open_boundary_dealloc(obc)
1365 end subroutine open_boundary_end
1368 subroutine open_boundary_impose_normal_slope(OBC, G, depth)
1370 type(dyn_horgrid_type),
intent(in) :: g
1371 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: depth
1376 if (.not.
associated(obc))
return
1378 if (.not.(obc%specified_u_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally .or. &
1379 obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
1382 do n=1,obc%number_of_segments
1383 segment=>obc%segment(n)
1384 if (.not. segment%on_pe) cycle
1385 if (segment%direction == obc_direction_e)
then
1387 do j=segment%HI%jsd,segment%HI%jed
1388 depth(i+1,j) = depth(i,j)
1390 elseif (segment%direction == obc_direction_w)
then
1392 do j=segment%HI%jsd,segment%HI%jed
1393 depth(i,j) = depth(i+1,j)
1395 elseif (segment%direction == obc_direction_n)
then
1397 do i=segment%HI%isd,segment%HI%ied
1398 depth(i,j+1) = depth(i,j)
1400 elseif (segment%direction == obc_direction_s)
then
1402 do i=segment%HI%isd,segment%HI%ied
1403 depth(i,j) = depth(i,j+1)
1408 end subroutine open_boundary_impose_normal_slope
1413 subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv)
1415 type(dyn_horgrid_type),
intent(inout) :: g
1416 real,
dimension(SZIB_(G),SZJ_(G)),
intent(inout) :: areacu
1417 real,
dimension(SZI_(G),SZJB_(G)),
intent(inout) :: areacv
1421 logical :: any_u, any_v
1423 if (.not.
associated(obc))
return
1425 do n=1,obc%number_of_segments
1426 segment=>obc%segment(n)
1427 if (.not. segment%on_pe) cycle
1428 if (segment%is_E_or_W)
then
1432 do j=segment%HI%jsd,segment%HI%jed
1433 if (g%mask2dCu(i,j) == 0) obc%segnum_u(i,j) = obc_none
1434 if (segment%direction == obc_direction_w)
then
1437 g%mask2dT(i+1,j) = 0
1440 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1441 if (segment%direction == obc_direction_w)
then
1444 g%mask2dCv(i+1,j) = 0
1450 do i=segment%HI%isd,segment%HI%ied
1451 if (g%mask2dCv(i,j) == 0) obc%segnum_v(i,j) = obc_none
1452 if (segment%direction == obc_direction_s)
then
1455 g%mask2dT(i,j+1) = 0
1458 do i=segment%HI%IsdB+1,segment%HI%IedB-1
1459 if (segment%direction == obc_direction_s)
then
1462 g%mask2dCu(i,j+1) = 0
1468 do n=1,obc%number_of_segments
1469 segment=>obc%segment(n)
1470 if (.not. segment%on_pe .or. .not. segment%specified) cycle
1471 if (segment%is_E_or_W)
then
1474 do j=segment%HI%jsd,segment%HI%jed
1475 if (segment%direction == obc_direction_e)
then
1476 areacu(i,j) = g%areaT(i,j)
1478 areacu(i,j) = g%areaT(i+1,j)
1484 do i=segment%HI%isd,segment%HI%ied
1485 if (segment%direction == obc_direction_s)
then
1486 areacv(i,j) = g%areaT(i,j+1)
1488 areacu(i,j) = g%areaT(i,j)
1500 do n=1,obc%number_of_segments
1501 segment=>obc%segment(n)
1502 if (.not. segment%on_pe) cycle
1503 if (segment%is_E_or_W)
then
1505 do j=segment%HI%jsd,segment%HI%jed
1506 if (obc%segnum_u(i,j) /= obc_none) any_u = .true.
1510 do i=segment%HI%isd,segment%HI%ied
1511 if (obc%segnum_v(i,j) /= obc_none) any_v = .true.
1517 if (.not.(any_u .or. any_v)) obc%OBC_pe = .false.
1519 end subroutine open_boundary_impose_land_mask
1522 subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt)
1523 type(ocean_grid_type),
intent(inout) :: g
1525 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u_new
1528 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_old
1529 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v_new
1532 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_old
1533 real,
intent(in) :: dt
1535 real :: dhdt, dhdx, dhdy, gamma_u, gamma_v, gamma_2
1536 real :: cff, cx, cy, tau
1537 real :: rx_max, ry_max
1538 real :: rx_new, rx_avg
1539 real :: ry_new, ry_avg
1540 real :: cff_new, cff_avg
1541 real,
pointer,
dimension(:,:,:) :: rx_tangential=>null()
1542 real,
pointer,
dimension(:,:,:) :: ry_tangential=>null()
1543 real,
pointer,
dimension(:,:,:) :: cff_tangential=>null()
1544 real,
parameter :: eps = 1.0e-20
1546 integer :: i, j, k, is, ie, js, je, nz, n
1547 integer :: is_obc, ie_obc, js_obc, je_obc
1549 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1551 if (.not.
associated(obc))
return
1553 if (.not.(obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
1559 do n=1,obc%number_of_segments
1560 segment=>obc%segment(n)
1561 if (.not. segment%on_pe) cycle
1562 if (segment%is_E_or_W .and. segment%radiation)
then
1565 do j=segment%HI%jsd,segment%HI%jed
1566 segment%rx_normal(i,j,k) = obc%rx_normal(i,j,k)
1569 elseif (segment%is_N_or_S .and. segment%radiation)
then
1572 do i=segment%HI%isd,segment%HI%ied
1573 segment%ry_normal(i,j,k) = obc%ry_normal(i,j,k)
1577 if (segment%is_E_or_W .and. segment%oblique)
then
1580 do j=segment%HI%jsd,segment%HI%jed
1581 segment%rx_normal(i,j,k) = obc%rx_normal(i,j,k)
1582 segment%ry_normal(i,j,k) = obc%ry_normal(i,j,k)
1583 segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
1586 elseif (segment%is_N_or_S .and. segment%oblique)
then
1589 do i=segment%HI%isd,segment%HI%ied
1590 segment%rx_normal(i,j,k) = obc%rx_normal(i,j,k)
1591 segment%ry_normal(i,j,k) = obc%ry_normal(i,j,k)
1592 segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
1598 gamma_u = obc%gamma_uv ; gamma_v = obc%gamma_uv
1599 rx_max = obc%rx_max ; ry_max = obc%rx_max
1600 do n=1,obc%number_of_segments
1601 segment=>obc%segment(n)
1602 if (.not. segment%on_pe) cycle
1603 if (segment%oblique)
call gradient_at_q_points(g,segment,u_new,v_new)
1604 if (segment%direction == obc_direction_e)
then
1606 if (i<g%HI%IscB) cycle
1607 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
1608 if (segment%radiation)
then
1609 dhdt = u_old(i-1,j,k)-u_new(i-1,j,k)
1610 dhdx = u_new(i-1,j,k)-u_new(i-2,j,k)
1612 if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
1613 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
1614 segment%rx_normal(i,j,k) = rx_avg
1618 segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) / (1.0+rx_avg)
1621 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
1622 elseif (segment%oblique)
then
1623 dhdt = u_old(i-1,j,k)-u_new(i-1,j,k)
1624 dhdx = u_new(i-1,j,k)-u_new(i-2,j,k)
1625 if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0)
then
1626 dhdy = segment%grad_normal(j-1,1,k)
1627 elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0)
then
1630 dhdy = segment%grad_normal(j,1,k)
1632 if (dhdt*dhdx < 0.0) dhdt = 0.0
1634 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
1635 ry_new = min(cff,max(dhdt*dhdy,-cff))
1636 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
1637 ry_avg = (1.0-gamma_u)*segment%ry_normal(i,j,k) + gamma_u*ry_new
1638 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
1639 segment%rx_normal(i,j,k) = rx_avg
1640 segment%ry_normal(i,j,k) = ry_avg
1641 segment%cff_normal(i,j,k) = cff_avg
1642 segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) - &
1643 (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
1647 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
1648 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
1649 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
1650 elseif (segment%gradient)
then
1651 segment%normal_vel(i,j,k) = u_new(i-1,j,k)
1653 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then
1655 if (dhdt*dhdx <= 0.0)
then
1656 tau = segment%Velocity_nudging_timescale_in
1658 tau = segment%Velocity_nudging_timescale_out
1660 gamma_2 = dt / (tau + dt)
1661 segment%normal_vel(i,j,k) = (1 - gamma_2) * segment%normal_vel(i,j,k) + &
1662 gamma_2 * segment%nudged_normal_vel(i,j,k)
1665 if (segment%radiation_tan .or. segment%radiation_grad)
then
1667 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1669 rx_tangential(i,segment%HI%JsdB,k) = segment%rx_normal(i,segment%HI%jsd,k)
1670 rx_tangential(i,segment%HI%JedB,k) = segment%rx_normal(i,segment%HI%jed,k)
1671 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1672 rx_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i,j+1,k))
1675 if (segment%radiation_tan)
then
1676 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1677 rx_avg = rx_tangential(i,j,k)
1678 segment%tangential_vel(i,j,k) = (v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) / (1.0+rx_avg)
1681 if (segment%nudged_tan)
then
1682 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1684 if (rx_tangential(i,j,k) <= 0.0)
then
1685 tau = segment%Velocity_nudging_timescale_in
1687 tau = segment%Velocity_nudging_timescale_out
1689 gamma_2 = dt / (tau + dt)
1690 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
1691 gamma_2 * segment%nudged_tangential_vel(i,j,k)
1694 if (segment%radiation_grad)
then
1695 js_obc = max(segment%HI%JsdB,g%jsd+1)
1696 je_obc = min(segment%HI%JedB,g%jed-1)
1697 do k=1,nz ;
do j=js_obc,je_obc
1698 rx_avg = rx_tangential(i,j,k)
1708 segment%tangential_grad(i,j,k) = ((v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
1709 rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) / (1.0+rx_avg)
1712 if (segment%nudged_grad)
then
1713 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1715 if (rx_tangential(i,j,k) <= 0.0)
then
1716 tau = segment%Velocity_nudging_timescale_in
1718 tau = segment%Velocity_nudging_timescale_out
1720 gamma_2 = dt / (tau + dt)
1721 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
1722 gamma_2 * segment%nudged_tangential_grad(i,j,k)
1725 deallocate(rx_tangential)
1727 if (segment%oblique_tan .or. segment%oblique_grad)
then
1729 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1730 allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1731 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1733 rx_tangential(i,segment%HI%JsdB,k) = segment%rx_normal(i,segment%HI%jsd,k)
1734 rx_tangential(i,segment%HI%JedB,k) = segment%rx_normal(i,segment%HI%jed,k)
1735 ry_tangential(i,segment%HI%JsdB,k) = segment%ry_normal(i,segment%HI%jsd,k)
1736 ry_tangential(i,segment%HI%JedB,k) = segment%ry_normal(i,segment%HI%jed,k)
1737 cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
1738 cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
1739 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1740 rx_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i,j+1,k))
1741 ry_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i,j+1,k))
1742 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
1745 if (segment%oblique_tan)
then
1746 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
1747 rx_avg = rx_tangential(i,j,k)
1748 ry_avg = ry_tangential(i,j,k)
1749 cff_avg = cff_tangential(i,j,k)
1750 segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) - &
1751 (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
1755 if (segment%nudged_tan)
then
1756 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1758 if (rx_tangential(i,j,k) <= 0.0)
then
1759 tau = segment%Velocity_nudging_timescale_in
1761 tau = segment%Velocity_nudging_timescale_out
1763 gamma_2 = dt / (tau + dt)
1764 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
1765 gamma_2 * segment%nudged_tangential_vel(i,j,k)
1768 if (segment%oblique_grad)
then
1769 js_obc = max(segment%HI%JsdB,g%jsd+1)
1770 je_obc = min(segment%HI%JedB,g%jed-1)
1771 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
1772 rx_avg = rx_tangential(i,j,k)
1773 ry_avg = ry_tangential(i,j,k)
1774 cff_avg = cff_tangential(i,j,k)
1775 segment%tangential_grad(i,j,k) = ((cff_avg*(v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) &
1776 + rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) - &
1777 (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k))) / &
1781 if (segment%nudged_grad)
then
1782 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1784 if (rx_tangential(i,j,k) <= 0.0)
then
1785 tau = segment%Velocity_nudging_timescale_in
1787 tau = segment%Velocity_nudging_timescale_out
1789 gamma_2 = dt / (tau + dt)
1790 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
1791 gamma_2 * segment%nudged_tangential_grad(i,j,k)
1794 deallocate(rx_tangential)
1795 deallocate(ry_tangential)
1796 deallocate(cff_tangential)
1800 if (segment%direction == obc_direction_w)
then
1802 if (i>g%HI%IecB) cycle
1803 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
1804 if (segment%radiation)
then
1805 dhdt = u_old(i+1,j,k)-u_new(i+1,j,k)
1806 dhdx = u_new(i+1,j,k)-u_new(i+2,j,k)
1808 if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
1809 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
1810 segment%rx_normal(i,j,k) = rx_avg
1814 segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) / (1.0+rx_avg)
1817 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
1818 elseif (segment%oblique)
then
1819 dhdt = u_old(i+1,j,k)-u_new(i+1,j,k)
1820 dhdx = u_new(i+1,j,k)-u_new(i+2,j,k)
1821 if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0)
then
1822 dhdy = segment%grad_normal(j-1,1,k)
1823 elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0)
then
1826 dhdy = segment%grad_normal(j,1,k)
1828 if (dhdt*dhdx < 0.0) dhdt = 0.0
1830 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
1831 ry_new = min(cff,max(dhdt*dhdy,-cff))
1832 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
1833 ry_avg = (1.0-gamma_u)*segment%ry_normal(i,j,k) + gamma_u*ry_new
1834 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
1835 segment%rx_normal(i,j,k) = rx_avg
1836 segment%ry_normal(i,j,k) = ry_avg
1837 segment%cff_normal(i,j,k) = cff_avg
1838 segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) - &
1839 (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
1843 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
1844 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
1845 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
1846 elseif (segment%gradient)
then
1847 segment%normal_vel(i,j,k) = u_new(i+1,j,k)
1849 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then
1851 if (dhdt*dhdx <= 0.0)
then
1852 tau = segment%Velocity_nudging_timescale_in
1854 tau = segment%Velocity_nudging_timescale_out
1856 gamma_2 = dt / (tau + dt)
1857 segment%normal_vel(i,j,k) = (1 - gamma_2) * segment%normal_vel(i,j,k) + &
1858 gamma_2 * segment%nudged_normal_vel(i,j,k)
1861 if (segment%radiation_tan .or. segment%radiation_grad)
then
1863 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1865 rx_tangential(i,segment%HI%JsdB,k) = segment%rx_normal(i,segment%HI%jsd,k)
1866 rx_tangential(i,segment%HI%JedB,k) = segment%rx_normal(i,segment%HI%jed,k)
1867 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1868 rx_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i,j+1,k))
1871 if (segment%radiation_tan)
then
1872 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1873 rx_avg = rx_tangential(i,j,k)
1874 segment%tangential_vel(i,j,k) = (v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) / (1.0+rx_avg)
1877 if (segment%nudged_tan)
then
1878 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1880 if (rx_tangential(i,j,k) <= 0.0)
then
1881 tau = segment%Velocity_nudging_timescale_in
1883 tau = segment%Velocity_nudging_timescale_out
1885 gamma_2 = dt / (tau + dt)
1886 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
1887 gamma_2 * segment%nudged_tangential_vel(i,j,k)
1890 if (segment%radiation_grad)
then
1891 js_obc = max(segment%HI%JsdB,g%jsd+1)
1892 je_obc = min(segment%HI%JedB,g%jed-1)
1893 do k=1,nz ;
do j=js_obc,je_obc
1894 rx_avg = rx_tangential(i,j,k)
1904 segment%tangential_grad(i,j,k) = ((v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
1905 rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) / (1.0+rx_avg)
1908 if (segment%nudged_grad)
then
1909 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1911 if (rx_tangential(i,j,k) <= 0.0)
then
1912 tau = segment%Velocity_nudging_timescale_in
1914 tau = segment%Velocity_nudging_timescale_out
1916 gamma_2 = dt / (tau + dt)
1917 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
1918 gamma_2 * segment%nudged_tangential_grad(i,j,k)
1921 deallocate(rx_tangential)
1923 if (segment%oblique_tan .or. segment%oblique_grad)
then
1925 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1926 allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1927 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1929 rx_tangential(i,segment%HI%JsdB,k) = segment%rx_normal(i,segment%HI%jsd,k)
1930 rx_tangential(i,segment%HI%JedB,k) = segment%rx_normal(i,segment%HI%jed,k)
1931 ry_tangential(i,segment%HI%JsdB,k) = segment%ry_normal(i,segment%HI%jsd,k)
1932 ry_tangential(i,segment%HI%JedB,k) = segment%ry_normal(i,segment%HI%jed,k)
1933 cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
1934 cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
1935 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1936 rx_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i,j+1,k))
1937 ry_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i,j+1,k))
1938 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
1941 if (segment%oblique_tan)
then
1942 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
1943 rx_avg = rx_tangential(i,j,k)
1944 ry_avg = ry_tangential(i,j,k)
1945 cff_avg = cff_tangential(i,j,k)
1946 segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) - &
1947 (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
1951 if (segment%nudged_tan)
then
1952 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1954 if (rx_tangential(i,j,k) <= 0.0)
then
1955 tau = segment%Velocity_nudging_timescale_in
1957 tau = segment%Velocity_nudging_timescale_out
1959 gamma_2 = dt / (tau + dt)
1960 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
1961 gamma_2 * segment%nudged_tangential_vel(i,j,k)
1964 if (segment%oblique_grad)
then
1965 js_obc = max(segment%HI%JsdB,g%jsd+1)
1966 je_obc = min(segment%HI%JedB,g%jed-1)
1967 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
1968 rx_avg = rx_tangential(i,j,k)
1969 ry_avg = ry_tangential(i,j,k)
1970 cff_avg = cff_tangential(i,j,k)
1971 segment%tangential_grad(i,j,k) = ((cff_avg*(v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) &
1972 + rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) - &
1973 (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k))) / &
1977 if (segment%nudged_grad)
then
1978 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1980 if (rx_tangential(i,j,k) <= 0.0)
then
1981 tau = segment%Velocity_nudging_timescale_in
1983 tau = segment%Velocity_nudging_timescale_out
1985 gamma_2 = dt / (tau + dt)
1986 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
1987 gamma_2 * segment%nudged_tangential_grad(i,j,k)
1990 deallocate(rx_tangential)
1991 deallocate(ry_tangential)
1992 deallocate(cff_tangential)
1996 if (segment%direction == obc_direction_n)
then
1998 if (j<g%HI%JscB) cycle
1999 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
2000 if (segment%radiation)
then
2001 dhdt = v_old(i,j-1,k)-v_new(i,j-1,k)
2002 dhdy = v_new(i,j-1,k)-v_new(i,j-2,k)
2004 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2005 ry_avg = (1.0-gamma_v)*segment%ry_normal(i,j,k) + gamma_v*ry_new
2006 segment%ry_normal(i,j,k) = ry_avg
2010 segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) / (1.0+ry_avg)
2013 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
2014 elseif (segment%oblique)
then
2015 dhdt = v_old(i,j-1,k)-v_new(i,j-1,k)
2016 dhdy = v_new(i,j-1,k)-v_new(i,j-2,k)
2017 segment%ry_normal(i,j,k) = ry_avg
2018 if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0)
then
2019 dhdx = segment%grad_normal(i-1,1,k)
2020 elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0)
then
2023 dhdx = segment%grad_normal(i,1,k)
2025 if (dhdt*dhdy < 0.0) dhdt = 0.0
2027 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2028 rx_new = min(cff,max(dhdt*dhdx,-cff))
2029 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
2030 ry_avg = (1.0-gamma_u)*segment%ry_normal(i,j,k) + gamma_u*ry_new
2031 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2032 segment%rx_normal(i,j,k) = rx_avg
2033 segment%ry_normal(i,j,k) = ry_avg
2034 segment%cff_normal(i,j,k) = cff_avg
2035 segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) - &
2036 (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) + min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
2040 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
2041 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
2042 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2043 elseif (segment%gradient)
then
2044 segment%normal_vel(i,j,k) = v_new(i,j-1,k)
2046 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then
2048 if (dhdt*dhdy <= 0.0)
then
2049 tau = segment%Velocity_nudging_timescale_in
2051 tau = segment%Velocity_nudging_timescale_out
2053 gamma_2 = dt / (tau + dt)
2054 segment%normal_vel(i,j,k) = (1 - gamma_2) * segment%normal_vel(i,j,k) + &
2055 gamma_2 * segment%nudged_normal_vel(i,j,k)
2058 if (segment%radiation_tan .or. segment%radiation_grad)
then
2060 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2062 rx_tangential(segment%HI%IsdB,j,k) = segment%ry_normal(segment%HI%isd,j,k)
2063 rx_tangential(segment%HI%IedB,j,k) = segment%ry_normal(segment%HI%ied,j,k)
2064 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2065 rx_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i+1,j,k))
2068 if (segment%radiation_tan)
then
2069 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2070 rx_avg = rx_tangential(i,j,k)
2071 segment%tangential_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i,j-1,k)) / (1.0+rx_avg)
2074 if (segment%nudged_tan)
then
2075 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2077 if (rx_tangential(i,j,k) <= 0.0)
then
2078 tau = segment%Velocity_nudging_timescale_in
2080 tau = segment%Velocity_nudging_timescale_out
2082 gamma_2 = dt / (tau + dt)
2083 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
2084 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2087 if (segment%radiation_grad)
then
2088 is_obc = max(segment%HI%IsdB,g%isd+1)
2089 ie_obc = min(segment%HI%IedB,g%ied-1)
2090 do k=1,nz ;
do i=is_obc,ie_obc
2091 rx_avg = rx_tangential(i,j,k)
2101 segment%tangential_grad(i,j,k) = ((u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
2102 rx_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) / (1.0+rx_avg)
2105 if (segment%nudged_grad)
then
2106 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2108 if (ry_tangential(i,j,k) <= 0.0)
then
2109 tau = segment%Velocity_nudging_timescale_in
2111 tau = segment%Velocity_nudging_timescale_out
2113 gamma_2 = dt / (tau + dt)
2114 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
2115 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2118 deallocate(rx_tangential)
2120 if (segment%oblique_tan .or. segment%oblique_grad)
then
2122 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2123 allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2124 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2126 rx_tangential(segment%HI%IsdB,j,k) = segment%ry_normal(segment%HI%isd,j,k)
2127 rx_tangential(segment%HI%IedB,j,k) = segment%ry_normal(segment%HI%ied,j,k)
2128 ry_tangential(segment%HI%IsdB,j,k) = segment%rx_normal(segment%HI%isd,j,k)
2129 ry_tangential(segment%HI%IedB,j,k) = segment%rx_normal(segment%HI%ied,j,k)
2130 cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
2131 cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
2132 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2133 rx_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i+1,j,k))
2134 ry_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i+1,j,k))
2135 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
2138 if (segment%oblique_tan)
then
2139 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2140 rx_avg = rx_tangential(i,j,k)
2141 ry_avg = ry_tangential(i,j,k)
2142 cff_avg = cff_tangential(i,j,k)
2143 segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i,j-1,k)) - &
2144 (max(ry_avg,0.0)*segment%grad_tan(i,2,k) + min(ry_avg,0.0)*segment%grad_tan(i+1,2,k))) / &
2148 if (segment%nudged_tan)
then
2149 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2151 if (ry_tangential(i,j,k) <= 0.0)
then
2152 tau = segment%Velocity_nudging_timescale_in
2154 tau = segment%Velocity_nudging_timescale_out
2156 gamma_2 = dt / (tau + dt)
2157 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
2158 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2161 if (segment%oblique_grad)
then
2162 is_obc = max(segment%HI%IsdB,g%isd+1)
2163 ie_obc = min(segment%HI%IedB,g%ied-1)
2164 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2165 rx_avg = rx_tangential(i,j,k)
2166 ry_avg = ry_tangential(i,j,k)
2167 cff_avg = cff_tangential(i,j,k)
2168 segment%tangential_grad(i,j,k) = ((cff_avg*(u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) &
2169 + rx_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) - &
2170 (max(ry_avg,0.0)*segment%grad_gradient(i,2,k) + min(ry_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
2174 if (segment%nudged_grad)
then
2175 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2177 if (ry_tangential(i,j,k) <= 0.0)
then
2178 tau = segment%Velocity_nudging_timescale_in
2180 tau = segment%Velocity_nudging_timescale_out
2182 gamma_2 = dt / (tau + dt)
2183 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
2184 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2187 deallocate(rx_tangential)
2188 deallocate(ry_tangential)
2189 deallocate(cff_tangential)
2193 if (segment%direction == obc_direction_s)
then
2195 if (j>g%HI%JecB) cycle
2196 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
2197 if (segment%radiation)
then
2198 dhdt = v_old(i,j+1,k)-v_new(i,j+1,k)
2199 dhdy = v_new(i,j+1,k)-v_new(i,j+2,k)
2201 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2202 ry_avg = (1.0-gamma_v)*segment%ry_normal(i,j,k) + gamma_v*ry_new
2203 segment%ry_normal(i,j,k) = ry_avg
2207 segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) / (1.0+ry_avg)
2210 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
2211 elseif (segment%oblique)
then
2212 dhdt = v_old(i,j+1,k)-v_new(i,j+1,k)
2213 dhdy = v_new(i,j+1,k)-v_new(i,j+2,k)
2214 if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0)
then
2215 dhdx = segment%grad_normal(i-1,1,k)
2216 elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0)
then
2219 dhdx = segment%grad_normal(i,1,k)
2221 if (dhdt*dhdy < 0.0) dhdt = 0.0
2223 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2224 rx_new = min(cff,max(dhdt*dhdx,-cff))
2225 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
2226 ry_avg = (1.0-gamma_u)*segment%ry_normal(i,j,k) + gamma_u*ry_new
2227 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2228 segment%rx_normal(i,j,k) = rx_avg
2229 segment%ry_normal(i,j,k) = ry_avg
2230 segment%cff_normal(i,j,k) = cff_avg
2231 segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) - &
2232 (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) + min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
2236 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
2237 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
2238 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2239 elseif (segment%gradient)
then
2240 segment%normal_vel(i,j,k) = v_new(i,j+1,k)
2242 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then
2244 if (dhdt*dhdy <= 0.0)
then
2245 tau = segment%Velocity_nudging_timescale_in
2247 tau = segment%Velocity_nudging_timescale_out
2249 gamma_2 = dt / (tau + dt)
2250 segment%normal_vel(i,j,k) = (1 - gamma_2) * segment%normal_vel(i,j,k) + &
2251 gamma_2 * segment%nudged_normal_vel(i,j,k)
2254 if (segment%radiation_tan .or. segment%radiation_grad)
then
2256 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2258 rx_tangential(segment%HI%IsdB,j,k) = segment%ry_normal(segment%HI%isd,j,k)
2259 rx_tangential(segment%HI%IedB,j,k) = segment%ry_normal(segment%HI%ied,j,k)
2260 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2261 rx_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i+1,j,k))
2264 if (segment%radiation_tan)
then
2265 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2266 rx_avg = rx_tangential(i,j,k)
2267 segment%tangential_vel(i,j,k) = (u_new(i,j+1,k) + rx_avg*u_new(i,j+2,k)) / (1.0+rx_avg)
2270 if (segment%nudged_tan)
then
2271 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2273 if (rx_tangential(i,j,k) <= 0.0)
then
2274 tau = segment%Velocity_nudging_timescale_in
2276 tau = segment%Velocity_nudging_timescale_out
2278 gamma_2 = dt / (tau + dt)
2279 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
2280 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2283 if (segment%radiation_grad)
then
2284 is_obc = max(segment%HI%IsdB,g%isd+1)
2285 ie_obc = min(segment%HI%IedB,g%ied-1)
2286 do k=1,nz ;
do i=is_obc,ie_obc
2287 rx_avg = rx_tangential(i,j,k)
2297 segment%tangential_grad(i,j,k) = ((u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
2298 rx_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) / (1.0+rx_avg)
2301 if (segment%nudged_grad)
then
2302 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2304 if (ry_tangential(i,j,k) <= 0.0)
then
2305 tau = segment%Velocity_nudging_timescale_in
2307 tau = segment%Velocity_nudging_timescale_out
2309 gamma_2 = dt / (tau + dt)
2310 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
2311 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2314 deallocate(rx_tangential)
2316 if (segment%oblique_tan .or. segment%oblique_grad)
then
2318 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2319 allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2320 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2322 rx_tangential(segment%HI%IsdB,j,k) = segment%ry_normal(segment%HI%isd,j,k)
2323 rx_tangential(segment%HI%IedB,j,k) = segment%ry_normal(segment%HI%ied,j,k)
2324 ry_tangential(segment%HI%IsdB,j,k) = segment%rx_normal(segment%HI%isd,j,k)
2325 ry_tangential(segment%HI%IedB,j,k) = segment%rx_normal(segment%HI%ied,j,k)
2326 cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
2327 cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
2328 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2329 rx_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i+1,j,k))
2330 ry_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i+1,j,k))
2331 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
2334 if (segment%oblique_tan)
then
2335 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2336 rx_avg = rx_tangential(i,j,k)
2337 ry_avg = ry_tangential(i,j,k)
2338 cff_avg = cff_tangential(i,j,k)
2339 segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j+1,k) + rx_avg*u_new(i,j+2,k)) - &
2340 (max(ry_avg,0.0)*segment%grad_tan(i,2,k) + min(ry_avg,0.0)*segment%grad_tan(i+1,2,k))) / &
2344 if (segment%nudged_tan)
then
2345 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2347 if (ry_tangential(i,j,k) <= 0.0)
then
2348 tau = segment%Velocity_nudging_timescale_in
2350 tau = segment%Velocity_nudging_timescale_out
2352 gamma_2 = dt / (tau + dt)
2353 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
2354 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2357 if (segment%oblique_grad)
then
2358 is_obc = max(segment%HI%IsdB,g%isd+1)
2359 ie_obc = min(segment%HI%IedB,g%ied-1)
2360 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2361 rx_avg = rx_tangential(i,j,k)
2362 ry_avg = ry_tangential(i,j,k)
2363 cff_avg = cff_tangential(i,j,k)
2364 segment%tangential_grad(i,j,k) = ((cff_avg*(u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) &
2365 + rx_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) - &
2366 (max(ry_avg,0.0)*segment%grad_gradient(i,2,k) + min(ry_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
2370 if (segment%nudged_grad)
then
2371 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2373 if (ry_tangential(i,j,k) <= 0.0)
then
2374 tau = segment%Velocity_nudging_timescale_in
2376 tau = segment%Velocity_nudging_timescale_out
2378 gamma_2 = dt / (tau + dt)
2379 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
2380 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2383 deallocate(rx_tangential)
2384 deallocate(ry_tangential)
2385 deallocate(cff_tangential)
2391 call open_boundary_apply_normal_flow(obc, g, u_new, v_new)
2393 call pass_vector(u_new, v_new, g%Domain, clock=id_clock_pass)
2395 end subroutine radiation_open_bdry_conds
2398 subroutine open_boundary_apply_normal_flow(OBC, G, u, v)
2401 type(ocean_grid_type),
intent(inout) :: g
2402 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
2403 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
2405 integer :: i, j, k, n
2408 if (.not.
associated(obc))
return
2410 do n=1,obc%number_of_segments
2411 segment => obc%segment(n)
2412 if (.not. segment%on_pe)
then
2414 elseif (segment%radiation .or. segment%oblique .or. segment%gradient)
then
2415 if (segment%is_E_or_W)
then
2417 do k=1,g%ke ;
do j=segment%HI%jsd,segment%HI%jed
2418 u(i,j,k) = segment%normal_vel(i,j,k)
2420 elseif (segment%is_N_or_S)
then
2422 do k=1,g%ke ;
do i=segment%HI%isd,segment%HI%ied
2423 v(i,j,k) = segment%normal_vel(i,j,k)
2429 end subroutine open_boundary_apply_normal_flow
2432 subroutine open_boundary_zero_normal_flow(OBC, G, u, v)
2435 type(ocean_grid_type),
intent(inout) :: g
2436 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
2437 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
2439 integer :: i, j, k, n
2442 if (.not.
associated(obc))
return
2444 do n=1,obc%number_of_segments
2445 segment => obc%segment(n)
2446 if (.not. segment%on_pe)
then
2448 elseif (segment%is_E_or_W)
then
2450 do k=1,g%ke ;
do j=segment%HI%jsd,segment%HI%jed
2453 elseif (segment%is_N_or_S)
then
2455 do k=1,g%ke ;
do i=segment%HI%isd,segment%HI%ied
2461 end subroutine open_boundary_zero_normal_flow
2464 subroutine gradient_at_q_points(G, segment, uvel, vvel)
2465 type(ocean_grid_type),
intent(in) :: G
2467 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: uvel
2468 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: vvel
2471 if (.not. segment%on_pe)
return
2473 if (segment%is_E_or_W)
then
2474 if (segment%direction == obc_direction_e)
then
2477 do j=max(segment%HI%JsdB, g%HI%JsdB+1),min(segment%HI%JedB, g%HI%JedB-1)
2478 segment%grad_normal(j,1,k) = (uvel(i-1,j+1,k)-uvel(i-1,j,k)) * g%mask2dBu(i-1,j)
2479 segment%grad_normal(j,2,k) = (uvel(i,j+1,k)-uvel(i,j,k)) * g%mask2dBu(i,j)
2482 if (segment%oblique_tan)
then
2484 do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
2485 segment%grad_tan(j,1,k) = (vvel(i-1,j,k)-vvel(i-1,j-1,k)) * g%mask2dT(i-1,j)
2486 segment%grad_tan(j,2,k) = (vvel(i,j,k)-vvel(i,j-1,k)) * g%mask2dT(i,j)
2490 if (segment%oblique_grad)
then
2492 do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
2493 segment%grad_gradient(j,1,k) = (((vvel(i-1,j,k) - vvel(i-2,j,k))*g%IdxBu(i-2,j)) - &
2494 (vvel(i-1,j-1,k) - vvel(i-2,j-1,k))*g%IdxBu(i-2,j-1)) * g%mask2dCu(i-1,j)
2495 segment%grad_gradient(j,2,k) = (((vvel(i,j,k) - vvel(i-1,j,k))*g%IdxBu(i-1,j)) - &
2496 (vvel(i,j-1,k) - vvel(i-1,j-1,k))*g%IdxBu(i-1,j-1)) * g%mask2dCu(i,j)
2503 do j=max(segment%HI%JsdB, g%HI%JsdB+1),min(segment%HI%JedB, g%HI%JedB-1)
2504 segment%grad_normal(j,1,k) = (uvel(i+1,j+1,k)-uvel(i+1,j,k)) * g%mask2dBu(i+1,j)
2505 segment%grad_normal(j,2,k) = (uvel(i,j+1,k)-uvel(i,j,k)) * g%mask2dBu(i,j)
2508 if (segment%oblique_tan)
then
2510 do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
2511 segment%grad_tan(j,1,k) = (vvel(i+2,j,k)-vvel(i+2,j-1,k)) * g%mask2dT(i+2,j)
2512 segment%grad_tan(j,2,k) = (vvel(i+1,j,k)-vvel(i+1,j-1,k)) * g%mask2dT(i+1,j)
2516 if (segment%oblique_grad)
then
2518 do j=max(segment%HI%jsd, g%HI%jsd+1),min(segment%HI%jed, g%HI%jed-1)
2519 segment%grad_gradient(j,1,k) = (((vvel(i+3,j,k) - vvel(i+2,j,k))*g%IdxBu(i+2,j)) - &
2520 (vvel(i+3,j-1,k) - vvel(i+2,j-1,k))*g%IdxBu(i+2,j-1)) * g%mask2dCu(i+2,j)
2521 segment%grad_gradient(j,2,k) = (((vvel(i+2,j,k) - vvel(i+1,j,k))*g%IdxBu(i+1,j)) - &
2522 (vvel(i+2,j-1,k) - vvel(i+1,j-1,k))*g%IdxBu(i+1,j-1)) * g%mask2dCu(i+1,j)
2527 elseif (segment%is_N_or_S)
then
2528 if (segment%direction == obc_direction_n)
then
2531 do i=max(segment%HI%IsdB, g%HI%IsdB+1),min(segment%HI%IedB, g%HI%IedB-1)
2532 segment%grad_normal(i,1,k) = (vvel(i+1,j-1,k)-vvel(i,j-1,k)) * g%mask2dBu(i,j-1)
2533 segment%grad_normal(i,2,k) = (vvel(i+1,j,k)-vvel(i,j,k)) * g%mask2dBu(i,j)
2536 if (segment%oblique_tan)
then
2538 do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
2539 segment%grad_tan(i,1,k) = (uvel(i,j-1,k)-uvel(i-1,j-1,k)) * g%mask2dT(i,j-1)
2540 segment%grad_tan(i,2,k) = (uvel(i,j,k)-uvel(i-1,j,k)) * g%mask2dT(i,j)
2544 if (segment%oblique_grad)
then
2546 do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
2547 segment%grad_gradient(i,1,k) = (((uvel(i,j-1,k) - uvel(i,j-2,k))*g%IdxBu(i,j-2)) - &
2548 (uvel(i-1,j-1,k) - uvel(i-1,j-2,k))*g%IdxBu(i-1,j-2)) * g%mask2dCv(i,j-1)
2549 segment%grad_gradient(i,2,k) = (((uvel(i,j,k) - uvel(i,j-1,k))*g%IdyBu(i,j-1)) - &
2550 (uvel(i-1,j,k) - uvel(i-1,j-1,k))*g%IdyBu(i-1,j-1)) * g%mask2dCv(i,j)
2557 do i=max(segment%HI%IsdB, g%HI%IsdB+1),min(segment%HI%IedB, g%HI%IedB-1)
2558 segment%grad_normal(i,1,k) = (vvel(i+1,j+1,k)-vvel(i,j+1,k)) * g%mask2dBu(i,j+1)
2559 segment%grad_normal(i,2,k) = (vvel(i+1,j,k)-vvel(i,j,k)) * g%mask2dBu(i,j)
2562 if (segment%oblique_tan)
then
2564 do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
2565 segment%grad_tan(i,1,k) = (uvel(i,j+2,k)-uvel(i-1,j+2,k)) * g%mask2dT(i,j+2)
2566 segment%grad_tan(i,2,k) = (uvel(i,j+1,k)-uvel(i-1,j+1,k)) * g%mask2dT(i,j+1)
2570 if (segment%oblique_grad)
then
2572 do i=max(segment%HI%isd, g%HI%isd+1),min(segment%HI%ied, g%HI%ied-1)
2573 segment%grad_gradient(i,1,k) = (((uvel(i,j+3,k) - uvel(i,j+2,k))*g%IdxBu(i,j+2)) - &
2574 (uvel(i-1,j+3,k) - uvel(i-1,j+2,k))*g%IdyBu(i-1,j+2)) * g%mask2dCv(i,j+2)
2575 segment%grad_gradient(i,2,k) = (((uvel(i,j+2,k) - uvel(i,j+1,k))*g%IdxBu(i,j+1)) - &
2576 (uvel(i-1,j+2,k) - uvel(i-1,j+1,k))*g%IdyBu(i-1,j+1)) * g%mask2dCv(i,j+1)
2583 end subroutine gradient_at_q_points
2588 subroutine set_tracer_data(OBC, tv, h, G, PF, tracer_Reg)
2589 type(ocean_grid_type),
intent(inout) :: g
2591 type(thermo_var_ptrs),
intent(inout) :: tv
2592 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(inout) :: h
2593 type(param_file_type),
intent(in) :: pf
2594 type(tracer_registry_type),
pointer :: tracer_reg
2596 integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, n
2597 integer :: isd_off, jsd_off
2598 integer :: isdb, iedb, jsdb, jedb
2600 character(len=40) :: mdl =
"set_tracer_data"
2601 character(len=200) :: filename, obc_file, inputdir
2603 real :: temp_u(g%domain%niglobal+1,g%domain%njglobal)
2604 real :: temp_v(g%domain%niglobal,g%domain%njglobal+1)
2606 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2607 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2608 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
2614 if (
associated(tv%T))
then
2616 call pass_var(tv%T, g%Domain)
2617 call pass_var(tv%S, g%Domain)
2619 do n=1,obc%number_of_segments
2620 segment => obc%segment(n)
2621 if (.not. segment%on_pe) cycle
2623 if (segment%direction == obc_direction_e)
then
2625 do k=1,g%ke ;
do j=segment%HI%jsd,segment%HI%jed
2626 tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k)
2628 elseif (segment%direction == obc_direction_w)
then
2630 do k=1,g%ke ;
do j=segment%HI%jsd,segment%HI%jed
2631 tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k)
2633 elseif (segment%direction == obc_direction_n)
then
2635 do k=1,g%ke ;
do i=segment%HI%isd,segment%HI%ied
2636 tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k)
2638 elseif (segment%direction == obc_direction_s)
then
2640 do k=1,g%ke ;
do i=segment%HI%isd,segment%HI%ied
2641 tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k)
2674 end subroutine set_tracer_data
2677 function lookup_seg_field(OBC_seg,field)
2679 character(len=32),
intent(in) :: field
2680 integer :: lookup_seg_field
2685 do n=1,obc_seg%num_fields
2686 if (trim(field) == obc_seg%field(n)%name)
then
2692 end function lookup_seg_field
2696 subroutine allocate_obc_segment_data(OBC, segment)
2700 integer :: isd, ied, jsd, jed
2701 integer :: IsdB, IedB, JsdB, JedB
2702 integer :: IscB, IecB, JscB, JecB
2703 character(len=40) :: mdl =
"allocate_OBC_segment_data"
2705 isd = segment%HI%isd ; ied = segment%HI%ied
2706 jsd = segment%HI%jsd ; jed = segment%HI%jed
2707 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
2708 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
2709 iscb = segment%HI%IscB ; iecb = segment%HI%IecB
2710 jscb = segment%HI%JscB ; jecb = segment%HI%JecB
2712 if (.not. segment%on_pe)
return
2714 if (segment%is_E_or_W)
then
2716 allocate(segment%Cg(isdb:iedb,jsd:jed)); segment%Cg(:,:)=0.
2717 allocate(segment%Htot(isdb:iedb,jsd:jed)); segment%Htot(:,:)=0.0
2718 allocate(segment%h(isdb:iedb,jsd:jed,obc%ke)); segment%h(:,:,:)=0.0
2719 allocate(segment%eta(isdb:iedb,jsd:jed)); segment%eta(:,:)=0.0
2720 if (segment%radiation)
then
2721 allocate(segment%rx_normal(isdb:iedb,jsd:jed,obc%ke)); segment%rx_normal(:,:,:)=0.0
2723 allocate(segment%normal_vel(isdb:iedb,jsd:jed,obc%ke)); segment%normal_vel(:,:,:)=0.0
2724 allocate(segment%normal_vel_bt(isdb:iedb,jsd:jed)); segment%normal_vel_bt(:,:)=0.0
2725 allocate(segment%normal_trans(isdb:iedb,jsd:jed,obc%ke)); segment%normal_trans(:,:,:)=0.0
2726 if (segment%nudged)
then
2727 allocate(segment%nudged_normal_vel(isdb:iedb,jsd:jed,obc%ke)); segment%nudged_normal_vel(:,:,:)=0.0
2729 if (segment%radiation_tan .or. segment%nudged_tan .or. segment%specified_tan .or. &
2730 obc%computed_vorticity .or. obc%computed_strain)
then
2731 allocate(segment%tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_vel(:,:,:)=0.0
2733 if (segment%nudged_tan)
then
2734 allocate(segment%nudged_tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_vel(:,:,:)=0.0
2736 if (segment%nudged_grad)
then
2737 allocate(segment%nudged_tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
2739 if (obc%specified_vorticity .or. obc%specified_strain .or. segment%radiation_grad .or. &
2740 segment%oblique_grad)
then
2741 allocate(segment%tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_grad(:,:,:)=0.0
2743 if (segment%oblique)
then
2744 allocate(segment%grad_normal(jsdb:jedb,2,obc%ke)); segment%grad_normal(:,:,:) = 0.0
2745 allocate(segment%rx_normal(isdb:iedb,jsd:jed,obc%ke)); segment%rx_normal(:,:,:)=0.0
2746 allocate(segment%ry_normal(isdb:iedb,jsd:jed,obc%ke)); segment%ry_normal(:,:,:)=0.0
2747 allocate(segment%cff_normal(isdb:iedb,jsd:jed,obc%ke)); segment%cff_normal(:,:,:)=0.0
2749 if (segment%oblique_tan)
then
2750 allocate(segment%grad_tan(jsd:jed,2,obc%ke)); segment%grad_tan(:,:,:) = 0.0
2752 if (segment%oblique_grad)
then
2753 allocate(segment%grad_gradient(jsd:jed,2,obc%ke)); segment%grad_gradient(:,:,:) = 0.0
2757 if (segment%is_N_or_S)
then
2759 allocate(segment%Cg(isd:ied,jsdb:jedb)); segment%Cg(:,:)=0.
2760 allocate(segment%Htot(isd:ied,jsdb:jedb)); segment%Htot(:,:)=0.0
2761 allocate(segment%h(isd:ied,jsdb:jedb,obc%ke)); segment%h(:,:,:)=0.0
2762 allocate(segment%eta(isd:ied,jsdb:jedb)); segment%eta(:,:)=0.0
2763 if (segment%radiation)
then
2764 allocate(segment%ry_normal(isd:ied,jsdb:jedb,obc%ke)); segment%ry_normal(:,:,:)=0.0
2766 allocate(segment%normal_vel(isd:ied,jsdb:jedb,obc%ke)); segment%normal_vel(:,:,:)=0.0
2767 allocate(segment%normal_vel_bt(isd:ied,jsdb:jedb)); segment%normal_vel_bt(:,:)=0.0
2768 allocate(segment%normal_trans(isd:ied,jsdb:jedb,obc%ke)); segment%normal_trans(:,:,:)=0.0
2769 if (segment%nudged)
then
2770 allocate(segment%nudged_normal_vel(isd:ied,jsdb:jedb,obc%ke)); segment%nudged_normal_vel(:,:,:)=0.0
2772 if (segment%radiation_tan .or. segment%nudged_tan .or. segment%specified_tan .or. &
2773 obc%computed_vorticity .or. obc%computed_strain)
then
2774 allocate(segment%tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_vel(:,:,:)=0.0
2776 if (segment%nudged_tan)
then
2777 allocate(segment%nudged_tangential_vel(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_vel(:,:,:)=0.0
2779 if (segment%nudged_grad)
then
2780 allocate(segment%nudged_tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%nudged_tangential_grad(:,:,:)=0.0
2782 if (obc%specified_vorticity .or. obc%specified_strain .or. segment%radiation_grad .or. &
2783 segment%oblique_grad)
then
2784 allocate(segment%tangential_grad(isdb:iedb,jsdb:jedb,obc%ke)); segment%tangential_grad(:,:,:)=0.0
2786 if (segment%oblique)
then
2787 allocate(segment%grad_normal(isdb:iedb,2,obc%ke)); segment%grad_normal(:,:,:) = 0.0
2788 allocate(segment%rx_normal(isd:ied,jsdb:jedb,obc%ke)); segment%rx_normal(:,:,:)=0.0
2789 allocate(segment%ry_normal(isd:ied,jsdb:jedb,obc%ke)); segment%ry_normal(:,:,:)=0.0
2790 allocate(segment%cff_normal(isd:ied,jsdb:jedb,obc%ke)); segment%cff_normal(:,:,:)=0.0
2792 if (segment%oblique_tan)
then
2793 allocate(segment%grad_tan(isd:ied,2,obc%ke)); segment%grad_tan(:,:,:) = 0.0
2795 if (segment%oblique_grad)
then
2796 allocate(segment%grad_gradient(isd:ied,2,obc%ke)); segment%grad_gradient(:,:,:) = 0.0
2800 end subroutine allocate_obc_segment_data
2803 subroutine deallocate_obc_segment_data(OBC, segment)
2807 character(len=40) :: mdl =
"deallocate_OBC_segment_data"
2809 if (.not. segment%on_pe)
return
2811 if (
associated (segment%Cg))
deallocate(segment%Cg)
2812 if (
associated (segment%Htot))
deallocate(segment%Htot)
2813 if (
associated (segment%h))
deallocate(segment%h)
2814 if (
associated (segment%eta))
deallocate(segment%eta)
2815 if (
associated (segment%rx_normal))
deallocate(segment%rx_normal)
2816 if (
associated (segment%ry_normal))
deallocate(segment%ry_normal)
2817 if (
associated (segment%cff_normal))
deallocate(segment%cff_normal)
2818 if (
associated (segment%grad_normal))
deallocate(segment%grad_normal)
2819 if (
associated (segment%grad_tan))
deallocate(segment%grad_tan)
2820 if (
associated (segment%grad_gradient))
deallocate(segment%grad_gradient)
2821 if (
associated (segment%normal_vel))
deallocate(segment%normal_vel)
2822 if (
associated (segment%normal_vel_bt))
deallocate(segment%normal_vel_bt)
2823 if (
associated (segment%normal_trans))
deallocate(segment%normal_trans)
2824 if (
associated (segment%nudged_normal_vel))
deallocate(segment%nudged_normal_vel)
2825 if (
associated (segment%tangential_vel))
deallocate(segment%tangential_vel)
2826 if (
associated (segment%nudged_tangential_vel))
deallocate(segment%nudged_tangential_vel)
2827 if (
associated (segment%nudged_tangential_grad))
deallocate(segment%nudged_tangential_grad)
2828 if (
associated (segment%tangential_grad))
deallocate(segment%tangential_grad)
2829 if (
associated (segment%tr_Reg))
call segment_tracer_registry_end(segment%tr_Reg)
2832 end subroutine deallocate_obc_segment_data
2837 subroutine open_boundary_test_extern_uv(G, OBC, u, v)
2838 type(ocean_grid_type),
intent(in) :: g
2840 real,
dimension(SZIB_(G),SZJ_(G), SZK_(G)),
intent(inout) :: u
2841 real,
dimension(SZI_(G),SZJB_(G), SZK_(G)),
intent(inout) :: v
2843 integer :: i, j, k, n
2845 if (.not.
associated(obc))
return
2847 do n = 1, obc%number_of_segments
2849 if (obc%segment(n)%is_N_or_S)
then
2850 j = obc%segment(n)%HI%JsdB
2851 if (obc%segment(n)%direction == obc_direction_n)
then
2852 do i = obc%segment(n)%HI%IsdB, obc%segment(n)%HI%IedB
2853 u(i,j+1,k) = obc%silly_u
2856 do i = obc%segment(n)%HI%IsdB, obc%segment(n)%HI%IedB
2857 u(i,j,k) = obc%silly_u
2860 elseif (obc%segment(n)%is_E_or_W)
then
2861 i = obc%segment(n)%HI%IsdB
2862 if (obc%segment(n)%direction == obc_direction_e)
then
2863 do j = obc%segment(n)%HI%JsdB, obc%segment(n)%HI%JedB
2864 v(i+1,j,k) = obc%silly_u
2867 do j = obc%segment(n)%HI%JsdB, obc%segment(n)%HI%JedB
2868 v(i,j,k) = obc%silly_u
2875 end subroutine open_boundary_test_extern_uv
2880 subroutine open_boundary_test_extern_h(G, OBC, h)
2881 type(ocean_grid_type),
intent(in) :: g
2883 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(inout) :: h
2885 integer :: i, j, k, n
2887 if (.not.
associated(obc))
return
2889 do n = 1, obc%number_of_segments
2891 if (obc%segment(n)%is_N_or_S)
then
2892 j = obc%segment(n)%HI%JsdB
2893 if (obc%segment(n)%direction == obc_direction_n)
then
2894 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
2895 h(i,j+1,k) = obc%silly_h
2898 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
2899 h(i,j,k) = obc%silly_h
2902 elseif (obc%segment(n)%is_E_or_W)
then
2903 i = obc%segment(n)%HI%IsdB
2904 if (obc%segment(n)%direction == obc_direction_e)
then
2905 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
2906 h(i+1,j,k) = obc%silly_h
2909 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
2910 h(i,j,k) = obc%silly_h
2917 end subroutine open_boundary_test_extern_h
2920 subroutine update_obc_segment_data(G, GV, US, OBC, tv, h, Time)
2921 type(ocean_grid_type),
intent(in) :: g
2922 type(verticalgrid_type),
intent(in) :: gv
2923 type(unit_scale_type),
intent(in) :: us
2925 type(thermo_var_ptrs),
intent(in) :: tv
2926 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(inout) :: h
2927 type(time_type),
intent(in) :: time
2929 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed
2930 integer :: isdb, iedb, jsdb, jedb, n, m, nz
2931 character(len=40) :: mdl =
"set_OBC_segment_data"
2932 character(len=200) :: filename, obc_file, inputdir
2934 integer,
dimension(4) :: siz,siz2
2936 integer :: ni_seg, nj_seg
2938 integer :: is_obc, ie_obc, js_obc, je_obc
2939 integer :: ishift, jshift
2940 real,
dimension(:,:),
pointer :: seg_vel => null()
2941 real,
dimension(:,:),
pointer :: seg_trans => null()
2942 real,
dimension(:,:,:),
allocatable :: tmp_buffer
2943 real,
dimension(:),
allocatable :: h_stack
2944 integer :: is_obc2, js_obc2
2945 real :: net_h_src, net_h_int, scl_fac
2946 real,
pointer,
dimension(:,:) :: normal_trans_bt=>null()
2948 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2949 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2950 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
2953 if (.not.
associated(obc))
return
2955 do n = 1, obc%number_of_segments
2956 segment => obc%segment(n)
2958 if (.not. segment%on_pe) cycle
2960 ni_seg = segment%ie_obc-segment%is_obc+1
2961 nj_seg = segment%je_obc-segment%js_obc+1
2962 is_obc = max(segment%is_obc,isd-1)
2963 ie_obc = min(segment%ie_obc,ied)
2964 js_obc = max(segment%js_obc,jsd-1)
2965 je_obc = min(segment%je_obc,jed)
2978 if (segment%is_E_or_W)
then
2979 allocate(normal_trans_bt(segment%HI%IsdB:segment%HI%IedB,segment%HI%jsd:segment%HI%jed))
2980 normal_trans_bt(:,:)=0.0
2981 if (segment%direction == obc_direction_w) ishift=1
2983 do j=segment%HI%jsd,segment%HI%jed
2984 segment%Cg(i,j) = us%L_T_to_m_s*sqrt(gv%g_prime(1)*g%bathyT(i+ishift,j))
2985 segment%Htot(i,j)=0.0
2987 segment%h(i,j,k) = h(i+ishift,j,k)
2988 segment%Htot(i,j)=segment%Htot(i,j)+segment%h(i,j,k)
2992 allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB))
2993 normal_trans_bt(:,:)=0.0
2994 if (segment%direction == obc_direction_s) jshift=1
2996 do i=segment%HI%isd,segment%HI%ied
2997 segment%Cg(i,j) = us%L_T_to_m_s*sqrt(gv%g_prime(1)*g%bathyT(i,j+jshift))
2998 segment%Htot(i,j)=0.0
3000 segment%h(i,j,k) = h(i,j+jshift,k)
3001 segment%Htot(i,j)=segment%Htot(i,j)+segment%h(i,j,k)
3006 allocate(h_stack(g%ke))
3008 do m = 1,segment%num_fields
3009 if (segment%field(m)%fid > 0)
then
3010 siz(1)=
size(segment%field(m)%buffer_src,1)
3011 siz(2)=
size(segment%field(m)%buffer_src,2)
3012 siz(3)=
size(segment%field(m)%buffer_src,3)
3013 if (.not.
associated(segment%field(m)%buffer_dst))
then
3014 if (siz(3) /= segment%field(m)%nk_src)
call mom_error(fatal,
'nk_src inconsistency')
3015 if (segment%field(m)%nk_src > 1)
then
3016 if (segment%is_E_or_W)
then
3017 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
3018 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3020 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
3022 if (segment%field(m)%name ==
'U')
then
3023 allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
3024 segment%field(m)%bt_vel(:,:)=0.0
3027 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
3028 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3030 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
3032 if (segment%field(m)%name ==
'V')
then
3033 allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
3034 segment%field(m)%bt_vel(:,:)=0.0
3038 if (segment%is_E_or_W)
then
3039 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
3040 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
3042 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,1))
3044 if (segment%field(m)%name ==
'U')
then
3045 allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
3046 segment%field(m)%bt_vel(:,:)=0.0
3049 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
3050 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1))
3052 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,1))
3054 if (segment%field(m)%name ==
'V')
then
3055 allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
3056 segment%field(m)%bt_vel(:,:)=0.0
3060 segment%field(m)%buffer_dst(:,:,:)=0.0
3064 if (obc%brushcutter_mode)
then
3065 allocate(tmp_buffer(1,nj_seg*2-1,segment%field(m)%nk_src))
3067 allocate(tmp_buffer(1,nj_seg,segment%field(m)%nk_src))
3070 if (obc%brushcutter_mode)
then
3071 allocate(tmp_buffer(ni_seg*2-1,1,segment%field(m)%nk_src))
3073 allocate(tmp_buffer(ni_seg,1,segment%field(m)%nk_src))
3077 call time_interp_external(segment%field(m)%fid,time, tmp_buffer)
3078 if (obc%brushcutter_mode)
then
3079 if (segment%is_E_or_W)
then
3080 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
3081 segment%field(m)%buffer_src(is_obc,:,:) = &
3082 tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset)+1:2,:)
3084 segment%field(m)%buffer_src(is_obc,:,:) = &
3085 tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset):2,:)
3088 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
3089 segment%field(m)%buffer_src(:,js_obc,:) = &
3090 tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset)+1:2,1,:)
3092 segment%field(m)%buffer_src(:,js_obc,:) = &
3093 tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset):2,1,:)
3097 if (segment%is_E_or_W)
then
3098 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
3099 segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset+1,:)
3101 segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset,:)
3104 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
3105 segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset+1,1,:)
3107 segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset,1,:)
3111 if (segment%field(m)%nk_src > 1)
then
3112 call time_interp_external(segment%field(m)%fid_dz,time, tmp_buffer)
3113 if (obc%brushcutter_mode)
then
3114 if (segment%is_E_or_W)
then
3115 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
3116 segment%field(m)%dz_src(is_obc,:,:) = &
3117 tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset)+1:2,:)
3119 segment%field(m)%dz_src(is_obc,:,:) = &
3120 tmp_buffer(1,2*(js_obc+g%jdg_offset)+1:2*(je_obc+g%jdg_offset):2,:)
3123 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
3124 segment%field(m)%dz_src(:,js_obc,:) = &
3125 tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset)+1:2,1,:)
3127 segment%field(m)%dz_src(:,js_obc,:) = &
3128 tmp_buffer(2*(is_obc+g%idg_offset)+1:2*(ie_obc+g%idg_offset):2,1,:)
3132 if (segment%is_E_or_W)
then
3133 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
3134 segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset+1,:)
3136 segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+g%jdg_offset+1:je_obc+g%jdg_offset,:)
3139 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
3140 segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset+1,1,:)
3142 segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+g%idg_offset+1:ie_obc+g%idg_offset,1,:)
3147 if (segment%is_E_or_W)
then
3149 if (segment%direction == obc_direction_e) ishift=0
3151 if (segment%field(m)%name ==
'V' .or. segment%field(m)%name ==
'DVDX')
then
3153 do j=max(js_obc,jsd),min(je_obc,jed-1)
3156 segment%field(m)%buffer_dst(i,j,:)=0.0
3157 if (g%mask2dCu(i,j)>0. .and. g%mask2dCu(i,j+1)>0.)
then
3158 h_stack(:) = 0.5*(h(i+ishift,j,:) + h(i+ishift,j+1,:))
3159 call remapping_core_h(obc%remap_CS, &
3160 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3161 segment%field(m)%buffer_src(i,j,:), &
3162 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3163 elseif (g%mask2dCu(i,j)>0.)
then
3164 h_stack(:) = h(i+ishift,j,:)
3165 call remapping_core_h(obc%remap_CS, &
3166 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3167 segment%field(m)%buffer_src(i,j,:), &
3168 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3169 elseif (g%mask2dCu(i,j+1)>0.)
then
3170 h_stack(:) = h(i+ishift,j+1,:)
3171 call remapping_core_h(obc%remap_CS, &
3172 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3173 segment%field(m)%buffer_src(i,j,:), &
3174 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3178 do j=js_obc+1,je_obc
3181 segment%field(m)%buffer_dst(i,j,:)=0.0
3182 if (g%mask2dCu(i,j)>0.)
then
3183 net_h_src = sum( segment%field(m)%dz_src(i,j,:) )
3184 net_h_int = sum( h(i+ishift,j,:) )
3185 scl_fac = net_h_int / net_h_src
3186 call remapping_core_h(obc%remap_CS, &
3187 segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,j,:), &
3188 segment%field(m)%buffer_src(i,j,:), &
3189 g%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(i,j,:))
3195 if (segment%direction == obc_direction_n) jshift=0
3197 if (segment%field(m)%name ==
'U' .or. segment%field(m)%name ==
'DUDY')
then
3199 do i=max(is_obc,isd),min(ie_obc,ied-1)
3200 segment%field(m)%buffer_dst(i,j,:)=0.0
3201 if (g%mask2dCv(i,j)>0. .and. g%mask2dCv(i+1,j)>0.)
then
3204 h_stack(:) = 0.5*(h(i,j+jshift,:) + h(i+1,j+jshift,:))
3205 call remapping_core_h(obc%remap_CS, &
3206 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3207 segment%field(m)%buffer_src(i,j,:), &
3208 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3209 elseif (g%mask2dCv(i,j)>0.)
then
3210 h_stack(:) = h(i,j+jshift,:)
3211 call remapping_core_h(obc%remap_CS, &
3212 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3213 segment%field(m)%buffer_src(i,j,:), &
3214 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3215 elseif (g%mask2dCv(i+1,j)>0.)
then
3216 h_stack(:) = h(i+1,j+jshift,:)
3217 call remapping_core_h(obc%remap_CS, &
3218 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3219 segment%field(m)%buffer_src(i,j,:), &
3220 g%ke, h_stack, segment%field(m)%buffer_dst(i,j,:))
3224 do i=is_obc+1,ie_obc
3227 segment%field(m)%buffer_dst(i,j,:)=0.0
3228 if (g%mask2dCv(i,j)>0.)
then
3229 net_h_src = sum( segment%field(m)%dz_src(i,j,:) )
3230 net_h_int = sum( h(i,j+jshift,:) )
3231 scl_fac = net_h_int / net_h_src
3232 call remapping_core_h(obc%remap_CS, &
3233 segment%field(m)%nk_src,segment%field(m)%dz_src(i,j,:), &
3234 segment%field(m)%buffer_src(i,j,:), &
3235 g%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,j,:))
3241 segment%field(m)%buffer_dst(:,:,1)=segment%field(m)%buffer_src(:,:,1)
3243 deallocate(tmp_buffer)
3245 if (.not.
associated(segment%field(m)%buffer_dst))
then
3246 if (segment%is_E_or_W)
then
3247 if (segment%field(m)%name ==
'V')
then
3248 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3249 allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc))
3250 elseif (segment%field(m)%name ==
'U')
then
3251 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
3252 allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc))
3253 elseif (segment%field(m)%name ==
'DVDX')
then
3254 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3256 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,g%ke))
3259 if (segment%field(m)%name ==
'U')
then
3260 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3261 allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc:je_obc))
3262 elseif (segment%field(m)%name ==
'V')
then
3263 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
3264 allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc))
3265 elseif (segment%field(m)%name ==
'DUDY')
then
3266 allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,g%ke))
3268 allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,g%ke))
3271 segment%field(m)%buffer_dst(:,:,:)=segment%field(m)%value
3272 if (trim(segment%field(m)%name) ==
'U' .or. trim(segment%field(m)%name) ==
'V')
then
3273 segment%field(m)%bt_vel(:,:)=segment%field(m)%value
3278 if (trim(segment%field(m)%name) ==
'U' .or. trim(segment%field(m)%name) ==
'V')
then
3279 if (segment%field(m)%fid>0)
then
3280 if (trim(segment%field(m)%name) ==
'U' .and. segment%is_E_or_W)
then
3282 do j=js_obc+1,je_obc
3283 normal_trans_bt(i,j) = 0.0
3285 segment%normal_vel(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3286 segment%normal_trans(i,j,k) = segment%field(m)%buffer_dst(i,j,k)*segment%h(i,j,k) * &
3288 normal_trans_bt(i,j) = normal_trans_bt(i,j)+segment%normal_trans(i,j,k)
3290 segment%normal_vel_bt(i,j) = normal_trans_bt(i,j)/(max(segment%Htot(i,j),1.e-12) * &
3292 if (
associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,j,:) = segment%normal_vel(i,j,:)
3294 elseif (trim(segment%field(m)%name) ==
'V' .and. segment%is_N_or_S)
then
3296 do i=is_obc+1,ie_obc
3297 normal_trans_bt(i,j) = 0.0
3299 segment%normal_vel(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3300 segment%normal_trans(i,j,k) = segment%field(m)%buffer_dst(i,j,k)*segment%h(i,j,k) * &
3302 normal_trans_bt(i,j) = normal_trans_bt(i,j)+segment%normal_trans(i,j,k)
3304 segment%normal_vel_bt(i,j) = normal_trans_bt(i,j)/(max(segment%Htot(i,j),1.e-12) * &
3306 if (
associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,j,:) = segment%normal_vel(i,j,:)
3308 elseif (trim(segment%field(m)%name) ==
'V' .and. segment%is_E_or_W .and. &
3309 associated(segment%tangential_vel))
then
3313 segment%tangential_vel(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3315 if (
associated(segment%nudged_tangential_vel)) &
3316 segment%nudged_tangential_vel(i,j,:) = segment%tangential_vel(i,j,:)
3318 elseif (trim(segment%field(m)%name) ==
'U' .and. segment%is_N_or_S .and. &
3319 associated(segment%tangential_vel))
then
3323 segment%tangential_vel(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3325 if (
associated(segment%nudged_tangential_vel)) &
3326 segment%nudged_tangential_vel(i,j,:) = segment%tangential_vel(i,j,:)
3328 elseif (trim(segment%field(m)%name) ==
'DVDX' .and. segment%is_E_or_W .and. &
3329 associated(segment%tangential_grad))
then
3333 segment%tangential_grad(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3336 elseif (trim(segment%field(m)%name) ==
'DUDY' .and. segment%is_N_or_S .and. &
3337 associated(segment%tangential_grad))
then
3341 segment%tangential_grad(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3350 if (segment%is_E_or_W)
then
3357 if (segment%is_N_or_S)
then
3365 if (trim(segment%field(m)%name) ==
'SSH')
then
3368 segment%eta(i,j) = segment%field(m)%buffer_dst(i,j,1)
3373 if (trim(segment%field(m)%name) ==
'TEMP')
then
3374 if (
associated(segment%field(m)%buffer_dst))
then
3375 do k=1,nz;
do j=js_obc2, je_obc;
do i=is_obc2,ie_obc
3376 segment%tr_Reg%Tr(1)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3377 enddo ;
enddo ;
enddo
3378 if (.not. segment%tr_Reg%Tr(1)%is_initialized)
then
3380 do k=1,nz;
do j=js_obc2, je_obc;
do i=is_obc2,ie_obc
3381 segment%tr_Reg%Tr(1)%tres(i,j,k) = segment%tr_Reg%Tr(1)%t(i,j,k)
3382 enddo ;
enddo ;
enddo
3383 segment%tr_Reg%Tr(1)%is_initialized=.true.
3386 segment%tr_Reg%Tr(1)%OBC_inflow_conc = segment%field(m)%value
3388 elseif (trim(segment%field(m)%name) ==
'SALT')
then
3389 if (
associated(segment%field(m)%buffer_dst))
then
3390 do k=1,nz;
do j=js_obc2, je_obc;
do i=is_obc2,ie_obc
3391 segment%tr_Reg%Tr(2)%t(i,j,k) = segment%field(m)%buffer_dst(i,j,k)
3392 enddo ;
enddo ;
enddo
3393 if (.not. segment%tr_Reg%Tr(1)%is_initialized)
then
3395 do k=1,nz;
do j=js_obc2, je_obc;
do i=is_obc2,ie_obc
3396 segment%tr_Reg%Tr(2)%tres(i,j,k) = segment%tr_Reg%Tr(2)%t(i,j,k)
3397 enddo ;
enddo ;
enddo
3398 segment%tr_Reg%Tr(1)%is_initialized=.true.
3401 segment%tr_Reg%Tr(2)%OBC_inflow_conc = segment%field(m)%value
3407 deallocate(normal_trans_bt)
3411 end subroutine update_obc_segment_data
3414 subroutine register_obc(name, param_file, Reg)
3415 character(len=32),
intent(in) :: name
3416 type(param_file_type),
intent(in) :: param_file
3419 character(len=256) :: mesg
3421 if (.not.
associated(reg))
call obc_registry_init(param_file, reg)
3423 if (reg%nobc>=max_fields_)
then
3424 write(mesg,
'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for &
3425 &all the open boundaries being registered via register_OBC.")') reg%nobc+1
3426 call mom_error(fatal,
"MOM register_tracer: "//mesg)
3428 reg%nobc = reg%nobc + 1
3431 reg%OB(nobc)%name = name
3433 if (reg%locked)
call mom_error(fatal, &
3434 "MOM register_tracer was called for variable "//trim(reg%OB(nobc)%name)//&
3435 " with a locked tracer registry.")
3437 end subroutine register_obc
3440 subroutine obc_registry_init(param_file, Reg)
3441 type(param_file_type),
intent(in) :: param_file
3444 integer,
save :: init_calls = 0
3446 #include "version_variable.h"
3447 character(len=40) :: mdl =
"MOM_open_boundary"
3448 character(len=256) :: mesg
3450 if (.not.
associated(reg))
then ;
allocate(reg)
3451 else ;
return ;
endif
3456 init_calls = init_calls + 1
3457 if (init_calls > 1)
then
3458 write(mesg,
'("OBC_registry_init called ",I3, &
3459 &" times with different registry pointers.")') init_calls
3460 if (is_root_pe())
call mom_error(warning,
"MOM_open_boundary"//mesg)
3463 end subroutine obc_registry_init
3466 function register_file_obc(param_file, CS, OBC_Reg)
3467 type(param_file_type),
intent(in) :: param_file
3470 logical :: register_file_obc
3471 character(len=32) :: casename =
"OBC file"
3473 if (
associated(cs))
then
3474 call mom_error(warning,
"register_file_OBC called with an "// &
3475 "associated control structure.")
3481 call register_obc(casename, param_file, obc_reg)
3482 register_file_obc = .true.
3484 end function register_file_obc
3487 subroutine file_obc_end(CS)
3490 if (
associated(cs))
then
3493 end subroutine file_obc_end
3496 subroutine segment_tracer_registry_init(param_file, segment)
3497 type(param_file_type),
intent(in) :: param_file
3500 integer,
save :: init_calls = 0
3503 #include "version_variable.h"
3504 character(len=40) :: mdl =
"segment_tracer_registry_init"
3505 character(len=256) :: mesg
3507 if (.not.
associated(segment%tr_Reg))
then
3508 allocate(segment%tr_Reg)
3513 init_calls = init_calls + 1
3516 if (init_calls == 1)
call log_version(param_file, mdl, version,
"")
3525 end subroutine segment_tracer_registry_init
3527 subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, &
3528 OBC_scalar, OBC_array)
3529 type(verticalgrid_type),
intent(in) :: gv
3530 type(tracer_type),
target :: tr_ptr
3537 type(param_file_type),
intent(in) :: param_file
3539 real,
optional,
intent(in) :: obc_scalar
3541 logical,
optional,
intent(in) :: obc_array
3547 integer :: isd, ied, jsd, jed
3548 integer :: isdb, iedb, jsdb, jedb
3549 character(len=256) :: mesg
3551 call segment_tracer_registry_init(param_file, segment)
3553 if (segment%tr_Reg%ntseg>=max_fields_)
then
3554 write(mesg,
'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for &
3555 &all the tracers being registered via register_tracer.")') segment%tr_Reg%ntseg+1
3556 call mom_error(fatal,
"MOM register_tracer: "//mesg)
3558 segment%tr_Reg%ntseg = segment%tr_Reg%ntseg + 1
3559 ntseg = segment%tr_Reg%ntseg
3561 isd = segment%HI%isd ; ied = segment%HI%ied
3562 jsd = segment%HI%jsd ; jed = segment%HI%jed
3563 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
3564 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
3566 segment%tr_Reg%Tr(ntseg)%Tr => tr_ptr
3567 segment%tr_Reg%Tr(ntseg)%name = tr_ptr%name
3569 if (segment%tr_Reg%locked)
call mom_error(fatal, &
3570 "MOM register_tracer was called for variable "//trim(segment%tr_Reg%Tr(ntseg)%name)//&
3571 " with a locked tracer registry.")
3573 if (
present(obc_scalar)) segment%tr_Reg%Tr(ntseg)%OBC_inflow_conc = obc_scalar
3574 if (
present(obc_array))
then
3575 if (segment%is_E_or_W)
then
3576 allocate(segment%tr_Reg%Tr(ntseg)%t(isdb:iedb,jsd:jed,1:gv%ke));segment%tr_Reg%Tr(ntseg)%t(:,:,:)=0.0
3577 allocate(segment%tr_Reg%Tr(ntseg)%tres(isdb:iedb,jsd:jed,1:gv%ke));segment%tr_Reg%Tr(ntseg)%tres(:,:,:)=0.0
3578 segment%tr_Reg%Tr(ntseg)%is_initialized=.false.
3579 elseif (segment%is_N_or_S)
then
3580 allocate(segment%tr_Reg%Tr(ntseg)%t(isd:ied,jsdb:jedb,1:gv%ke));segment%tr_Reg%Tr(ntseg)%t(:,:,:)=0.0
3581 allocate(segment%tr_Reg%Tr(ntseg)%tres(isd:ied,jsdb:jedb,1:gv%ke));segment%tr_Reg%Tr(ntseg)%tres(:,:,:)=0.0
3582 segment%tr_Reg%Tr(ntseg)%is_initialized=.false.
3586 end subroutine register_segment_tracer
3589 subroutine segment_tracer_registry_end(Reg)
3595 if (
associated(reg))
then
3597 if (
associated(reg%Tr(n)%t))
deallocate(reg%Tr(n)%t)
3601 end subroutine segment_tracer_registry_end
3603 subroutine register_temp_salt_segments(GV, OBC, tr_Reg, param_file)
3604 type(verticalgrid_type),
intent(in) :: gv
3606 type(tracer_registry_type),
pointer :: tr_reg
3607 type(param_file_type),
intent(in) :: param_file
3610 integer :: isd, ied, isdb, iedb, jsd, jed, jsdb, jedb, nz, nf
3611 integer :: i, j, k, n
3612 character(len=32) :: name
3614 type(tracer_type),
pointer :: tr_ptr => null()
3616 if (.not.
associated(obc))
return
3618 do n=1, obc%number_of_segments
3619 segment=>obc%segment(n)
3620 if (.not. segment%on_pe) cycle
3622 if (
associated(segment%tr_Reg)) &
3623 call mom_error(fatal,
"register_temp_salt_segments: tracer array was previously allocated")
3626 call tracer_name_lookup(tr_reg, tr_ptr, name)
3627 call register_segment_tracer(tr_ptr, param_file, gv, segment, &
3628 obc_array=segment%temp_segment_data_exists)
3630 call tracer_name_lookup(tr_reg, tr_ptr, name)
3631 call register_segment_tracer(tr_ptr, param_file, gv, segment, &
3632 obc_array=segment%salt_segment_data_exists)
3635 end subroutine register_temp_salt_segments
3637 subroutine fill_temp_salt_segments(G, OBC, tv)
3638 type(ocean_grid_type),
intent(inout) :: g
3640 type(thermo_var_ptrs),
intent(inout) :: tv
3643 integer :: isd, ied, isdb, iedb, jsd, jed, jsdb, jedb, n, nz
3647 if (.not.
associated(obc))
return
3648 if (.not.
associated(tv%T) .and.
associated(tv%S))
return
3651 call pass_var(tv%T, g%Domain)
3652 call pass_var(tv%S, g%Domain)
3656 do n=1, obc%number_of_segments
3657 segment => obc%segment(n)
3658 if (.not. segment%on_pe) cycle
3660 isd = segment%HI%isd ; ied = segment%HI%ied
3661 jsd = segment%HI%jsd ; jed = segment%HI%jed
3662 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
3663 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
3666 if (segment%is_E_or_W)
then
3668 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
3669 if (segment%direction == obc_direction_w)
then
3670 segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i+1,j,k)
3671 segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i+1,j,k)
3673 segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j,k)
3674 segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j,k)
3679 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
3680 if (segment%direction == obc_direction_s)
then
3681 segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j+1,k)
3682 segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j+1,k)
3684 segment%tr_Reg%Tr(1)%t(i,j,k) = tv%T(i,j,k)
3685 segment%tr_Reg%Tr(2)%t(i,j,k) = tv%S(i,j,k)
3689 segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:)
3690 segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:)
3692 end subroutine fill_temp_salt_segments
3697 subroutine mask_outside_obcs(G, US, param_file, OBC)
3698 type(dyn_horgrid_type),
intent(inout) :: G
3699 type(param_file_type),
intent(in) :: param_file
3701 type(unit_scale_type),
intent(in) :: US
3704 integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n
3706 logical :: fatal_error = .false.
3708 integer,
parameter :: cin = 3, cout = 4, cland = -1, cedge = -2
3709 character(len=256) :: mesg
3711 real,
allocatable,
dimension(:,:) :: color, color2
3714 if (.not.
associated(obc))
return
3716 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
3717 units=
"m", default=0.0, scale=us%m_to_Z, do_not_log=.true.)
3719 allocate(color(g%isd:g%ied, g%jsd:g%jed)) ; color = 0
3720 allocate(color2(g%isd:g%ied, g%jsd:g%jed)) ; color2 = 0
3725 color(g%isd,j) = cedge
3726 color(g%ied,j) = cedge
3727 color2(g%isd,j) = cedge
3728 color2(g%ied,j) = cedge
3731 color(i,g%jsd) = cedge
3732 color(i,g%jed) = cedge
3733 color2(i,g%jsd) = cedge
3734 color2(i,g%jed) = cedge
3741 if (g%bathyT(i,j) <= min_depth) color(i,j) = cland
3742 if (g%bathyT(i,j) <= min_depth) color2(i,j) = cland
3746 do j=g%jsd,g%jed ;
do i=g%IsdB+1,g%IedB-1
3747 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
3748 if (color(i,j) == 0.0) color(i,j) = cout
3749 if (color(i+1,j) == 0.0) color(i+1,j) = cin
3750 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
3751 if (color(i,j) == 0.0) color(i,j) = cin
3752 if (color(i+1,j) == 0.0) color(i+1,j) = cout
3755 do j=g%JsdB+1,g%JedB-1 ;
do i=g%isd,g%ied
3756 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then
3757 if (color(i,j) == 0.0) color(i,j) = cout
3758 if (color(i,j+1) == 0.0) color(i,j+1) = cin
3759 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then
3760 if (color(i,j) == 0.0) color(i,j) = cin
3761 if (color(i,j+1) == 0.0) color(i,j+1) = cout
3765 do j=g%JsdB+1,g%JedB-1 ;
do i=g%isd,g%ied
3766 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then
3767 if (color2(i,j) == 0.0) color2(i,j) = cout
3768 if (color2(i,j+1) == 0.0) color2(i,j+1) = cin
3769 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then
3770 if (color2(i,j) == 0.0) color2(i,j) = cin
3771 if (color2(i,j+1) == 0.0) color2(i,j+1) = cout
3774 do j=g%jsd,g%jed ;
do i=g%IsdB+1,g%IedB-1
3775 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
3776 if (color2(i,j) == 0.0) color2(i,j) = cout
3777 if (color2(i+1,j) == 0.0) color2(i+1,j) = cin
3778 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
3779 if (color2(i,j) == 0.0) color2(i,j) = cin
3780 if (color2(i+1,j) == 0.0) color2(i+1,j) = cout
3785 call flood_fill(g, color, cin, cout, cland)
3786 call flood_fill2(g, color2, cin, cout, cland)
3789 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
3790 if (color(i,j) /= color2(i,j))
then
3791 fatal_error = .true.
3792 write(mesg,
'("MOM_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", &
3793 "the masking of the outside grid points.")') i, j
3794 call mom_error(warning,
"MOM register_tracer: "//mesg, all_print=.true.)
3796 if (color(i,j) == cout) g%bathyT(i,j) = min_depth
3798 if (fatal_error)
call mom_error(fatal, &
3799 "MOM_open_boundary: inconsistent OBC segments.")
3803 end subroutine mask_outside_obcs
3806 subroutine flood_fill(G, color, cin, cout, cland)
3807 type(dyn_horgrid_type),
intent(inout) :: G
3808 real,
dimension(:,:),
intent(inout) :: color
3809 integer,
intent(in) :: cin
3810 integer,
intent(in) :: cout
3811 integer,
intent(in) :: cland
3814 integer :: i, j, ncount
3817 do while (ncount > 0)
3819 do j=g%jsd+1,g%jed-1
3820 do i=g%isd+1,g%ied-1
3821 if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0)
then
3822 color(i,j) = color(i-1,j)
3825 if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0)
then
3826 color(i,j) = color(i+1,j)
3829 if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0)
then
3830 color(i,j) = color(i,j-1)
3833 if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0)
then
3834 color(i,j) = color(i,j+1)
3839 do j=g%jed-1,g%jsd+1,-1
3840 do i=g%ied-1,g%isd+1,-1
3841 if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0)
then
3842 color(i,j) = color(i-1,j)
3845 if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0)
then
3846 color(i,j) = color(i+1,j)
3849 if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0)
then
3850 color(i,j) = color(i,j-1)
3853 if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0)
then
3854 color(i,j) = color(i,j+1)
3859 call pass_var(color, g%Domain)
3860 call sum_across_pes(ncount)
3863 end subroutine flood_fill
3866 subroutine flood_fill2(G, color, cin, cout, cland)
3867 type(dyn_horgrid_type),
intent(inout) :: G
3868 real,
dimension(:,:),
intent(inout) :: color
3869 integer,
intent(in) :: cin
3870 integer,
intent(in) :: cout
3871 integer,
intent(in) :: cland
3874 integer :: i, j, ncount
3877 do while (ncount > 0)
3879 do i=g%isd+1,g%ied-1
3880 do j=g%jsd+1,g%jed-1
3881 if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0)
then
3882 color(i,j) = color(i-1,j)
3885 if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0)
then
3886 color(i,j) = color(i+1,j)
3889 if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0)
then
3890 color(i,j) = color(i,j-1)
3893 if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0)
then
3894 color(i,j) = color(i,j+1)
3899 do i=g%ied-1,g%isd+1,-1
3900 do j=g%jed-1,g%jsd+1,-1
3901 if (color(i,j) == 0.0 .and. color(i-1,j) > 0.0)
then
3902 color(i,j) = color(i-1,j)
3905 if (color(i,j) == 0.0 .and. color(i+1,j) > 0.0)
then
3906 color(i,j) = color(i+1,j)
3909 if (color(i,j) == 0.0 .and. color(i,j-1) > 0.0)
then
3910 color(i,j) = color(i,j-1)
3913 if (color(i,j) == 0.0 .and. color(i,j+1) > 0.0)
then
3914 color(i,j) = color(i,j+1)
3919 call pass_var(color, g%Domain)
3920 call sum_across_pes(ncount)
3923 end subroutine flood_fill2
3926 subroutine open_boundary_register_restarts(HI, GV, OBC_CS,restart_CSp)
3927 type(hor_index_type),
intent(in) :: hi
3928 type(verticalgrid_type),
pointer :: gv
3930 type(mom_restart_cs),
pointer :: restart_csp
3934 if (.not.
associated(obc_cs)) &
3935 call mom_error(fatal,
"open_boundary_register_restarts: Called with "//&
3936 "uninitialized OBC control structure")
3938 if (
associated(obc_cs%rx_normal) .or.
associated(obc_cs%ry_normal)) &
3939 call mom_error(fatal,
"open_boundary_register_restarts: Restart "//&
3940 "arrays were previously allocated")
3946 if (obc_cs%radiation_BCs_exist_globally .or. obc_cs%oblique_BCs_exist_globally)
then
3947 allocate(obc_cs%rx_normal(hi%isdB:hi%iedB,hi%jsd:hi%jed,gv%ke))
3948 obc_cs%rx_normal(:,:,:) = 0.0
3949 vd = var_desc(
"rx_normal",
"m s-1",
"Normal Phase Speed for EW OBCs",
'u',
'L')
3950 call register_restart_field(obc_cs%rx_normal, vd, .true., restart_csp)
3951 allocate(obc_cs%ry_normal(hi%isd:hi%ied,hi%jsdB:hi%jedB,gv%ke))
3952 obc_cs%ry_normal(:,:,:) = 0.0
3953 vd = var_desc(
"ry_normal",
"m s-1",
"Normal Phase Speed for NS OBCs",
'v',
'L')
3954 call register_restart_field(obc_cs%ry_normal, vd, .true., restart_csp)
3956 if (obc_cs%oblique_BCs_exist_globally)
then
3957 allocate(obc_cs%cff_normal(hi%IsdB:hi%IedB,hi%jsdB:hi%jedB,gv%ke))
3958 obc_cs%cff_normal(:,:,:) = 0.0
3959 vd = var_desc(
"cff_normal",
"m s-1",
"denominator for oblique OBCs",
'q',
'L')
3960 call register_restart_field(obc_cs%cff_normal, vd, .true., restart_csp)
3963 end subroutine open_boundary_register_restarts