forked from mom-ocean/MOM6
-
Notifications
You must be signed in to change notification settings - Fork 59
/
MOM_open_boundary.F90
6145 lines (5672 loc) · 298 KB
/
MOM_open_boundary.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!> Controls where open boundary conditions are applied
module MOM_open_boundary
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_array_transform, only : rotate_array, rotate_array_pair
use MOM_array_transform, only : allocate_rotated_array
use MOM_coms, only : sum_across_PEs, Set_PElist, Get_PElist, PE_here, num_PEs
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
use MOM_debugging, only : hchksum, uvchksum
use MOM_diag_mediator, only : diag_ctrl, time_type
use MOM_domains, only : pass_var, pass_vector
use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
use MOM_domains, only : To_All, EAST_FACE, NORTH_FACE, SCALAR_PAIR, CGRID_NE, CORNER
use MOM_dyn_horgrid, only : dyn_horgrid_type
use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, NOTE, is_root_pe
use MOM_file_parser, only : get_param, log_version, param_file_type, log_param
use MOM_grid, only : ocean_grid_type, hor_index_type
use MOM_interface_heights, only : thickness_to_dz
use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init
use MOM_interpolate, only : external_field
use MOM_io, only : slasher, field_size, SINGLE_FILE
use MOM_io, only : vardesc, query_vardesc, var_desc
use MOM_obsolete_params, only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char
use MOM_regridding, only : regridding_CS
use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS
use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping
use MOM_restart, only : register_restart_field, register_restart_pair
use MOM_restart, only : query_initialized, MOM_restart_CS
use MOM_string_functions, only : extract_word, remove_spaces, uppercase, lowercase
use MOM_tidal_forcing, only : astro_longitudes, astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency
use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-)
use MOM_tracer_registry, only : tracer_type, tracer_registry_type, tracer_name_lookup
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
implicit none ; private
#include <MOM_memory.h>
public open_boundary_apply_normal_flow
public open_boundary_config
public open_boundary_setup_vert
public open_boundary_init
public open_boundary_query
public open_boundary_end
public open_boundary_impose_normal_slope
public open_boundary_impose_land_mask
public radiation_open_bdry_conds
public set_tracer_data
public update_OBC_segment_data
public open_boundary_test_extern_uv
public open_boundary_test_extern_h
public open_boundary_zero_normal_flow
public parse_segment_str
public parse_segment_manifest_str
public parse_segment_data_str
public register_OBC, OBC_registry_init
public register_file_OBC, file_OBC_end
public segment_tracer_registry_init
public segment_tracer_registry_end
public register_segment_tracer
public register_temp_salt_segments
public register_obgc_segments
public fill_temp_salt_segments
public fill_obgc_segments
public set_obgc_segments_props
public setup_OBC_tracer_reservoirs
public open_boundary_register_restarts
public update_segment_tracer_reservoirs
public update_OBC_ramp
public remap_OBC_fields
public rotate_OBC_config
public rotate_OBC_init
public initialize_segment_data
public flood_fill
public flood_fill2
integer, parameter, public :: OBC_NONE = 0 !< Indicates the use of no open boundary
integer, parameter, public :: OBC_DIRECTION_N = 100 !< Indicates the boundary is an effective northern boundary
integer, parameter, public :: OBC_DIRECTION_S = 200 !< Indicates the boundary is an effective southern boundary
integer, parameter, public :: OBC_DIRECTION_E = 300 !< Indicates the boundary is an effective eastern boundary
integer, parameter, public :: OBC_DIRECTION_W = 400 !< Indicates the boundary is an effective western boundary
integer, parameter :: MAX_OBC_FIELDS = 100 !< Maximum number of data fields needed for OBC segments
!> Open boundary segment data from files (mostly).
type, public :: OBC_segment_data_type
type(external_field) :: handle !< handle from FMS associated with segment data on disk
type(external_field) :: dz_handle !< handle from FMS associated with segment thicknesses on disk
logical :: use_IO = .false. !< True if segment data is based on file input
character(len=32) :: name !< a name identifier for the segment data
character(len=8) :: genre !< an identifier for the segment data
real :: scale !< A scaling factor for converting input data to
!! the internal units of this field. For salinity this would
!! be in units of [S ppt-1 ~> 1]
real, allocatable :: buffer_src(:,:,:) !< buffer for segment data located at cell faces and on
!! the original vertical grid in the internally scaled
!! units for the field in question, such as [L T-1 ~> m s-1]
!! for a velocity or [S ~> ppt] for salinity.
integer :: nk_src !< Number of vertical levels in the source data
real, allocatable :: dz_src(:,:,:) !< vertical grid cell spacing of the incoming segment
!! data in [Z ~> m].
real, allocatable :: buffer_dst(:,:,:) !< buffer src data remapped to the target vertical grid
!! in the internally scaled units for the field in
!! question, such as [L T-1 ~> m s-1] for a velocity or
!! [S ~> ppt] for salinity.
real :: value !< A constant value for the inflow concentration if not read
!! from file, in the internal units of a field, such as [S ~> ppt]
!! for salinity.
real :: resrv_lfac_in = 1. !< The reservoir inverse length scale factor for the inward
!! direction per field [nondim]. The general 1/Lscale_in is
!! multiplied by this factor for a specific tracer.
real :: resrv_lfac_out= 1. !< The reservoir inverse length scale factor for the outward
!! direction per field [nondim]. The general 1/Lscale_out is
!! multiplied by this factor for a specific tracer.
end type OBC_segment_data_type
!> Tracer on OBC segment data structure, for putting into a segment tracer registry.
type, public :: OBC_segment_tracer_type
real, allocatable :: t(:,:,:) !< tracer concentration array in rescaled units,
!! like [S ~> ppt] for salinity.
real :: OBC_inflow_conc = 0.0 !< tracer concentration for generic inflows in rescaled units,
!! like [S ~> ppt] for salinity.
character(len=32) :: name !< tracer name used for error messages
type(tracer_type), pointer :: Tr => NULL() !< metadata describing the tracer
real, allocatable :: tres(:,:,:) !< tracer reservoir array in rescaled units,
!! like [S ~> ppt] for salinity.
real :: scale !< A scaling factor for converting the units of input
!! data, like [S ppt-1 ~> 1] for salinity.
logical :: is_initialized !< reservoir values have been set when True
integer :: ntr_index = -1 !< index of segment tracer in the global tracer registry
integer :: fd_index = -1 !< index of segment tracer in the input fields
end type OBC_segment_tracer_type
!> Registry type for tracers on segments
type, public :: segment_tracer_registry_type
integer :: ntseg = 0 !< number of registered tracer segments
type(OBC_segment_tracer_type) :: Tr(MAX_FIELDS_) !< array of registered tracers
logical :: locked = .false. !< New tracers may be registered if locked=.false.
!! When locked=.true.,no more tracers can be registered.
!! Not sure who should lock it or when...
end type segment_tracer_registry_type
!> Open boundary segment data structure. Unless otherwise noted, 2-d and 3-d arrays are discretized
!! at the same position as normal velocity points in the middle of the OBC segments.
type, public :: OBC_segment_type
logical :: Flather !< If true, applies Flather + Chapman radiation of barotropic gravity waves.
logical :: radiation !< If true, 1D Orlanksi radiation boundary conditions are applied.
!! If False, a gradient condition is applied.
logical :: radiation_tan !< If true, 1D Orlanksi radiation boundary conditions are applied to
!! tangential flows.
logical :: radiation_grad !< If true, 1D Orlanksi radiation boundary conditions are applied to
!! dudv and dvdx.
logical :: oblique !< Oblique waves supported at radiation boundary.
logical :: oblique_tan !< If true, 2D radiation boundary conditions are applied to
!! tangential flows.
logical :: oblique_grad !< If true, 2D radiation boundary conditions are applied to
!! dudv and dvdx.
logical :: nudged !< Optional supplement to radiation boundary.
logical :: nudged_tan !< Optional supplement to nudge tangential velocity.
logical :: nudged_grad !< Optional supplement to nudge normal gradient of tangential velocity.
logical :: specified !< Boundary normal velocity fixed to external value.
logical :: specified_tan !< Boundary tangential velocity fixed to external value.
logical :: specified_grad !< Boundary gradient of tangential velocity fixed to external value.
logical :: open !< Boundary is open for continuity solver, and there are no other
!! parameterized mass fluxes at the open boundary.
logical :: gradient !< Zero gradient at boundary.
logical :: values_needed !< Whether or not any external OBC fields are needed.
logical :: u_values_needed !< Whether or not external u OBC fields are needed.
logical :: uamp_values_needed !< Whether or not external u amplitude OBC fields are needed.
logical :: uphase_values_needed !< Whether or not external u phase OBC fields are needed.
logical :: v_values_needed !< Whether or not external v OBC fields are needed.
logical :: vamp_values_needed !< Whether or not external v amplitude OBC fields are needed.
logical :: vphase_values_needed !< Whether or not external v phase OBC fields are needed.
logical :: t_values_needed!< Whether or not external T OBC fields are needed.
logical :: s_values_needed!< Whether or not external S OBC fields are needed.
logical :: z_values_needed!< Whether or not external zeta OBC fields are needed.
logical :: zamp_values_needed !< Whether or not external zeta amplitude OBC fields are needed.
logical :: zphase_values_needed !< Whether or not external zeta phase OBC fields are needed.
logical :: g_values_needed!< Whether or not external gradient OBC fields are needed.
integer :: direction !< Boundary faces one of the four directions.
logical :: is_N_or_S !< True if the OB is facing North or South and exists on this PE.
logical :: is_E_or_W !< True if the OB is facing East or West and exists on this PE.
logical :: is_E_or_W_2 !< True if the OB is facing East or West anywhere.
type(OBC_segment_data_type), pointer :: field(:) => NULL() !< OBC data
integer :: num_fields !< number of OBC data fields (e.g. u_normal,u_parallel and eta for Flather)
integer :: Is_obc !< i-indices of boundary segment.
integer :: Ie_obc !< i-indices of boundary segment.
integer :: Js_obc !< j-indices of boundary segment.
integer :: Je_obc !< j-indices of boundary segment.
integer :: uamp_index !< Save where uamp is in segment%field.
integer :: uphase_index !< Save where uphase is in segment%field.
integer :: vamp_index !< Save where vamp is in segment%field.
integer :: vphase_index !< Save where vphase is in segment%field.
integer :: zamp_index !< Save where zamp is in segment%field.
integer :: zphase_index !< Save where zphase is in segment%field.
real :: Velocity_nudging_timescale_in !< Nudging timescale on inflow [T ~> s].
real :: Velocity_nudging_timescale_out !< Nudging timescale on outflow [T ~> s].
logical :: on_pe !< true if any portion of the segment is located in this PE's data domain
logical :: temp_segment_data_exists !< true if temperature data arrays are present
logical :: salt_segment_data_exists !< true if salinity data arrays are present
real, allocatable :: Cg(:,:) !< The external gravity wave speed [L T-1 ~> m s-1]
!! at OBC-points.
real, allocatable :: Htot(:,:) !< The total column thickness [H ~> m or kg m-2] at OBC-points.
real, allocatable :: dZtot(:,:) !< The total column vertical extent [Z ~> m] at OBC-points.
real, allocatable :: h(:,:,:) !< The cell thickness [H ~> m or kg m-2] at OBC-points.
real, allocatable :: normal_vel(:,:,:) !< The layer velocity normal to the OB
!! segment [L T-1 ~> m s-1].
real, allocatable :: tangential_vel(:,:,:) !< The layer velocity tangential to the OB segment
!! [L T-1 ~> m s-1], discretized at the corner points.
real, allocatable :: tangential_grad(:,:,:) !< The gradient of the velocity tangential to the OB
!! segment [T-1 ~> s-1], discretized at the corner points.
real, allocatable :: normal_trans(:,:,:) !< The layer transport normal to the OB
!! segment [H L2 T-1 ~> m3 s-1].
real, allocatable :: normal_vel_bt(:,:) !< The barotropic velocity normal to
!! the OB segment [L T-1 ~> m s-1].
real, allocatable :: SSH(:,:) !< The sea-surface elevation along the
!! segment [Z ~> m].
real, allocatable :: grad_normal(:,:,:) !< The gradient of the normal flow along the
!! segment times the grid spacing [L T-1 ~> m s-1],
!! with the first index being the corner-point index
!! along the segment, and the second index being 1 (for
!! values one point into the domain) or 2 (for values
!! along the OBC itself)
real, allocatable :: grad_tan(:,:,:) !< The gradient of the tangential flow along the
!! segment times the grid spacing [L T-1 ~> m s-1], with the
!! first index being the velocity/tracer point index along the
!! segment, and the second being 1 for the value 1.5 points
!! inside the domain and 2 for the value half a point
!! inside the domain.
real, allocatable :: grad_gradient(:,:,:) !< The gradient normal to the segment of the gradient
!! tangetial to the segment of tangential flow along the segment
!! times the grid spacing [T-1 ~> s-1], with the first
!! index being the velocity/tracer point index along the segment,
!! and the second being 1 for the value 2 points into the domain
!! and 2 for the value 1 point into the domain.
real, allocatable :: rx_norm_rad(:,:,:) !< The previous normal phase speed use for EW radiation
!! OBC, in grid points per timestep [nondim]
real, allocatable :: ry_norm_rad(:,:,:) !< The previous normal phase speed use for NS radiation
!! OBC, in grid points per timestep [nondim]
real, allocatable :: rx_norm_obl(:,:,:) !< The previous x-direction normalized radiation coefficient
!! for either EW or NS oblique OBCs [L2 T-2 ~> m2 s-2]
real, allocatable :: ry_norm_obl(:,:,:) !< The previous y-direction normalized radiation coefficient
!! for either EW or NS oblique OBCs [L2 T-2 ~> m2 s-2]
real, allocatable :: cff_normal(:,:,:) !< The denominator for oblique radiation of the normal
!! velocity [L2 T-2 ~> m2 s-2]
real, allocatable :: nudged_normal_vel(:,:,:) !< The layer velocity normal to the OB segment
!! that values should be nudged towards [L T-1 ~> m s-1].
real, allocatable :: nudged_tangential_vel(:,:,:) !< The layer velocity tangential to the OB segment
!! that values should be nudged towards [L T-1 ~> m s-1],
!! discretized at the corner (PV) points.
real, allocatable :: nudged_tangential_grad(:,:,:) !< The layer dvdx or dudy towards which nudging
!! can occur [T-1 ~> s-1].
type(segment_tracer_registry_type), pointer :: tr_Reg=> NULL()!< A pointer to the tracer registry for the segment.
type(hor_index_type) :: HI !< Horizontal index ranges
real :: Tr_InvLscale_out !< An effective inverse length scale for restoring
!! the tracer concentration in a fictitious
!! reservoir towards interior values when flow
!! is exiting the domain [L-1 ~> m-1]
real :: Tr_InvLscale_in !< An effective inverse length scale for restoring
!! the tracer concentration towards an externally
!! imposed value when flow is entering [L-1 ~> m-1]
end type OBC_segment_type
!> Open-boundary data
type, public :: ocean_OBC_type
integer :: number_of_segments = 0 !< The number of open-boundary segments.
integer :: ke = 0 !< The number of model layers
logical :: open_u_BCs_exist_globally = .false. !< True if any zonal velocity points
!! in the global domain use open BCs.
logical :: open_v_BCs_exist_globally = .false. !< True if any meridional velocity points
!! in the global domain use open BCs.
logical :: Flather_u_BCs_exist_globally = .false. !< True if any zonal velocity points
!! in the global domain use Flather BCs.
logical :: Flather_v_BCs_exist_globally = .false. !< True if any meridional velocity points
!! in the global domain use Flather BCs.
logical :: oblique_BCs_exist_globally = .false. !< True if any velocity points
!! in the global domain use oblique BCs.
logical :: nudged_u_BCs_exist_globally = .false. !< True if any velocity points in the
!! global domain use nudged BCs.
logical :: nudged_v_BCs_exist_globally = .false. !< True if any velocity points in the
!! global domain use nudged BCs.
logical :: specified_u_BCs_exist_globally = .false. !< True if any zonal velocity points
!! in the global domain use specified BCs.
logical :: specified_v_BCs_exist_globally = .false. !< True if any meridional velocity points
!! in the global domain use specified BCs.
logical :: radiation_BCs_exist_globally = .false. !< True if radiations BCs are in use anywhere.
logical :: user_BCs_set_globally = .false. !< True if any OBC_USER_CONFIG is set
!! for input from user directory.
logical :: update_OBC = .false. !< Is OBC data time-dependent
logical :: update_OBC_seg_data = .false. !< Is it the time for OBC segment data update for fields that
!! require less frequent update
logical :: needs_IO_for_data = .false. !< Is any i/o needed for OBCs on the current PE
logical :: any_needs_IO_for_data = .false. !< Is any i/o needed for OBCs globally
logical :: some_need_no_IO_for_data = .false. !< Are there any PEs with OBCs that do not need i/o.
logical :: zero_vorticity = .false. !< If True, sets relative vorticity to zero on open boundaries.
logical :: freeslip_vorticity = .false. !< If True, sets normal gradient of tangential velocity to zero
!! in the relative vorticity on open boundaries.
logical :: computed_vorticity = .false. !< If True, uses external data for tangential velocity
!! in the relative vorticity on open boundaries.
logical :: specified_vorticity = .false. !< If True, uses external data for tangential velocity
!! gradients in the relative vorticity on open boundaries.
logical :: zero_strain = .false. !< If True, sets strain to zero on open boundaries.
logical :: freeslip_strain = .false. !< If True, sets normal gradient of tangential velocity to zero
!! in the strain on open boundaries.
logical :: computed_strain = .false. !< If True, uses external data for tangential velocity to compute
!! normal gradient in the strain on open boundaries.
logical :: specified_strain = .false. !< If True, uses external data for tangential velocity gradients
!! to compute strain on open boundaries.
logical :: zero_biharmonic = .false. !< If True, zeros the Laplacian of flow on open boundaries for
!! use in the biharmonic viscosity term.
logical :: brushcutter_mode = .false. !< If True, read data on supergrid.
logical, allocatable :: tracer_x_reservoirs_used(:) !< Dimensioned by the number of tracers, set globally,
!! true for those with x reservoirs (needed for restarts).
logical, allocatable :: tracer_y_reservoirs_used(:) !< Dimensioned by the number of tracers, set globally,
!! true for those with y reservoirs (needed for restarts).
integer :: ntr = 0 !< number of tracers
integer :: n_tide_constituents = 0 !< Number of tidal constituents to add to the boundary.
logical :: add_tide_constituents = .false. !< If true, add tidal constituents to the boundary elevation
!! and velocity. Will be set to true if n_tide_constituents > 0.
character(len=2), allocatable, dimension(:) :: tide_names !< Names of tidal constituents to add to the boundary data.
real, allocatable, dimension(:) :: tide_frequencies !< Angular frequencies of chosen tidal constituents [T-1 ~> s-1].
real, allocatable, dimension(:) :: tide_eq_phases !< Equilibrium phases of chosen tidal constituents [rad].
real, allocatable, dimension(:) :: tide_fn !< Amplitude modulation of boundary tides by nodal cycle [nondim].
real, allocatable, dimension(:) :: tide_un !< Phase modulation of boundary tides by nodal cycle [rad].
logical :: add_eq_phase = .false. !< If true, add the equilibrium phase argument
!! to the specified boundary tidal phase.
logical :: add_nodal_terms = .false. !< If true, insert terms for the 18.6 year modulation when
!! calculating tidal boundary conditions.
type(time_type) :: time_ref !< Reference date (t = 0) for tidal forcing.
type(astro_longitudes) :: tidal_longitudes !< Lunar and solar longitudes used to calculate tidal forcing.
! Properties of the segments used.
type(OBC_segment_type), allocatable :: segment(:) !< List of segment objects.
! Which segment object describes the current point.
integer, allocatable :: segnum_u(:,:) !< Segment number of u-points.
integer, allocatable :: segnum_v(:,:) !< Segment number of v-points.
! Keep the OBC segment properties for external BGC tracers
type(external_tracers_segments_props), pointer :: obgc_segments_props => NULL() !< obgc segment properties
integer :: num_obgc_tracers = 0 !< The total number of obgc tracers
! The following parameters are used in the baroclinic radiation code:
real :: gamma_uv !< The relative weighting for the baroclinic radiation
!! velocities (or speed of characteristics) at the
!! new time level (1) or the running mean (0) for velocities [nondim].
!! Valid values range from 0 to 1, with a default of 0.3.
real :: rx_max !< The maximum magnitude of the baroclinic radiation velocity (or speed of
!! characteristics) in units of grid points per timestep [nondim].
logical :: OBC_pe !< Is there an open boundary on this tile?
type(remapping_CS), pointer :: remap_z_CS=> NULL() !< ALE remapping control structure for
!! z-space data on segments
type(remapping_CS), pointer :: remap_h_CS=> NULL() !< ALE remapping control structure for
!! thickness-based fields on segments
type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries
real, allocatable :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of
!! grid points per timestep [nondim]
real, allocatable :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of
!! grid points per timestep [nondim]
real, allocatable :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared
!! at u points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared
!! at u points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared
!! at v points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared
!! at v points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation
!! rates at u points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation
!! rates at v points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc]
real, allocatable :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc]
logical :: debug !< If true, write verbose checksums for debugging purposes.
real :: silly_h !< A silly value of thickness outside of the domain that can be used to test
!! the independence of the OBCs to this external data [Z ~> m].
real :: silly_u !< A silly value of velocity outside of the domain that can be used to test
!! the independence of the OBCs to this external data [L T-1 ~> m s-1].
logical :: ramp = .false. !< If True, ramp from zero to the external values for SSH.
logical :: ramping_is_activated = .false. !< True if the ramping has been initialized
real :: ramp_timescale !< If ramp is True, use this timescale for ramping [T ~> s].
real :: trunc_ramp_time !< If ramp is True, time after which ramp is done [T ~> s].
real :: ramp_value !< If ramp is True, where we are on the ramp from
!! zero to one [nondim].
type(time_type) :: ramp_start_time !< Time when model was started.
integer :: remap_answer_date !< The vintage of the order of arithmetic and expressions to use
!! for remapping. Values below 20190101 recover the remapping
!! answers from 2018, while higher values use more robust
!! forms of the same remapping expressions.
logical :: check_reconstruction !< Flag for remapping to run checks on reconstruction
logical :: check_remapping !< Flag for remapping to run internal checks
logical :: force_bounds_in_subcell !< Flag for remapping to hide overshoot using bounds
logical :: om4_remap_via_sub_cells !< If true, use the OM4 remapping algorithm
character(40) :: remappingScheme !< String selecting the vertical remapping scheme
type(group_pass_type) :: pass_oblique !< Structure for group halo pass
end type ocean_OBC_type
!> Control structure for open boundaries that read from files.
!! Probably lots to update here.
type, public :: file_OBC_CS ; private
real :: tide_flow = 3.0e6 !< Placeholder for now..., perhaps in [m3 s-1]?
end type file_OBC_CS
!> Type to carry something (what??) for the OBC registry.
type, public :: OBC_struct_type
character(len=32) :: name !< OBC name used for error messages
end type OBC_struct_type
!> Type to carry basic OBC information needed for updating values.
type, public :: OBC_registry_type
integer :: nobc = 0 !< number of registered open boundary types.
type(OBC_struct_type) :: OB(MAX_FIELDS_) !< array of registered boundary types.
logical :: locked = .false. !< New OBC types may be registered if locked=.false.
!! When locked=.true.,no more boundaries can be registered.
end type OBC_registry_type
!> Type to carry OBC information needed for setting segments for OBGC tracers
type, private :: external_tracers_segments_props
type(external_tracers_segments_props), pointer :: next => NULL() !< pointer to the next node
character(len=128) :: tracer_name !< tracer name
character(len=128) :: tracer_src_file !< tracer source file for BC
character(len=128) :: tracer_src_field !< name of the field in source file to extract BC
real :: lfac_in !< multiplicative factor for inbound tracer reservoir length scale [nondim]
real :: lfac_out !< multiplicative factor for outbound tracer reservoir length scale [nondim]
end type external_tracers_segments_props
integer :: id_clock_pass !< A CPU time clock
character(len=40) :: mdl = "MOM_open_boundary" !< This module's name.
contains
!> Enables OBC module and reads configuration parameters
!> This routine is called from MOM_initialize_fixed which
!> occurs before the initialization of the vertical coordinate
!> and ALE_init. Therefore segment data are not fully initialized
!> here. The remainder of the segment data are initialized in a
!> later call to update_open_boundary_data
subroutine open_boundary_config(G, US, param_file, OBC)
type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(param_file_type), intent(in) :: param_file !< Parameter file handle
type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
! Local variables
integer :: l ! For looping over segments
logical :: debug, debug_OBC, mask_outside, reentrant_x, reentrant_y
character(len=15) :: segment_param_str ! The run-time parameter name for each segment
character(len=1024) :: segment_str ! The contents (rhs) for parameter "segment_param_str"
character(len=200) :: config1 ! String for OBC_USER_CONFIG
real :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m]
integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
logical :: check_remapping, force_bounds_in_subcell
logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
! This include declares and sets the variable "version".
# include "version_variable.h"
allocate(OBC)
call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", OBC%number_of_segments, &
default=0, do_not_log=.true.)
call log_version(param_file, mdl, version, &
"Controls where open boundaries are located, what kind of boundary condition "//&
"to impose, and what data to apply, if any.", &
all_default=(OBC%number_of_segments<=0))
call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", OBC%number_of_segments, &
"The number of open boundary segments.", &
default=0)
call get_param(param_file, mdl, "OBC_USER_CONFIG", config1, &
"A string that sets how the open boundary conditions are "//&
" configured: \n", default="none", do_not_log=.true.)
call get_param(param_file, mdl, "NK", OBC%ke, &
"The number of model layers", default=0, do_not_log=.true.)
if (config1 /= "none" .and. config1 /= "dyed_obcs") OBC%user_BCs_set_globally = .true.
if (OBC%number_of_segments > 0) then
call get_param(param_file, mdl, "OBC_ZERO_VORTICITY", OBC%zero_vorticity, &
"If true, sets relative vorticity to zero on open boundaries.", &
default=.false.)
call get_param(param_file, mdl, "OBC_FREESLIP_VORTICITY", OBC%freeslip_vorticity, &
"If true, sets the normal gradient of tangential velocity to "//&
"zero in the relative vorticity on open boundaries. This cannot "//&
"be true if another OBC_XXX_VORTICITY option is True.", default=.true.)
call get_param(param_file, mdl, "OBC_COMPUTED_VORTICITY", OBC%computed_vorticity, &
"If true, uses the external values of tangential velocity "//&
"in the relative vorticity on open boundaries. This cannot "//&
"be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
call get_param(param_file, mdl, "OBC_SPECIFIED_VORTICITY", OBC%specified_vorticity, &
"If true, uses the external values of tangential velocity "//&
"in the relative vorticity on open boundaries. This cannot "//&
"be true if another OBC_XXX_VORTICITY option is True.", default=.false.)
if ((OBC%zero_vorticity .and. OBC%freeslip_vorticity) .or. &
(OBC%zero_vorticity .and. OBC%computed_vorticity) .or. &
(OBC%zero_vorticity .and. OBC%specified_vorticity) .or. &
(OBC%freeslip_vorticity .and. OBC%computed_vorticity) .or. &
(OBC%freeslip_vorticity .and. OBC%specified_vorticity) .or. &
(OBC%computed_vorticity .and. OBC%specified_vorticity)) &
call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"//&
"Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"//&
"and OBC_IMPORTED_VORTICITY can be True at once.")
call get_param(param_file, mdl, "OBC_ZERO_STRAIN", OBC%zero_strain, &
"If true, sets the strain used in the stress tensor to zero on open boundaries.", &
default=.false.)
call get_param(param_file, mdl, "OBC_FREESLIP_STRAIN", OBC%freeslip_strain, &
"If true, sets the normal gradient of tangential velocity to "//&
"zero in the strain use in the stress tensor on open boundaries. This cannot "//&
"be true if another OBC_XXX_STRAIN option is True.", default=.true.)
call get_param(param_file, mdl, "OBC_COMPUTED_STRAIN", OBC%computed_strain, &
"If true, sets the normal gradient of tangential velocity to "//&
"zero in the strain use in the stress tensor on open boundaries. This cannot "//&
"be true if another OBC_XXX_STRAIN option is True.", default=.false.)
call get_param(param_file, mdl, "OBC_SPECIFIED_STRAIN", OBC%specified_strain, &
"If true, sets the normal gradient of tangential velocity to "//&
"zero in the strain use in the stress tensor on open boundaries. This cannot "//&
"be true if another OBC_XXX_STRAIN option is True.", default=.false.)
if ((OBC%zero_strain .and. OBC%freeslip_strain) .or. &
(OBC%zero_strain .and. OBC%computed_strain) .or. &
(OBC%zero_strain .and. OBC%specified_strain) .or. &
(OBC%freeslip_strain .and. OBC%computed_strain) .or. &
(OBC%freeslip_strain .and. OBC%specified_strain) .or. &
(OBC%computed_strain .and. OBC%specified_strain)) &
call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config: \n"//&
"Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN, OBC_COMPUTED_STRAIN \n"//&
"and OBC_IMPORTED_STRAIN can be True at once.")
call get_param(param_file, mdl, "OBC_ZERO_BIHARMONIC", OBC%zero_biharmonic, &
"If true, zeros the Laplacian of flow on open boundaries in the biharmonic "//&
"viscosity term.", default=.false.)
call get_param(param_file, mdl, "MASK_OUTSIDE_OBCS", mask_outside, &
"If true, set the areas outside open boundaries to be land.", &
default=.false.)
call get_param(param_file, mdl, "RAMP_OBCS", OBC%ramp, &
"If true, ramps from zero to the external values over time, with"//&
"a ramping timescale given by RAMP_TIMESCALE. Ramping SSH only so far", &
default=.false.)
call get_param(param_file, mdl, "OBC_RAMP_TIMESCALE", OBC%ramp_timescale, &
"If RAMP_OBCS is true, this sets the ramping timescale.", &
units="days", default=1.0, scale=86400.0*US%s_to_T)
call get_param(param_file, mdl, "OBC_TIDE_N_CONSTITUENTS", OBC%n_tide_constituents, &
"Number of tidal constituents being added to the open boundary.", &
default=0)
if (OBC%n_tide_constituents > 0) then
OBC%add_tide_constituents = .true.
else
OBC%add_tide_constituents = .false.
endif
call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
! This extra get_param call is to enable logging if either DEBUG or DEBUG_OBC are true.
call get_param(param_file, mdl, "DEBUG_OBC", debug_OBC, default=debug)
call get_param(param_file, mdl, "DEBUG_OBC", OBC%debug, &
"If true, do additional calls to help debug the performance "//&
"of the open boundary condition code.", &
default=debug, do_not_log=.not.(debug_OBC.or.debug), debuggingParam=.true.)
call get_param(param_file, mdl, "OBC_SILLY_THICK", OBC%silly_h, &
"A silly value of thicknesses used outside of open boundary "//&
"conditions for debugging.", units="m", default=0.0, scale=US%m_to_Z, &
do_not_log=.not.OBC%debug, debuggingParam=.true.)
call get_param(param_file, mdl, "OBC_SILLY_VEL", OBC%silly_u, &
"A silly value of velocities used outside of open boundary "//&
"conditions for debugging.", units="m/s", default=0.0, scale=US%m_s_to_L_T, &
do_not_log=.not.OBC%debug, debuggingParam=.true.)
reentrant_x = .false.
call get_param(param_file, mdl, "REENTRANT_X", reentrant_x, default=.true.)
reentrant_y = .false.
call get_param(param_file, mdl, "REENTRANT_Y", reentrant_y, default=.false.)
! Allocate everything
allocate(OBC%segment(1:OBC%number_of_segments))
do l=1,OBC%number_of_segments
OBC%segment(l)%Flather = .false.
OBC%segment(l)%radiation = .false.
OBC%segment(l)%radiation_tan = .false.
OBC%segment(l)%radiation_grad = .false.
OBC%segment(l)%oblique = .false.
OBC%segment(l)%oblique_tan = .false.
OBC%segment(l)%oblique_grad = .false.
OBC%segment(l)%nudged = .false.
OBC%segment(l)%nudged_tan = .false.
OBC%segment(l)%nudged_grad = .false.
OBC%segment(l)%specified = .false.
OBC%segment(l)%specified_tan = .false.
OBC%segment(l)%specified_grad = .false.
OBC%segment(l)%open = .false.
OBC%segment(l)%gradient = .false.
OBC%segment(l)%values_needed = .false.
OBC%segment(l)%u_values_needed = .false.
OBC%segment(l)%uamp_values_needed = OBC%add_tide_constituents
OBC%segment(l)%uphase_values_needed = OBC%add_tide_constituents
OBC%segment(l)%v_values_needed = .false.
OBC%segment(l)%vamp_values_needed = OBC%add_tide_constituents
OBC%segment(l)%vphase_values_needed = OBC%add_tide_constituents
OBC%segment(l)%t_values_needed = .false.
OBC%segment(l)%s_values_needed = .false.
OBC%segment(l)%z_values_needed = .false.
OBC%segment(l)%zamp_values_needed = OBC%add_tide_constituents
OBC%segment(l)%zphase_values_needed = OBC%add_tide_constituents
OBC%segment(l)%g_values_needed = .false.
OBC%segment(l)%direction = OBC_NONE
OBC%segment(l)%is_N_or_S = .false.
OBC%segment(l)%is_E_or_W = .false.
OBC%segment(l)%is_E_or_W_2 = .false.
OBC%segment(l)%Velocity_nudging_timescale_in = 0.0
OBC%segment(l)%Velocity_nudging_timescale_out = 0.0
OBC%segment(l)%num_fields = 0
enddo
allocate(OBC%segnum_u(G%IsdB:G%IedB,G%jsd:G%jed), source=OBC_NONE)
allocate(OBC%segnum_v(G%isd:G%ied,G%JsdB:G%JedB), source=OBC_NONE)
do l = 1, OBC%number_of_segments
write(segment_param_str(1:15),"('OBC_SEGMENT_',i3.3)") l
call get_param(param_file, mdl, segment_param_str, segment_str, &
"Documentation needs to be dynamic?????", &
fail_if_missing=.true.)
segment_str = remove_spaces(segment_str)
if (segment_str(1:2) == 'I=') then
call setup_u_point_obc(OBC, G, US, segment_str, l, param_file, reentrant_y)
elseif (segment_str(1:2) == 'J=') then
call setup_v_point_obc(OBC, G, US, segment_str, l, param_file, reentrant_x)
else
call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config: "//&
"Unable to interpret "//segment_param_str//" = "//trim(segment_str))
endif
enddo
! Moved this earlier because time_interp_external_init needs to be called
! before anything that uses time_interp_external (such as initialize_segment_data)
if (OBC%specified_u_BCs_exist_globally .or. OBC%specified_v_BCs_exist_globally .or. &
OBC%open_u_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then
! Need this for ocean_only mode boundary interpolation.
call time_interp_external_init()
endif
! if (open_boundary_query(OBC, needs_ext_seg_data=.true.)) &
! call initialize_segment_data(G, OBC, param_file)
if (open_boundary_query(OBC, apply_open_OBC=.true.)) then
call get_param(param_file, mdl, "OBC_RADIATION_MAX", OBC%rx_max, &
"The maximum magnitude of the baroclinic radiation velocity (or speed of "//&
"characteristics), in gridpoints per timestep. This is only "//&
"used if one of the open boundary segments is using Orlanski.", &
units="nondim", default=1.0)
call get_param(param_file, mdl, "OBC_RAD_VEL_WT", OBC%gamma_uv, &
"The relative weighting for the baroclinic radiation "//&
"velocities (or speed of characteristics) at the new "//&
"time level (1) or the running mean (0) for velocities. "//&
"Valid values range from 0 to 1. This is only used if "//&
"one of the open boundary segments is using Orlanski.", &
units="nondim", default=0.3)
endif
Lscale_in = 0.
Lscale_out = 0.
if (open_boundary_query(OBC, apply_open_OBC=.true.)) then
call get_param(param_file, mdl, "OBC_TRACER_RESERVOIR_LENGTH_SCALE_OUT ", Lscale_out, &
"An effective length scale for restoring the tracer concentration "//&
"at the boundaries to externally imposed values when the flow "//&
"is exiting the domain.", units="m", default=0.0, scale=US%m_to_L)
call get_param(param_file, mdl, "OBC_TRACER_RESERVOIR_LENGTH_SCALE_IN ", Lscale_in, &
"An effective length scale for restoring the tracer concentration "//&
"at the boundaries to values from the interior when the flow "//&
"is entering the domain.", units="m", default=0.0, scale=US%m_to_L)
endif
if (mask_outside) call mask_outside_OBCs(G, US, param_file, OBC)
! All tracers are using the same restoring length scale for now, but we may want to make this
! tracer-specific in the future for example, in cases where certain tracers are poorly constrained
! by data while others are well constrained - MJH.
do l = 1, OBC%number_of_segments
OBC%segment(l)%Tr_InvLscale_in = 0.0
if (Lscale_in>0.) OBC%segment(l)%Tr_InvLscale_in = 1.0/Lscale_in
OBC%segment(l)%Tr_InvLscale_out = 0.0
if (Lscale_out>0.) OBC%segment(l)%Tr_InvLscale_out = 1.0/Lscale_out
enddo
call get_param(param_file, mdl, "REMAPPING_SCHEME", OBC%remappingScheme, &
default=remappingDefaultScheme, do_not_log=.true.)
call get_param(param_file, mdl, "OBC_REMAPPING_SCHEME", OBC%remappingScheme, &
"This sets the reconstruction scheme used "//&
"for OBC vertical remapping for all variables. "//&
"It can be one of the following schemes: \n"//&
trim(remappingSchemesDoc), default=OBC%remappingScheme)
call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", OBC%check_reconstruction, &
"If true, cell-by-cell reconstructions are checked for "//&
"consistency and if non-monotonicity or an inconsistency is "//&
"detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", OBC%check_remapping, &
"If true, the results of remapping are checked for "//&
"conservation and new extrema and if an inconsistency is "//&
"detected then a FATAL error is issued.", default=.false.,do_not_log=.true.)
call get_param(param_file, mdl, "BRUSHCUTTER_MODE", OBC%brushcutter_mode, &
"If true, read external OBC data on the supergrid.", &
default=.false.)
call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", OBC%force_bounds_in_subcell, &
"If true, the values on the intermediate grid used for remapping "//&
"are forced to be bounded, which might not be the case due to "//&
"round off.", default=.false.,do_not_log=.true.)
call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
"This sets the default value for the various _ANSWER_DATE parameters.", &
default=99991231)
call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", OBC%remap_answer_date, &
"The vintage of the expressions and order of arithmetic to use for remapping. "//&
"Values below 20190101 result in the use of older, less accurate expressions "//&
"that were in use at the end of 2018. Higher values result in the use of more "//&
"robust and accurate forms of mathematically equivalent expressions.", &
default=default_answer_date)
call get_param(param_file, mdl, "OBC_REMAPPING_USE_OM4_SUBCELLS", OBC%om4_remap_via_sub_cells, &
"If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//&
"See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
"We recommend setting this option to false.", default=.true.)
endif ! OBC%number_of_segments > 0
! Safety check
if ((OBC%open_u_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) .and. &
.not.G%symmetric ) call MOM_error(FATAL, &
"MOM_open_boundary, open_boundary_config: "//&
"Symmetric memory must be used when using Flather OBCs.")
! Need to do this last, because it depends on time_interp_external_init having already been called
if (OBC%add_tide_constituents) then
call initialize_obc_tides(OBC, US, param_file)
! Tide update is done within update_OBC_segment_data, so this should be true if tides are included.
OBC%update_OBC = .true.
endif
if (.not.(OBC%specified_u_BCs_exist_globally .or. OBC%specified_v_BCs_exist_globally .or. &
OBC%open_u_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally)) then
! No open boundaries have been requested
call open_boundary_dealloc(OBC)
endif
end subroutine open_boundary_config
!> Setup vertical remapping for open boundaries
subroutine open_boundary_setup_vert(GV, US, OBC)
type(verticalGrid_type), intent(in) :: GV !< Container for vertical grid information
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure
! Local variables
real :: dz_neglect, dz_neglect_edge ! Small thicknesses in vertical height units [Z ~> m]
if (associated(OBC)) then
if (OBC%number_of_segments > 0) then
if (GV%Boussinesq .and. (OBC%remap_answer_date < 20190101)) then
dz_neglect = US%m_to_Z * 1.0e-30 ; dz_neglect_edge = US%m_to_Z * 1.0e-10
elseif (GV%semi_Boussinesq .and. (OBC%remap_answer_date < 20190101)) then
dz_neglect = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-30 ; dz_neglect_edge = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-10
else
dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff
endif
allocate(OBC%remap_z_CS)
call initialize_remapping(OBC%remap_z_CS, OBC%remappingScheme, boundary_extrapolation=.false., &
check_reconstruction=OBC%check_reconstruction, check_remapping=OBC%check_remapping, &
om4_remap_via_sub_cells=OBC%om4_remap_via_sub_cells, &
force_bounds_in_subcell=OBC%force_bounds_in_subcell, answer_date=OBC%remap_answer_date, &
h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge)
allocate(OBC%remap_h_CS)
call initialize_remapping(OBC%remap_h_CS, OBC%remappingScheme, boundary_extrapolation=.false., &
check_reconstruction=OBC%check_reconstruction, check_remapping=OBC%check_remapping, &
om4_remap_via_sub_cells=OBC%om4_remap_via_sub_cells, &
force_bounds_in_subcell=OBC%force_bounds_in_subcell, answer_date=OBC%remap_answer_date, &
h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff)
endif
endif
end subroutine open_boundary_setup_vert
!> Allocate space for reading OBC data from files. It sets up the required vertical
!! remapping. In the process, it does funky stuff with the MPI processes.
subroutine initialize_segment_data(G, GV, US, OBC, PF)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Container for vertical grid information
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(ocean_OBC_type), target, intent(inout) :: OBC !< Open boundary control structure
type(param_file_type), intent(in) :: PF !< Parameter file handle
integer :: n, m, num_fields, mm
character(len=1024) :: segstr
character(len=256) :: filename
character(len=20) :: segnam, suffix
character(len=32) :: fieldname
real :: value ! A value that is parsed from the segment data string [various units]
character(len=32), dimension(MAX_OBC_FIELDS) :: fields ! segment field names
character(len=128) :: inputdir
type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list
character(len=256) :: mesg ! Message for error messages.
integer, dimension(4) :: siz,siz2
integer :: is, ie, js, je
integer :: isd, ied, jsd, jed
integer :: IsdB, IedB, JsdB, JedB
integer, dimension(:), allocatable :: saved_pelist
integer :: current_pe
integer, dimension(1) :: single_pelist
type(external_tracers_segments_props), pointer :: obgc_segments_props_list =>NULL()
!will be able to dynamically switch between sub-sampling refined grid data or model grid
integer :: IO_needs(3) ! Sums to determine global OBC data use and update patterns.
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
! There is a problem with the order of the OBC initialization
! with respect to ALE_init. Currently handling this by copying the
! param file so that I can use it later in step_MOM in order to finish
! initializing segments on the first step.
call get_param(PF, mdl, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
if (OBC%user_BCs_set_globally) return
! Try this here just for the documentation. It is repeated below.
do n=1, OBC%number_of_segments
write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n
call get_param(PF, mdl, segnam, segstr, 'OBC segment docs')
enddo
!< temporarily disable communication in order to read segment data independently
allocate(saved_pelist(0:num_PEs()-1))
call Get_PElist(saved_pelist)
current_pe = PE_here()
single_pelist(1) = current_pe
call Set_PElist(single_pelist)
do n=1, OBC%number_of_segments
segment => OBC%segment(n)
if (.not. segment%values_needed) cycle
write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n
write(suffix,"('_segment_',i3.3)") n
! needs documentation !! Yet, unsafe for now, causes grief for
! MOM_parameter_docs in circle_obcs on two processes.
! call get_param(PF, mdl, segnam, segstr, 'xyz')
! Clear out any old values
segstr = ''
call get_param(PF, mdl, segnam, segstr)
if (segstr == '') then
write(mesg,'("No OBC_SEGMENT_XXX_DATA string for OBC segment ",I3)') n
call MOM_error(FATAL, mesg)
endif
call parse_segment_manifest_str(trim(segstr), num_fields, fields)
if (num_fields == 0) then
call MOM_mesg('initialize_segment_data: num_fields = 0')
cycle ! cycle to next segment
endif
!There are OBC%num_obgc_tracers obgc tracers are there that are not listed in param file
segment%num_fields = num_fields + OBC%num_obgc_tracers
allocate(segment%field(segment%num_fields))
segment%temp_segment_data_exists = .false.
segment%salt_segment_data_exists = .false.
!!
! CODE HERE FOR OTHER OPTIONS (CLAMPED, NUDGED,..)
!!
isd = segment%HI%isd ; ied = segment%HI%ied
jsd = segment%HI%jsd ; jed = segment%HI%jed
IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB
JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB
obgc_segments_props_list => OBC%obgc_segments_props !pointer to the head node
do m=1,segment%num_fields
if (m <= num_fields) then
!These are tracers with segments specified in MOM6 style override files
call parse_segment_data_str(trim(segstr), m, trim(fields(m)), value, filename, fieldname)
else
!These are obgc tracers with segments specified by external modules.
!Set a flag so that these can be distinguished from native tracers as they may need
!extra steps for preparation and handling.
segment%field(m)%genre = 'obgc'
!Query the obgc segment properties by traversing the linkedlist
call get_obgc_segments_props(obgc_segments_props_list,fields(m),filename,fieldname,&
segment%field(m)%resrv_lfac_in,segment%field(m)%resrv_lfac_out)
!Make sure the obgc tracer is not specified in the MOM6 param file too.
do mm=1,num_fields
if (trim(fields(m)) == trim(fields(mm))) then
if (is_root_pe()) &
call MOM_error(FATAL,"MOM_open_boundary:initialize_segment_data(): obgc tracer " //trim(fields(m))// &
" appears in OBC_SEGMENT_XXX_DATA string in MOM6 param file. This is not supported!")
endif
enddo
endif
if (trim(filename) /= 'none') then
OBC%update_OBC = .true. ! Data is assumed to be time-dependent if we are reading from file
OBC%needs_IO_for_data = .true. ! At least one segment is using I/O for OBC data
! segment%values_needed = .true. ! Indicates that i/o will be needed for this segment
segment%field(m)%name = trim(fields(m))
! The scale factor for tracers may also be set in register_segment_tracer, and a constant input
! value is rescaled there.
segment%field(m)%scale = scale_factor_from_name(fields(m), GV, US, segment%tr_Reg)
segment%field(m)%use_IO = .true.
if (segment%field(m)%name == 'TEMP') then
segment%temp_segment_data_exists = .true.
segment%t_values_needed = .false.
endif
if (segment%field(m)%name == 'SALT') then
segment%salt_segment_data_exists = .true.
segment%s_values_needed = .false.
endif
filename = trim(inputdir)//trim(filename)
fieldname = trim(fieldname)//trim(suffix)
call field_size(filename,fieldname,siz,no_domain=.true.)
! if (siz(4) == 1) segment%values_needed = .false.
if (segment%on_pe) then
if (OBC%brushcutter_mode .and. (modulo(siz(1),2) == 0 .or. modulo(siz(2),2) == 0)) then
write(mesg,'("Brushcutter mode sizes ", I6, I6)') siz(1), siz(2)
call MOM_error(WARNING, mesg // " " // trim(filename) // " " // trim(fieldname))
call MOM_error(FATAL,'segment data are not on the supergrid')
endif
siz2(1) = 1
if (siz(1)>1) then
if (OBC%brushcutter_mode) then
siz2(1) = (siz(1)-1)/2
else
siz2(1) = siz(1)
endif
endif
siz2(2) = 1
if (siz(2)>1) then
if (OBC%brushcutter_mode) then
siz2(2) = (siz(2)-1)/2
else
siz2(2) = siz(2)
endif
endif
siz2(3) = siz(3)
if (segment%is_E_or_W) then
if (segment%field(m)%name == 'V') then
allocate(segment%field(m)%buffer_src(IsdB:IedB,JsdB:JedB,siz2(3)))
segment%v_values_needed = .false.
elseif (segment%field(m)%name == 'Vamp') then
allocate(segment%field(m)%buffer_src(IsdB:IedB,JsdB:JedB,siz2(3)))
segment%vamp_values_needed = .false.
segment%vamp_index = m
elseif (segment%field(m)%name == 'Vphase') then
allocate(segment%field(m)%buffer_src(IsdB:IedB,JsdB:JedB,siz2(3)))
segment%vphase_values_needed = .false.
segment%vphase_index = m
elseif (segment%field(m)%name == 'DVDX') then
allocate(segment%field(m)%buffer_src(IsdB:IedB,JsdB:JedB,siz2(3)))
segment%g_values_needed = .false.
else
allocate(segment%field(m)%buffer_src(IsdB:IedB,jsd:jed,siz2(3)))
if (segment%field(m)%name == 'U') then
segment%u_values_needed = .false.
elseif (segment%field(m)%name == 'Uamp') then
segment%uamp_values_needed = .false.
segment%uamp_index = m
elseif (segment%field(m)%name == 'Uphase') then
segment%uphase_values_needed = .false.
segment%uphase_index = m
elseif (segment%field(m)%name == 'SSH') then
segment%z_values_needed = .false.
elseif (segment%field(m)%name == 'SSHamp') then
segment%zamp_values_needed = .false.
segment%zamp_index = m
elseif (segment%field(m)%name == 'SSHphase') then
segment%zphase_values_needed = .false.
segment%zphase_index = m
elseif (segment%field(m)%name == 'TEMP') then
segment%t_values_needed = .false.
elseif (segment%field(m)%name == 'SALT') then
segment%s_values_needed = .false.
endif
endif
else
if (segment%field(m)%name == 'U') then
allocate(segment%field(m)%buffer_src(IsdB:IedB,JsdB:JedB,siz2(3)))
segment%u_values_needed = .false.
elseif (segment%field(m)%name == 'DUDY') then
allocate(segment%field(m)%buffer_src(IsdB:IedB,JsdB:JedB,siz2(3)))
segment%g_values_needed = .false.
elseif (segment%field(m)%name == 'Uamp') then
allocate(segment%field(m)%buffer_src(IsdB:IedB,JsdB:JedB,siz2(3)))
segment%uamp_values_needed = .false.
segment%uamp_index = m
elseif (segment%field(m)%name == 'Uphase') then
allocate(segment%field(m)%buffer_src(IsdB:IedB,JsdB:JedB,siz2(3)))
segment%uphase_values_needed = .false.
segment%uphase_index = m
else
allocate(segment%field(m)%buffer_src(isd:ied,JsdB:JedB,siz2(3)))
if (segment%field(m)%name == 'V') then
segment%v_values_needed = .false.
elseif (segment%field(m)%name == 'Vamp') then
segment%vamp_values_needed = .false.
segment%vamp_index = m
elseif (segment%field(m)%name == 'Vphase') then
segment%vphase_values_needed = .false.
segment%vphase_index = m
elseif (segment%field(m)%name == 'SSH') then
segment%z_values_needed = .false.
elseif (segment%field(m)%name == 'SSHamp') then
segment%zamp_values_needed = .false.
segment%zamp_index = m