forked from mom-ocean/MOM6
-
Notifications
You must be signed in to change notification settings - Fork 62
/
Copy pathMOM_state_initialization.F90
2520 lines (2214 loc) · 115 KB
/
MOM_state_initialization.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
!> Initialize state variables, u, v, h, T and S.
module MOM_state_initialization
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_checksums, only : hchksum, qchksum, uchksum, vchksum, chksum
use MOM_coms, only : max_across_PEs, min_across_PEs
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP
use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast
use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID
use MOM_EOS, only : find_depth_of_pressure_in_cell
use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : get_param, read_param, log_param, param_file_type
use MOM_file_parser, only : log_version
use MOM_get_input, only : directories
use MOM_grid, only : ocean_grid_type, isPointInCell
use MOM_interface_heights, only : find_eta
use MOM_io, only : close_file, fieldtype, file_exists
use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE
use MOM_io, only : slasher, vardesc, write_field
use MOM_io, only : EAST_FACE, NORTH_FACE
use MOM_grid_initialize, only : initialize_masks, set_grid_metrics
use MOM_restart, only : restore_state, MOM_restart_CS
use MOM_sponge, only : set_up_sponge_field, set_up_sponge_ML_density
use MOM_sponge, only : initialize_sponge, sponge_CS
use MOM_ALE_sponge, only : set_up_ALE_sponge_field, initialize_ALE_sponge
use MOM_ALE_sponge, only : ALE_sponge_CS
use MOM_string_functions, only : uppercase
use MOM_time_manager, only : time_type, set_time
use MOM_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type
use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type
use MOM_variables, only : OBC_NONE, OBC_SIMPLE, OBC_FLATHER_E, OBC_FLATHER_W
use MOM_variables, only : OBC_FLATHER_N, OBC_FLATHER_S
use MOM_verticalGrid, only : setVerticalGridAxes, verticalGrid_type
use MOM_ALE, only : pressure_gradient_plm
use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type
use MOM_EOS, only : int_specific_vol_dp
use user_initialization, only : user_initialize_thickness, user_initialize_velocity
use user_initialization, only : user_init_temperature_salinity
use user_initialization, only : user_set_Open_Bdry_Conds, user_initialize_sponges
use DOME_initialization, only : DOME_initialize_thickness
use DOME_initialization, only : DOME_set_Open_Bdry_Conds
use DOME_initialization, only : DOME_initialize_sponges
use ISOMIP_initialization, only : ISOMIP_initialize_thickness
use ISOMIP_initialization, only : ISOMIP_initialize_sponges
use ISOMIP_initialization, only : ISOMIP_initialize_temperature_salinity
use baroclinic_zone_initialization, only : baroclinic_zone_init_temperature_salinity
use benchmark_initialization, only : benchmark_initialize_thickness
use benchmark_initialization, only : benchmark_init_temperature_salinity
use circle_obcs_initialization, only : circle_obcs_initialize_thickness
use lock_exchange_initialization, only : lock_exchange_initialize_thickness
use external_gwave_initialization, only : external_gwave_initialize_thickness
use DOME2d_initialization, only : DOME2d_initialize_thickness
use DOME2d_initialization, only : DOME2d_initialize_temperature_salinity
use DOME2d_initialization, only : DOME2d_initialize_sponges
use adjustment_initialization, only : adjustment_initialize_thickness
use adjustment_initialization, only : adjustment_initialize_temperature_salinity
use sloshing_initialization, only : sloshing_initialize_thickness
use sloshing_initialization, only : sloshing_initialize_temperature_salinity
use seamount_initialization, only : seamount_initialize_thickness
use seamount_initialization, only : seamount_initialize_temperature_salinity
use Phillips_initialization, only : Phillips_initialize_thickness
use Phillips_initialization, only : Phillips_initialize_velocity
use Phillips_initialization, only : Phillips_initialize_sponges
use Rossby_front_2d_initialization, only : Rossby_front_initialize_thickness
use Rossby_front_2d_initialization, only : Rossby_front_initialize_temperature_salinity
use Rossby_front_2d_initialization, only : Rossby_front_initialize_velocity
use SCM_idealized_hurricane, only : SCM_idealized_hurricane_TS_init
use midas_vertmap, only : find_interfaces, tracer_Z_init
use midas_vertmap, only : determine_temperature
use MOM_ALE, only : ALE_initRegridding, ALE_CS, ALE_initThicknessToCoord
use MOM_ALE, only : ALE_remap_scalar, ALE_build_grid
use MOM_regridding, only : regridding_CS, set_regrid_params
use MOM_remapping, only : remapping_CS, initialize_remapping
use MOM_tracer_initialization_from_Z, only : horiz_interp_and_extrap_tracer
implicit none ; private
#include <MOM_memory.h>
public MOM_initialize_state
character(len=40) :: mod = "MOM_state_initialization" ! This module's name.
contains
! -----------------------------------------------------------------------------
subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, &
restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, &
ALE_sponge_CSp, OBC, Time_in)
type(ocean_grid_type), intent(inout) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: u
real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: v
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: h
type(thermo_var_ptrs), intent(inout) :: tv
type(time_type), intent(inout) :: Time
type(param_file_type), intent(in) :: PF
type(directories), intent(in) :: dirs
type(MOM_restart_CS), pointer :: restart_CS
type(ALE_CS), pointer :: ALE_CSp
type(tracer_registry_type), pointer :: tracer_Reg
type(sponge_CS), pointer :: sponge_CSp
type(ALE_sponge_CS), pointer :: ALE_sponge_CSp
type(ocean_OBC_type), pointer :: OBC
type(time_type), optional, intent(in) :: Time_in
! Arguments: u - Zonal velocity, in m s-1.
! (out) v - Meridional velocity, in m s-1.
! (out) h - Layer thickness, in m.
! (out) tv - A structure containing pointers to any available
! thermodynamic fields, including potential temperature and
! salinity or mixed layer density. Absent fields have NULL ptrs.
! (out) Time - Time at the start of the run segment.
! (inout) G - The ocean's grid structure.
! (in) GV - The ocean's vertical grid structure.
! (in) PF - A structure indicating the open file to parse for
! model parameter values.
! (in) dirs - A structure containing several relevant directory paths.
! (inout) restart_CS - A pointer to the restart control structure.
! (inout) CS - A structure of pointers to be exchanged with MOM.F90.
! (in) Time_in - Time at the start of the run segment. Time_in overrides
! any value set for Time.
character(len=200) :: filename ! The name of an input file.
character(len=200) :: filename2 ! The name of an input files.
character(len=200) :: inputdir ! The directory where NetCDF input files are.
character(len=200) :: config
logical :: from_Z_file, useALE
logical :: new_sim
integer :: write_geom
logical :: use_temperature, use_sponge
logical :: use_EOS ! If true, density is calculated from T & S using an
! equation of state.
logical :: depress_sfc ! If true, remove the mass that would be displaced
! by a large surface pressure by squeezing the column.
logical :: trim_ic_for_p_surf ! If true, remove the mass that would be displaced
! by a large surface pressure, such as with an ice sheet.
logical :: Analytic_FV_PGF, obsol_test
logical :: apply_OBC_u, apply_OBC_v
logical :: apply_OBC_u_flather_east, apply_OBC_u_flather_west
logical :: apply_OBC_v_flather_north, apply_OBC_v_flather_south
logical :: convert
type(EOS_type), pointer :: eos => NULL()
logical :: debug ! indicates whether to write debugging output
! This include declares and sets the variable "version".
#include "version_variable.h"
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
call callTree_enter("MOM_initialize_state(), MOM_state_initialization.F90")
call log_version(PF, mod, version)
call get_param(PF, mod, "DEBUG", debug, default=.false.)
new_sim = .false.
if ((dirs%input_filename(1:1) == 'n') .and. &
(LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true.
call get_param(PF, mod, "INPUTDIR", inputdir, &
"The directory in which input files are found.", default=".")
inputdir = slasher(inputdir)
use_temperature = ASSOCIATED(tv%T)
useALE = associated(ALE_CSp)
use_EOS = associated(tv%eqn_of_state)
if (use_EOS) eos => tv%eqn_of_state
!====================================================================
! Initialize temporally evolving fields, either as initial
! conditions or by reading them from a restart (or saves) file.
!====================================================================
if (new_sim) then
! This block initializes all of the fields internally. !
call MOM_mesg("Run initialized internally.", 3)
if (present(Time_in)) Time = Time_in
! Otherwise leave Time at its input value.
call get_param(PF, mod, "INIT_LAYERS_FROM_Z_FILE", from_Z_file, &
"If true, intialize the layer thicknesses, temperatures, \n"//&
"and salnities from a Z-space file on a latitude- \n"//&
"longitude grid.", default=.false.)
! h will be converted from m to H below
h(:,:,:) = GV%Angstrom_z
if (from_Z_file) then
! Initialize thickness and T/S from z-coordinate data in a file.
if (.NOT.use_temperature) call MOM_error(FATAL,"MOM_initialize_state : "//&
"use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")
call MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs)
call pass_var(h, G%Domain)
else
! Initialize thickness, h.
call get_param(PF, mod, "THICKNESS_CONFIG", config, &
"A string that determines how the initial layer \n"//&
"thicknesses are specified for a new run: \n"//&
" \t file - read interface heights from the file specified \n"//&
" \t thickness_file - read thicknesses from the file specified \n"//&
" \t\t by (THICKNESS_FILE).\n"//&
" \t coord - determined by ALE coordinate.\n"//&
" \t uniform - uniform thickness layers evenly distributed \n"//&
" \t\t between the surface and MAXIMUM_DEPTH. \n"//&
" \t DOME - use a slope and channel configuration for the \n"//&
" \t\t DOME sill-overflow test case. \n"//&
" \t ISOMIP - use a configuration for the \n"//&
" \t\t ISOMIP test case. \n"//&
" \t benchmark - use the benchmark test case thicknesses. \n"//&
" \t search - search a density profile for the interface \n"//&
" \t\t densities. This is not yet implemented. \n"//&
" \t circle_obcs - the circle_obcs test case is used. \n"//&
" \t DOME2D - 2D version of DOME initialization. \n"//&
" \t adjustment2d - TBD AJA. \n"//&
" \t sloshing - TBD AJA. \n"//&
" \t seamount - TBD AJA. \n"//&
" \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
" \t USER - call a user modified routine.", &
fail_if_missing=.true.)
select case (trim(config))
case ("file"); call initialize_thickness_from_file(h, G, GV, PF, .false.)
case ("thickness_file"); call initialize_thickness_from_file(h, G, GV, PF, .true.)
case ("coord")
if (useALE) then
call ALE_initThicknessToCoord( ALE_CSp, G, h )
else
call MOM_error(FATAL, "MOM_initialize_state: USE_REGRIDDING must be True "//&
"for THICKNESS_CONFIG of 'coord'")
endif
case ("uniform"); call initialize_thickness_uniform(h, G, GV, PF)
case ("DOME"); call DOME_initialize_thickness(h, G, GV, PF)
case ("ISOMIP"); call ISOMIP_initialize_thickness(h, G, GV, PF, tv)
case ("benchmark"); call benchmark_initialize_thickness(h, G, GV, PF, &
tv%eqn_of_state, tv%P_Ref)
case ("search"); call initialize_thickness_search
case ("circle_obcs"); call circle_obcs_initialize_thickness(h, G, GV, PF)
case ("lock_exchange"); call lock_exchange_initialize_thickness(h, G, GV, PF)
case ("external_gwave"); call external_gwave_initialize_thickness(h, G, PF)
case ("DOME2D"); call DOME2d_initialize_thickness(h, G, GV, PF)
case ("adjustment2d"); call adjustment_initialize_thickness(h, G, GV, PF)
case ("sloshing"); call sloshing_initialize_thickness(h, G, GV, PF)
case ("seamount"); call seamount_initialize_thickness(h, G, GV, PF)
case ("phillips"); call Phillips_initialize_thickness(h, G, GV, PF)
case ("rossby_front"); call Rossby_front_initialize_thickness(h, G, GV, PF)
case ("USER"); call user_initialize_thickness(h, G, PF, tv%T)
case default ; call MOM_error(FATAL, "MOM_initialize_state: "//&
"Unrecognized layer thickness configuration "//trim(config))
end select
call pass_var(h, G%Domain)
! Initialize temperature and salinity (T and S).
if ( use_temperature ) then
call get_param(PF, mod, "TS_CONFIG", config, &
"A string that determines how the initial tempertures \n"//&
"and salinities are specified for a new run: \n"//&
" \t file - read velocities from the file specified \n"//&
" \t\t by (TS_FILE). \n"//&
" \t fit - find the temperatures that are consistent with \n"//&
" \t\t the layer densities and salinity S_REF. \n"//&
" \t TS_profile - use temperature and salinity profiles \n"//&
" \t\t (read from TS_FILE) to set layer densities. \n"//&
" \t benchmark - use the benchmark test case T & S. \n"//&
" \t linear - linear in logical layer space. \n"//&
" \t DOME2D - 2D DOME initialization. \n"//&
" \t ISOMIP - ISOMIP initialization. \n"//&
" \t adjustment2d - TBD AJA. \n"//&
" \t sloshing - TBD AJA. \n"//&
" \t seamount - TBD AJA. \n"//&
" \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
" \t SCM_ideal_hurr - used in the SCM idealized hurricane test.\n"//&
" \t USER - call a user modified routine.", &
fail_if_missing=.true.)
! " \t baroclinic_zone - an analytic baroclinic zone. \n"//&
select case (trim(config))
case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, G, GV, PF, eos, tv%P_Ref)
case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, G, PF)
case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, &
G, GV, PF, eos, tv%P_Ref)
case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, G, PF)
case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, G, PF)
case ("DOME2D"); call DOME2d_initialize_temperature_salinity ( tv%T, &
tv%S, h, G, PF, eos)
case ("ISOMIP"); call ISOMIP_initialize_temperature_salinity ( tv%T, &
tv%S, h, G, GV, PF, eos)
case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, &
tv%S, h, G, PF, eos)
case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, &
tv%S, h, G, PF)
case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, &
tv%S, h, G, PF, eos)
case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, &
tv%S, h, G, GV, PF, eos)
case ("rossby_front"); call Rossby_front_initialize_temperature_salinity ( tv%T, &
tv%S, h, G, PF, eos)
case ("SCM_ideal_hurr"); call SCM_idealized_hurricane_TS_init ( tv%T, &
tv%S, h, G, GV, PF)
case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, G, PF, eos)
case default ; call MOM_error(FATAL, "MOM_initialize_state: "//&
"Unrecognized Temp & salt configuration "//trim(config))
end select
endif
endif ! not from_Z_file.
! Initialize velocity components, u and v
call get_param(PF, mod, "VELOCITY_CONFIG", config, &
"A string that determines how the initial velocities \n"//&
"are specified for a new run: \n"//&
" \t file - read velocities from the file specified \n"//&
" \t\t by (VELOCITY_FILE). \n"//&
" \t zero - the fluid is initially at rest. \n"//&
" \t uniform - the flow is uniform (determined by\n"//&
" \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//&
" \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
" \t USER - call a user modified routine.", default="zero")
select case (trim(config))
case ("file"); call initialize_velocity_from_file(u, v, G, PF)
case ("zero"); call initialize_velocity_zero(u, v, G, PF)
case ("uniform"); call initialize_velocity_uniform(u, v, G, PF)
case ("circular"); call initialize_velocity_circular(u, v, G, PF)
case ("phillips"); call Phillips_initialize_velocity(u, v, G, GV, PF)
case ("rossby_front"); call Rossby_front_initialize_velocity(u, v, h, G, GV, PF)
case ("USER"); call user_initialize_velocity(u, v, G, PF)
case default ; call MOM_error(FATAL, "MOM_initialize_state: "//&
"Unrecognized velocity configuration "//trim(config))
end select
call pass_vector(u, v, G%Domain)
if (debug) call uchksum(u, "MOM_initialize_state: u ", G, haloshift=1)
if (debug) call vchksum(v, "MOM_initialize_state: v ", G, haloshift=1)
! Optionally convert the thicknesses from m to kg m-2. This is particularly
! useful in a non-Boussinesq model.
call get_param(PF, mod, "CONVERT_THICKNESS_UNITS", convert, &
"If true, convert the thickness initial conditions from \n"//&
"units of m to kg m-2 or vice versa, depending on whether \n"//&
"BOUSSINESQ is defined. This does not apply if a restart \n"//&
"file is read.", default=.false.)
if (convert .and. .not. GV%Boussinesq) then
! Convert h from m to kg m-2 then to thickness units (H)
call convert_thickness(h, G, GV, PF, tv)
elseif (GV%Boussinesq) then
! Convert h from m to thickness units (H)
h(:,:,:) = h(:,:,:)*GV%m_to_H
else
h(:,:,:) = h(:,:,:)*GV%kg_m2_to_H
endif
! Remove the mass that would be displaced by an ice shelf or inverse barometer.
call get_param(PF, mod, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
"If true, depress the initial surface to avoid huge \n"//&
"tsunamis when a large surface pressure is applied.", &
default=.false.)
call get_param(PF, mod, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
"If true, cuts way the top of the column for initial conditions\n"//&
"at the depth where the hydrostatic presure matches the imposed\n"//&
"surface pressure which is read from file.", default=.false.)
if (depress_sfc .and. trim_ic_for_p_surf) call MOM_error(FATAL, "MOM_initialize_state: "//&
"DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
if (depress_sfc) call depress_surface(h, G, GV, PF, tv)
if (trim_ic_for_p_surf) call trim_for_ice(PF, G, GV, ALE_CSp, tv, h)
else ! Previous block for new_sim=.T., this block restores state
! This line calls a subroutine that reads the initial conditions !
! from a previously generated file. !
call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, &
G, restart_CS)
if (present(Time_in)) Time = Time_in
endif
if ( use_temperature ) then
call pass_var(tv%T, G%Domain, complete=.false.)
call pass_var(tv%S, G%Domain, complete=.false.)
endif
call pass_var(h, G%Domain)
if (debug) then
call hchksum(h*GV%H_to_m, "MOM_initialize_state: h ", G, haloshift=1)
if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", G, haloshift=1)
if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", G, haloshift=1)
endif
call get_param(PF, mod, "SPONGE", use_sponge, &
"If true, sponges may be applied anywhere in the domain. \n"//&
"The exact location and properties of those sponges are \n"//&
"specified via SPONGE_CONFIG.", default=.false.)
if ( use_sponge ) then
! The 3 arguments here are (1) a flag indicating whether the sponge !
! values are to be read from files, (2) the name of a file containing!
! the state toward which the model is damped, and (3) the file in !
! which the 2-D damping rate field can be found. !
call get_param(PF, mod, "SPONGE_CONFIG", config, &
"A string that sets how the sponges are configured: \n"//&
" \t file - read sponge properties from the file \n"//&
" \t\t specified by (SPONGE_FILE).\n"//&
" \t ISOMIP - apply ale sponge in the ISOMIP case \n"//&
" \t DOME - use a slope and channel configuration for the \n"//&
" \t\t DOME sill-overflow test case. \n"//&
" \t USER - call a user modified routine.", default="file")
select case (trim(config))
case ("DOME"); call DOME_initialize_sponges(G, GV, tv, PF, sponge_CSp)
case ("DOME2D"); call DOME2d_initialize_sponges(G, GV, tv, PF, useALE, &
sponge_CSp, ALE_sponge_CSp)
case ("ISOMIP"); call ISOMIP_initialize_sponges(G, GV, tv, PF, ALE_sponge_CSp)
case ("USER"); call user_initialize_sponges(G, use_temperature, tv, &
PF, sponge_CSp, h)
case ("phillips"); call Phillips_initialize_sponges(G, use_temperature, tv, &
PF, sponge_CSp, h)
case ("file"); call initialize_sponges_file(G, GV, use_temperature, tv, &
PF, sponge_CSp)
case default ; call MOM_error(FATAL, "MOM_initialize_state: "//&
"Unrecognized sponge configuration "//trim(config))
end select
endif
! This subroutine call sets optional open boundary conditions.
call get_param(PF, mod, "APPLY_OBC_U", apply_OBC_u, &
"If true, open boundary conditions may be set at some \n"//&
"u-points, with the configuration controlled by OBC_CONFIG", &
default=.false.)
call get_param(PF, mod, "APPLY_OBC_V", apply_OBC_v, &
"If true, open boundary conditions may be set at some \n"//&
"v-points, with the configuration controlled by OBC_CONFIG", &
default=.false.)
if (apply_OBC_u .or. apply_OBC_v) then
call get_param(PF, mod, "OBC_CONFIG", config, &
"A string that sets how the open boundary conditions are \n"//&
" configured: \n"//&
" \t DOME - use a slope and channel configuration for the \n"//&
" \t\t DOME sill-overflow test case. \n"//&
" \t USER - call a user modified routine.", default="file", &
fail_if_missing=.true.)
if (trim(config) == "DOME") then
call DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, PF, tracer_Reg)
elseif (trim(config) == "USER") then
call user_set_Open_Bdry_Conds(OBC, tv, G, PF, tracer_Reg)
else
call MOM_error(FATAL, "The open boundary conditions specified by "//&
"OBC_CONFIG = "//trim(config)//" have not been fully implemented.")
call set_Open_Bdry_Conds(OBC, tv, G, GV, PF, tracer_Reg)
endif
endif
call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", apply_OBC_u_flather_east,&
"Apply a Flather open boundary condition on the eastern \n"//&
"side of the global domain", default=.false.)
call get_param(PF, mod, "APPLY_OBC_U_FLATHER_WEST", apply_OBC_u_flather_west,&
"Apply a Flather open boundary condition on the western \n"//&
"side of the global domain", default=.false.)
call get_param(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", apply_OBC_v_flather_north,&
"Apply a Flather open boundary condition on the northern \n"//&
"side of the global domain", default=.false.)
call get_param(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,&
"Apply a Flather open boundary condition on the southern \n"//&
"side of the global domain", default=.false.)
if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west .or. apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) then
call set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg)
endif
call callTree_leave('MOM_initialize_state()')
end subroutine MOM_initialize_state
! -----------------------------------------------------------------------------
! -----------------------------------------------------------------------------
subroutine initialize_thickness_from_file(h, G, GV, param_file, file_has_thickness)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h
type(param_file_type), intent(in) :: param_file
logical, intent(in) :: file_has_thickness
! Arguments: h - The thickness that is being initialized.
! (in) G - The ocean's grid structure.
! (in) GV - The ocean's vertical grid structure.
! (in) param_file - A structure indicating the open file to parse for
! model parameter values.
! (in) file_has_thickness - If true, this file contains thicknesses;
! otherwise it contains interface heights.
! This subroutine reads the layer thicknesses from file.
real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1)
integer :: inconsistent = 0
real :: dilate ! The amount by which each layer is dilated to agree
! with the bottom depth and free surface height, nondim.
logical :: correct_thickness
character(len=40) :: mod = "initialize_thickness_from_file" ! This subroutine's name.
character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path
integer :: i, j, k, is, ie, js, je, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90")
call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
call get_param(param_file, mod, "THICKNESS_FILE", thickness_file, &
"The name of the thickness file.", fail_if_missing=.true.)
filename = trim(inputdir)//trim(thickness_file)
call log_param(param_file, mod, "INPUTDIR/THICKNESS_FILE", filename)
if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
" initialize_thickness_from_file: Unable to open "//trim(filename))
if (file_has_thickness) then
call read_data(filename,"h",h(:,:,:),domain=G%Domain%mpp_domain)
else
call get_param(param_file, mod, "ADJUST_THICKNESS", correct_thickness, &
"If true, all mass below the bottom removed if the \n"//&
"topography is shallower than the thickness input file \n"//&
"would indicate.", default=.false.)
call read_data(filename,"eta",eta(:,:,:),domain=G%Domain%mpp_domain)
if (correct_thickness) then
call adjustEtaToFitBathymetry(G, GV, eta, h)
else
do k=nz,1,-1 ; do j=js,je ; do i=is,ie
if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_z)) then
eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_z
h(i,j,k) = GV%Angstrom_z
else
h(i,j,k) = eta(i,j,K) - eta(i,j,K+1)
endif
enddo ; enddo ; enddo
do j=js,je ; do i=is,ie
if (abs(eta(i,j,nz+1) + G%bathyT(i,j)) > 1.0) &
inconsistent = inconsistent + 1
enddo ; enddo
call sum_across_PEs(inconsistent)
if ((inconsistent > 0) .and. (is_root_pe())) then
write(mesg,'("Thickness initial conditions are inconsistent ",'// &
'"with topography in ",I8," places.")') inconsistent
call MOM_error(WARNING, mesg)
endif
endif
endif
call callTree_leave(trim(mod)//'()')
end subroutine initialize_thickness_from_file
! -----------------------------------------------------------------------------
! -----------------------------------------------------------------------------
!> Adjust interface heights to fit the bathymetry and diagnose layer thickness.
!! If the bottom most interface is below the topography then the bottom-most
!! layers are contracted to GV%Angstrom_z.
!! If the bottom most interface is above the topography then the entire column
!! is dilated (expanded) to fill the void.
!! @remark{There is a (hard-wired) "tolerance" parameter such that the
!! criteria for adjustment must equal or exceed 10cm.}
!! @param[in] G Grid type
!! @param[in,out] eta Interface heights
!! @param[out] h Layer thicknesses
subroutine adjustEtaToFitBathymetry(G, GV, eta, h)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G), SZK_(G)+1), intent(inout) :: eta
real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h
! Local variables
integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
real, parameter :: hTolerance = 0.1 !< Tolerance to exceed adjustment criteria (m)
real :: hTmp, eTmp, dilate
character(len=100) :: mesg
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
contractions = 0
do j=js,je ; do i=is,ie
if (-eta(i,j,nz+1) > G%bathyT(i,j) + hTolerance) then
eta(i,j,nz+1) = -G%bathyT(i,j)
contractions = contractions + 1
endif
enddo ; enddo
call sum_across_PEs(contractions)
if ((contractions > 0) .and. (is_root_pe())) then
write(mesg,'("Thickness initial conditions were contracted ",'// &
'"to fit topography in ",I8," places.")') contractions
call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
endif
do k=nz,1,-1 ; do j=js,je ; do i=is,ie
! Collapse layers to thinnest possible if the thickness less than
! the thinnest possible (or negative).
if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_z)) then
eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_z
h(i,j,k) = GV%Angstrom_z
else
h(i,j,k) = eta(i,j,K) - eta(i,j,K+1)
endif
enddo ; enddo ; enddo
dilations = 0
do j=js,je ; do i=is,ie
! The whole column is dilated to accommodate deeper topography than
! the bathymetry would indicate.
! This should be... if ((G%mask2dt(i,j)*(eta(i,j,1)-eta(i,j,nz+1)) > 0.0) .and. &
if (-eta(i,j,nz+1) < G%bathyT(i,j) - hTolerance) then
dilations = dilations + 1
if (eta(i,j,1) <= eta(i,j,nz+1)) then
do k=1,nz ; h(i,j,k) = (eta(i,j,1)+G%bathyT(i,j)) / real(nz) ; enddo
else
dilate = (eta(i,j,1)+G%bathyT(i,j)) / (eta(i,j,1)-eta(i,j,nz+1))
do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
endif
do k=nz, 2, -1; eta(i,j,K) = eta(i,j,K+1) + h(i,j,k); enddo
endif
enddo ; enddo
call sum_across_PEs(dilations)
if ((dilations > 0) .and. (is_root_pe())) then
write(mesg,'("Thickness initial conditions were dilated ",'// &
'"to fit topography in ",I8," places.")') dilations
call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
endif
end subroutine adjustEtaToFitBathymetry
! -----------------------------------------------------------------------------
! -----------------------------------------------------------------------------
subroutine initialize_thickness_uniform(h, G, GV, param_file)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h
type(param_file_type), intent(in) :: param_file
! Arguments: h - The thickness that is being initialized.
! (in) G - The ocean's grid structure.
! (in) GV - The ocean's vertical grid structure.
! (in) param_file - A structure indicating the open file to parse for
! model parameter values.
! This subroutine initializes the layer thicknesses to be uniform.
character(len=40) :: mod = "initialize_thickness_uniform" ! This subroutine's name.
real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually !
! negative because it is positive upward. !
real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface !
! positive upward, in m. !
integer :: i, j, k, is, ie, js, je, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90")
if (G%max_depth<=0.) call MOM_error(FATAL,"initialize_thickness_uniform: "// &
"MAXIMUM_DEPTH has a non-sensical value! Was it set?")
do k=1,nz
e0(K) = -G%max_depth * real(k-1) / real(nz)
enddo
do j=js,je ; do i=is,ie
! This sets the initial thickness (in m) of the layers. The !
! thicknesses are set to insure that: 1. each layer is at least an !
! Angstrom thick, and 2. the interfaces are where they should be !
! based on the resting depths and interface height perturbations, !
! as long at this doesn't interfere with 1. !
eta1D(nz+1) = -1.0*G%bathyT(i,j)
do k=nz,1,-1
eta1D(K) = e0(K)
if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_z)) then
eta1D(K) = eta1D(K+1) + GV%Angstrom_z
h(i,j,k) = GV%Angstrom_z
else
h(i,j,k) = eta1D(K) - eta1D(K+1)
endif
enddo
enddo ; enddo
call callTree_leave(trim(mod)//'()')
end subroutine initialize_thickness_uniform
! -----------------------------------------------------------------------------
! -----------------------------------------------------------------------------
subroutine initialize_thickness_search
! search density space for location of layers
call MOM_error(FATAL," MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")
end subroutine initialize_thickness_search
! -----------------------------------------------------------------------------
subroutine convert_thickness(h, G, GV, param_file, tv)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h
type(param_file_type), intent(in) :: param_file
type(thermo_var_ptrs), intent(in) :: tv
! Arguments: h - The thickness that is being initialized.
! (in) G - The ocean's grid structure.
! (in) GV - The ocean's vertical grid structure.
! (in) param_file - A structure indicating the open file to parse for
! model parameter values.
real, dimension(SZI_(G),SZJ_(G)) :: &
p_top, p_bot
real :: dz_geo(SZI_(G),SZJ_(G)) ! The change in geopotential height
! across a layer, in m2 s-2.
real :: rho(SZI_(G))
real :: I_gEarth
logical :: Boussinesq
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
integer :: itt, max_itt
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
max_itt = 10
Boussinesq = GV%Boussinesq
I_gEarth = 1.0 / G%g_Earth
if (Boussinesq) then
call MOM_error(FATAL,"Not yet converting thickness with Boussinesq approx.")
else
if (associated(tv%eqn_of_state)) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
enddo ; enddo
do k=1,nz
do j=js,je
do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, &
is, ie-is+1, tv%eqn_of_state)
do i=is,ie
p_bot(i,j) = p_top(i,j) + G%g_Earth * h(i,j,k) * rho(i)
enddo
enddo
do itt=1,max_itt
call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, &
0.0, G%HI, tv%eqn_of_state, dz_geo)
if (itt < max_itt) then ; do j=js,je
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_bot(:,j), rho, &
is, ie-is+1, tv%eqn_of_state)
! Use Newton's method to correct the bottom value.
! The hydrostatic equation is linear to such a
! high degree that no bounds-checking is needed.
do i=is,ie
p_bot(i,j) = p_bot(i,j) + rho(i) * (G%g_Earth*h(i,j,k) - dz_geo(i,j))
enddo
enddo ; endif
enddo
do j=js,je ; do i=is,ie
h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * GV%kg_m2_to_H * I_gEarth
enddo ; enddo
enddo
else
do k=1,nz ; do j=js,je ; do i=is,ie
h(i,j,k) = h(i,j,k) * GV%Rlay(k) * GV%kg_m2_to_H
enddo ; enddo ; enddo
endif
endif
end subroutine convert_thickness
subroutine depress_surface(h, G, GV, param_file, tv)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h
type(param_file_type), intent(in) :: param_file
type(thermo_var_ptrs), intent(in) :: tv
! Arguments: h - The thickness that is being initialized.
! (in) G - The ocean's grid structure.
! (in) GV - The ocean's vertical grid structure.
! (in) param_file - A structure indicating the open file to parse for
! model parameter values.
real, dimension(SZI_(G),SZJ_(G)) :: &
eta_sfc ! The free surface height that the model should use, in m.
real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
eta ! The free surface height that the model should use, in m.
real :: dilate ! A ratio by which layers are dilated, nondim.
real :: scale_factor ! A scaling factor for the eta_sfc values that are read
! in, which can be used to change units, for example.
character(len=40) :: mod = "depress_surface" ! This subroutine's name.
character(len=200) :: inputdir, eta_srf_file ! Strings for file/path
character(len=200) :: filename, eta_srf_var ! Strings for file/path
integer :: i, j, k, is, ie, js, je, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
! Read the surface height (or pressure) from a file.
call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
call get_param(param_file, mod, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
"The initial condition file for the surface height.", &
fail_if_missing=.true.)
call get_param(param_file, mod, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
"The initial condition variable for the surface height.",&
default="SSH")
filename = trim(inputdir)//trim(eta_srf_file)
call log_param(param_file, mod, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
call read_data(filename,eta_srf_var,eta_sfc,domain=G%Domain%mpp_domain)
call get_param(param_file, mod, "SURFACE_HEIGHT_IC_SCALE", scale_factor, &
"A scaling factor to convert SURFACE_HEIGHT_IC_VAR into \n"//&
"units of m", units="variable", default=1.0)
if (scale_factor /= 1.0) then ; do j=js,je ; do i=is,ie
eta_sfc(i,j) = eta_sfc(i,j) * scale_factor
enddo ; enddo ; endif
! Convert thicknesses to interface heights.
call find_eta(h, tv, G%g_Earth, G, GV, eta)
do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
! if (eta_sfc(i,j) < eta(i,j,nz+1)) then
! Issue a warning?
! endif
if (eta_sfc(i,j) > eta(i,j,1)) then
! Dilate the water column to agree, but only up to 10-fold.
if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1))) then
dilate = 10.0
call MOM_error(WARNING, "Free surface height dilation attempted "//&
"to exceed 10-fold.", all_print=.true.)
else
dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
endif
do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
elseif (eta(i,j,1) > eta_sfc(i,j)) then
! Remove any mass that is above the target free surface.
do k=1,nz
if (eta(i,j,K) <= eta_sfc(i,j)) exit
if (eta(i,j,K+1) >= eta_sfc(i,j)) then
h(i,j,k) = GV%Angstrom
else
h(i,j,k) = max(GV%Angstrom, h(i,j,k) * &
(eta_sfc(i,j) - eta(i,j,K+1)) / (eta(i,j,K) - eta(i,j,K+1)) )
endif
enddo
endif
endif ; enddo ; enddo
end subroutine depress_surface
!> Adjust the layer thicknesses by cutting away the top of each model column at the depth
!! where the hydrostatic pressure matches an imposed surface pressure read from file.
subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h)
type(param_file_type), intent(in) :: PF !< Parameter file structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness (H units, m or Pa)
! Local variables
character(len=200) :: mod = "trim_for_ice"
real, dimension(SZI_(G),SZJ_(G)) :: p_surf ! Imposed pressure on ocean at surface (Pa)
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_t, S_b, T_t, T_b ! Top and bottom edge values for reconstructions
! of salinity and temperature within each layer.
character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path
real :: scale_factor, min_thickness
integer :: i, j, k
call get_param(PF, mod, "SURFACE_PRESSURE_FILE", p_surf_file, &
"The initial condition file for the surface height.", &
fail_if_missing=.true.)
call get_param(PF, mod, "SURFACE_PRESSURE_VAR", p_surf_var, &
"The initial condition variable for the surface height.", &
units="kg m-2", default="")
call get_param(PF, mod, "INPUTDIR", inputdir, default=".", do_not_log=.true.)
filename = trim(slasher(inputdir))//trim(p_surf_file)
call log_param(PF, mod, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
call read_data(filename, p_surf_var, p_surf, domain=G%Domain%mpp_domain)
call get_param(PF, mod, "SURFACE_PRESSURE_SCALE", scale_factor, &
"A scaling factor to convert SURFACE_PRESSURE_VAR from\n"//&
"file SURFACE_PRESSURE_FILE into a surface pressure.", &
units="file dependent", default=1.)
if (scale_factor /= 1.) p_surf(:,:) = scale_factor * p_surf(:,:)
call get_param(PF, mod, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
units='m', default=1.e-3)
! Find edge values of T and S used in reconstructions
if ( associated(ALE_CSp) ) then ! This should only be associated if we are in ALE mode
! if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then
call pressure_gradient_plm(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h)
! elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then
! call pressure_gradient_ppm(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h)
! endif
else
! call MOM_error(FATAL, "trim_for_ice: Does not work without ALE mode")
do k=1,G%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
T_t(i,j,k) = tv%T(i,j,k) ; T_b(i,j,k) = tv%T(i,j,k)
S_t(i,j,k) = tv%S(i,j,k) ; S_b(i,j,k) = tv%S(i,j,k)
enddo ; enddo ; enddo
endif
do j=G%jsc,G%jec ; do i=G%isc,G%iec
call cut_off_column_top(GV%ke, tv, GV%Rho0, G%G_earth, G%bathyT(i,j), min_thickness, &
tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), p_surf(i,j), h(i,j,:))
enddo ; enddo
end subroutine trim_for_ice
!> Adjust the layer thicknesses by cutting away the top at the depth where the hydrostatic
!! pressure matches p_surf
subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, T, T_t, T_b, S, S_t, S_b, p_surf, h)
integer, intent(in) :: nk !< Number of layers
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
real, intent(in) :: Rho0 !< Reference density (kg/m3)
real, intent(in) :: G_earth !< Gravitational acceleration (m/s2)
real, intent(in) :: depth !< Depth of ocean column (m)
real, intent(in) :: min_thickness !< Smallest thickness allowed (m)
real, dimension(nk), intent(in) :: T !<
real, dimension(nk), intent(in) :: T_t !<
real, dimension(nk), intent(in) :: T_b !<
real, dimension(nk), intent(in) :: S !<
real, dimension(nk), intent(in) :: S_t !<
real, dimension(nk), intent(in) :: S_b !<
real, intent(in) :: p_surf !< Imposed pressure on ocean at surface (Pa)
real, dimension(nk), intent(inout) :: h !< Layer thickness (H units, m or Pa)
! Local variables
real, dimension(nk+1) :: e ! Top and bottom edge values for reconstructions
real :: P_t, P_b, z_out, e_top
integer :: k
! Calculate original interface positions
e(nk+1) = -depth
do k=nk,1,-1
e(K) = e(K+1) + h(k)
enddo
P_t = 0.
e_top = e(1)
do k=1,nk
call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), &
P_t, p_surf, Rho0, G_earth, tv%eqn_of_state, P_b, z_out)
if (z_out>=e(K)) then
! Imposed pressure was less that pressure at top of cell
exit
elseif (z_out<=e(K+1)) then
! Imposed pressure was greater than pressure at bottom of cell
e_top = e(K+1)
else
! Imposed pressure was fell between pressures at top and bottom of cell
e_top = z_out
exit
endif
P_t = P_b
enddo
if (e_top<e(1)) then
! Clip layers from the top down, if at all
do K=1,nk
if (e(K)>e_top) then
! Original e(K) is too high
e(K) = e_top
e_top = e_top - min_thickness ! Next interface must be at least this deep
endif
! This layer needs trimming
h(k) = max( min_thickness, e(K) - e(K+1) )
if (e(K)<e_top) exit ! No need to go further
enddo
endif
end subroutine cut_off_column_top
! -----------------------------------------------------------------------------
subroutine initialize_velocity_from_file(u, v, G, param_file)
type(ocean_grid_type), intent(in) :: G
real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u
real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v
type(param_file_type), intent(in) :: param_file
! Arguments: u - The zonal velocity that is being initialized.
! (out) v - The meridional velocity that is being initialized.
! (in) G - The ocean's grid structure.
! (in) param_file - parameter file type
! This subroutine reads the initial velocity components from file
character(len=40) :: mod = "initialize_velocity_from_file" ! This subroutine's name.
character(len=200) :: filename,velocity_file,inputdir ! Strings for file/path
call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90")
call get_param(param_file, mod, "VELOCITY_FILE", velocity_file, &
"The name of the velocity initial condition file.", &
fail_if_missing=.true.)
call get_param(param_file, mod, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
filename = trim(inputdir)//trim(velocity_file)
call log_param(param_file, mod, "INPUTDIR/VELOCITY_FILE", filename)
if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
" initialize_velocity_from_file: Unable to open "//trim(filename))
! Read the velocities from a netcdf file.
call read_data(filename,"u",u(:,:,:),domain=G%Domain%mpp_domain,position=EAST_FACE)
call read_data(filename,"v",v(:,:,:),domain=G%Domain%mpp_domain,position=NORTH_FACE)
call callTree_leave(trim(mod)//'()')
end subroutine initialize_velocity_from_file
! -----------------------------------------------------------------------------
! -----------------------------------------------------------------------------
subroutine initialize_velocity_zero(u, v, G, param_file)