-
Notifications
You must be signed in to change notification settings - Fork 246
/
Copy pathMOM_barotropic.F90
5182 lines (4805 loc) · 266 KB
/
MOM_barotropic.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
!> Baropotric solver
module MOM_barotropic
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_debugging, only : hchksum, uvchksum
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field
use MOM_diag_mediator, only : safe_alloc_ptr, diag_ctrl, enable_averaging
use MOM_domains, only : min_across_PEs, clone_MOM_domain, deallocate_MOM_domain
use MOM_domains, only : To_All, Scalar_Pair, AGRID, CORNER, MOM_domain_type
use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
use MOM_domains, only : start_group_pass, complete_group_pass, pass_var, pass_vector
use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_forcing_type, only : mech_forcing
use MOM_grid, only : ocean_grid_type
use MOM_hor_index, only : hor_index_type
use MOM_io, only : vardesc, var_desc, MOM_read_data, slasher
use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE, open_boundary_query
use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_segment_type
use MOM_restart, only : register_restart_field, register_restart_pair
use MOM_restart, only : query_initialized, MOM_restart_CS
use MOM_tidal_forcing, only : tidal_forcing_sensitivity, tidal_forcing_CS
use MOM_time_manager, only : time_type, real_to_time, operator(+), operator(-)
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : BT_cont_type, alloc_bt_cont_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_variables, only : accel_diag_ptrs
implicit none ; private
#include <MOM_memory.h>
#ifdef STATIC_MEMORY_
# ifndef BTHALO_
# define BTHALO_ 0
# endif
# define WHALOI_ MAX(BTHALO_-NIHALO_,0)
# define WHALOJ_ MAX(BTHALO_-NJHALO_,0)
# define NIMEMW_ 1-WHALOI_:NIMEM_+WHALOI_
# define NJMEMW_ 1-WHALOJ_:NJMEM_+WHALOJ_
# define NIMEMBW_ -WHALOI_:NIMEM_+WHALOI_
# define NJMEMBW_ -WHALOJ_:NJMEM_+WHALOJ_
# define SZIW_(G) NIMEMW_
# define SZJW_(G) NJMEMW_
# define SZIBW_(G) NIMEMBW_
# define SZJBW_(G) NJMEMBW_
#else
# define NIMEMW_ :
# define NJMEMW_ :
# define NIMEMBW_ :
# define NJMEMBW_ :
# define SZIW_(G) G%isdw:G%iedw
# define SZJW_(G) G%jsdw:G%jedw
# define SZIBW_(G) G%isdw-1:G%iedw
# define SZJBW_(G) G%jsdw-1:G%jedw
#endif
public btcalc, bt_mass_source, btstep, barotropic_init, barotropic_end
public register_barotropic_restarts, set_dtbt, barotropic_get_tav
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.
!> The barotropic stepping open boundary condition type
type, private :: BT_OBC_type
real, dimension(:,:), pointer :: Cg_u => NULL() !< The external wave speed at u-points [L T-1 ~> m s-1].
real, dimension(:,:), pointer :: Cg_v => NULL() !< The external wave speed at u-points [L T-1 ~> m s-1].
real, dimension(:,:), pointer :: H_u => NULL() !< The total thickness at the u-points [H ~> m or kg m-2].
real, dimension(:,:), pointer :: H_v => NULL() !< The total thickness at the v-points [H ~> m or kg m-2].
real, dimension(:,:), pointer :: uhbt => NULL() !< The zonal barotropic thickness fluxes specified
!! for open boundary conditions (if any) [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(:,:), pointer :: vhbt => NULL() !< The meridional barotropic thickness fluxes specified
!! for open boundary conditions (if any) [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(:,:), pointer :: ubt_outer => NULL() !< The zonal velocities just outside the domain,
!! as set by the open boundary conditions [L T-1 ~> m s-1].
real, dimension(:,:), pointer :: vbt_outer => NULL() !< The meridional velocities just outside the domain,
!! as set by the open boundary conditions [L T-1 ~> m s-1].
real, dimension(:,:), pointer :: eta_outer_u => NULL() !< The surface height outside of the domain
!! at a u-point with an open boundary condition [H ~> m or kg m-2].
real, dimension(:,:), pointer :: eta_outer_v => NULL() !< The surface height outside of the domain
!! at a v-point with an open boundary condition [H ~> m or kg m-2].
logical :: apply_u_OBCs !< True if this PE has an open boundary at a u-point.
logical :: apply_v_OBCs !< True if this PE has an open boundary at a v-point.
!>@{ Index ranges for the open boundary conditions
integer :: is_u_obc, ie_u_obc, js_u_obc, je_u_obc
integer :: is_v_obc, ie_v_obc, js_v_obc, je_v_obc
!>@}
logical :: is_alloced = .false. !< True if BT_OBC is in use and has been allocated
type(group_pass_type) :: pass_uv !< Structure for group halo pass
type(group_pass_type) :: pass_uhvh !< Structure for group halo pass
type(group_pass_type) :: pass_h !< Structure for group halo pass
type(group_pass_type) :: pass_cg !< Structure for group halo pass
type(group_pass_type) :: pass_eta_outer !< Structure for group halo pass
end type BT_OBC_type
!> The barotropic stepping control stucture
type, public :: barotropic_CS ; private
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: frhatu
!< The fraction of the total column thickness interpolated to u grid points in each layer [nondim].
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: frhatv
!< The fraction of the total column thickness interpolated to v grid points in each layer [nondim].
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: IDatu
!< Inverse of the basin depth at u grid points [Z-1 ~> m-1].
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: lin_drag_u
!< A spatially varying linear drag coefficient acting on the zonal barotropic flow
!! [H T-1 ~> m s-1 or kg m-2 s-1].
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: ubt_IC
!< The barotropic solvers estimate of the zonal velocity that will be the initial
!! condition for the next call to btstep [L T-1 ~> m s-1].
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: ubtav
!< The barotropic zonal velocity averaged over the baroclinic time step [L T-1 ~> m s-1].
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: IDatv
!< Inverse of the basin depth at v grid points [Z-1 ~> m-1].
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: lin_drag_v
!< A spatially varying linear drag coefficient acting on the zonal barotropic flow
!! [H T-1 ~> m s-1 or kg m-2 s-1].
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: vbt_IC
!< The barotropic solvers estimate of the zonal velocity that will be the initial
!! condition for the next call to btstep [L T-1 ~> m s-1].
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: vbtav
!< The barotropic meridional velocity averaged over the baroclinic time step [L T-1 ~> m s-1].
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: eta_cor
!< The difference between the free surface height from the barotropic calculation and the sum
!! of the layer thicknesses. This difference is imposed as a forcing term in the barotropic
!! calculation over a baroclinic timestep [H ~> m or kg m-2].
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: eta_cor_bound
!< A limit on the rate at which eta_cor can be applied while avoiding instability
!! [H T-1 ~> m s-1 or kg m-2 s-1]. This is only used if CS%bound_BT_corr is true.
real ALLOCABLE_, dimension(NIMEMW_,NJMEMW_) :: &
ua_polarity, & !< Test vector components for checking grid polarity.
va_polarity, & !< Test vector components for checking grid polarity.
bathyT !< A copy of bathyT (ocean bottom depth) with wide halos [Z ~> m]
real ALLOCABLE_, dimension(NIMEMW_,NJMEMW_) :: IareaT
!< This is a copy of G%IareaT with wide halos, but will
!! still utilize the macro IareaT when referenced, [L-2 ~> m-2].
real ALLOCABLE_, dimension(NIMEMBW_,NJMEMW_) :: &
D_u_Cor, & !< A simply averaged depth at u points [Z ~> m].
dy_Cu, & !< A copy of G%dy_Cu with wide halos [L ~> m].
IdxCu !< A copy of G%IdxCu with wide halos [L-1 ~> m-1].
real ALLOCABLE_, dimension(NIMEMW_,NJMEMBW_) :: &
D_v_Cor, & !< A simply averaged depth at v points [Z ~> m].
dx_Cv, & !< A copy of G%dx_Cv with wide halos [L ~> m].
IdyCv !< A copy of G%IdyCv with wide halos [L-1 ~> m-1].
real ALLOCABLE_, dimension(NIMEMBW_,NJMEMBW_) :: &
q_D !< f / D at PV points [Z-1 T-1 ~> m-1 s-1].
real, dimension(:,:,:), pointer :: frhatu1 => NULL() !< Predictor step values of frhatu stored for diagnostics.
real, dimension(:,:,:), pointer :: frhatv1 => NULL() !< Predictor step values of frhatv stored for diagnostics.
type(BT_OBC_type) :: BT_OBC !< A structure with all of this modules fields
!! for applying open boundary conditions.
real :: dtbt !< The barotropic time step [T ~> s].
real :: dtbt_fraction !< The fraction of the maximum time-step that
!! should used. The default is 0.98.
real :: dtbt_max !< The maximum stable barotropic time step [T ~> s].
real :: dt_bt_filter !< The time-scale over which the barotropic mode solutions are
!! filtered [T ~> s] if positive, or as a fraction of DT if
!! negative [nondim]. This can never be taken to be longer than 2*dt.
!! Set this to 0 to apply no filtering.
integer :: nstep_last = 0 !< The number of barotropic timesteps per baroclinic
!! time step the last time btstep was called.
real :: bebt !< A nondimensional number, from 0 to 1, that
!! determines the gravity wave time stepping scheme.
!! 0.0 gives a forward-backward scheme, while 1.0
!! give backward Euler. In practice, bebt should be
!! of order 0.2 or greater.
logical :: split !< If true, use the split time stepping scheme.
logical :: bound_BT_corr !< If true, the magnitude of the fake mass source
!! in the barotropic equation that drives the two
!! estimates of the free surface height toward each
!! other is bounded to avoid driving corrective
!! velocities that exceed MAXCFL_BT_CONT.
logical :: gradual_BT_ICs !< If true, adjust the initial conditions for the
!! barotropic solver to the values from the layered
!! solution over a whole timestep instead of
!! instantly. This is a decent approximation to the
!! inclusion of sum(u dh_dt) while also correcting
!! for truncation errors.
logical :: Sadourny !< If true, the Coriolis terms are discretized
!! with Sadourny's energy conserving scheme,
!! otherwise the Arakawa & Hsu scheme is used. If
!! the deformation radius is not resolved Sadourny's
!! scheme should probably be used.
logical :: integral_bt_cont !< If true, use the time-integrated velocity over the barotropic steps
!! to determine the integrated transports used to update the continuity
!! equation. Otherwise the transports are the sum of the transports
!! based on ]a series of instantaneous velocities and the BT_CONT_TYPE
!! for transports. This is only valid if a BT_CONT_TYPE is used.
logical :: Nonlinear_continuity !< If true, the barotropic continuity equation
!! uses the full ocean thickness for transport.
integer :: Nonlin_cont_update_period !< The number of barotropic time steps
!! between updates to the face area, or 0 only to
!! update at the start of a call to btstep. The
!! default is 1.
logical :: BT_project_velocity !< If true, step the barotropic velocity first
!! and project out the velocity tendency by 1+BEBT
!! when calculating the transport. The default
!! (false) is to use a predictor continuity step to
!! find the pressure field, and then do a corrector
!! continuity step using a weighted average of the
!! old and new velocities, with weights of (1-BEBT) and BEBT.
logical :: nonlin_stress !< If true, use the full depth of the ocean at the start of the
!! barotropic step when calculating the surface stress contribution to
!! the barotropic acclerations. Otherwise use the depth based on bathyT.
real :: BT_Coriolis_scale !< A factor by which the barotropic Coriolis acceleration anomaly
!! terms are scaled.
logical :: answers_2018 !< If true, use expressions for the barotropic solver that recover
!! the answers from the end of 2018. Otherwise, use more efficient
!! or general expressions.
logical :: dynamic_psurf !< If true, add a dynamic pressure due to a viscous
!! ice shelf, for instance.
real :: Dmin_dyn_psurf !< The minimum depth to use in limiting the size
!! of the dynamic surface pressure for stability [Z ~> m].
real :: ice_strength_length !< The length scale at which the damping rate
!! due to the ice strength should be the same as if
!! a Laplacian were applied [L ~> m].
real :: const_dyn_psurf !< The constant that scales the dynamic surface
!! pressure [nondim]. Stable values are < ~1.0.
!! The default is 0.9.
logical :: tides !< If true, apply tidal momentum forcing.
logical :: tidal_sal_bug !< If true, the tidal self-attraction and loading anomaly in the
!! barotropic solver has the wrong sign, replicating a long-standing
!! bug.
real :: G_extra !< A nondimensional factor by which gtot is enhanced.
integer :: hvel_scheme !< An integer indicating how the thicknesses at
!! velocity points are calculated. Valid values are
!! given by the parameters defined below:
!! HARMONIC, ARITHMETIC, HYBRID, and FROM_BT_CONT
logical :: strong_drag !< If true, use a stronger estimate of the retarding
!! effects of strong bottom drag.
logical :: linear_wave_drag !< If true, apply a linear drag to the barotropic
!! velocities, using rates set by lin_drag_u & _v
!! divided by the depth of the ocean.
logical :: linearized_BT_PV !< If true, the PV and interface thicknesses used
!! in the barotropic Coriolis calculation is time
!! invariant and linearized.
logical :: use_wide_halos !< If true, use wide halos and march in during the
!! barotropic time stepping for efficiency.
logical :: clip_velocity !< If true, limit any velocity components that are
!! are large enough for a CFL number to exceed
!! CFL_trunc. This should only be used as a
!! desperate debugging measure.
logical :: debug !< If true, write verbose checksums for debugging purposes.
logical :: debug_bt !< If true, write verbose checksums for debugging purposes.
real :: vel_underflow !< Velocity components smaller than vel_underflow
!! are set to 0 [L T-1 ~> m s-1].
real :: maxvel !< Velocity components greater than maxvel are
!! truncated to maxvel [L T-1 ~> m s-1].
real :: CFL_trunc !< If clip_velocity is true, velocity components will
!! be truncated when they are large enough that the
!! corresponding CFL number exceeds this value, nondim.
real :: maxCFL_BT_cont !< The maximum permitted CFL number associated with the
!! barotropic accelerations from the summed velocities
!! times the time-derivatives of thicknesses. The
!! default is 0.1, and there will probably be real
!! problems if this were set close to 1.
logical :: BT_cont_bounds !< If true, use the BT_cont_type variables to set limits
!! on the magnitude of the corrective mass fluxes.
logical :: visc_rem_u_uh0 !< If true, use the viscous remnants when estimating
!! the barotropic velocities that were used to
!! calculate uh0 and vh0. False is probably the
!! better choice.
logical :: adjust_BT_cont !< If true, adjust the curve fit to the BT_cont type
!! that is used by the barotropic solver to match the
!! transport about which the flow is being linearized.
logical :: use_old_coriolis_bracket_bug !< If True, use an order of operations
!! that is not bitwise rotationally symmetric in the
!! meridional Coriolis term of the barotropic solver.
type(time_type), pointer :: Time => NULL() !< A pointer to the ocean models clock.
type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate
!! the timing of diagnostic output.
type(MOM_domain_type), pointer :: BT_Domain => NULL() !< Barotropic MOM domain
type(hor_index_type), pointer :: debug_BT_HI => NULL() !< debugging copy of horizontal index_type
type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< Control structure for tides
logical :: module_is_initialized = .false. !< If true, module has been initialized
integer :: isdw !< The lower i-memory limit for the wide halo arrays.
integer :: iedw !< The upper i-memory limit for the wide halo arrays.
integer :: jsdw !< The lower j-memory limit for the wide halo arrays.
integer :: jedw !< The upper j-memory limit for the wide halo arrays.
type(group_pass_type) :: pass_q_DCor !< Handle for a group halo pass
type(group_pass_type) :: pass_gtot !< Handle for a group halo pass
type(group_pass_type) :: pass_tmp_uv !< Handle for a group halo pass
type(group_pass_type) :: pass_eta_bt_rem !< Handle for a group halo pass
type(group_pass_type) :: pass_force_hbt0_Cor_ref !< Handle for a group halo pass
type(group_pass_type) :: pass_Dat_uv !< Handle for a group halo pass
type(group_pass_type) :: pass_eta_ubt !< Handle for a group halo pass
type(group_pass_type) :: pass_etaav !< Handle for a group halo pass
type(group_pass_type) :: pass_ubt_Cor !< Handle for a group halo pass
type(group_pass_type) :: pass_ubta_uhbta !< Handle for a group halo pass
type(group_pass_type) :: pass_e_anom !< Handle for a group halo pass
!>@{ Diagnostic IDs
integer :: id_PFu_bt = -1, id_PFv_bt = -1, id_Coru_bt = -1, id_Corv_bt = -1
integer :: id_ubtforce = -1, id_vbtforce = -1, id_uaccel = -1, id_vaccel = -1
integer :: id_visc_rem_u = -1, id_visc_rem_v = -1, id_eta_cor = -1
integer :: id_ubt = -1, id_vbt = -1, id_eta_bt = -1, id_ubtav = -1, id_vbtav = -1
integer :: id_ubt_st = -1, id_vbt_st = -1, id_eta_st = -1
integer :: id_ubtdt = -1, id_vbtdt = -1
integer :: id_ubt_hifreq = -1, id_vbt_hifreq = -1, id_eta_hifreq = -1
integer :: id_uhbt_hifreq = -1, id_vhbt_hifreq = -1, id_eta_pred_hifreq = -1
integer :: id_gtotn = -1, id_gtots = -1, id_gtote = -1, id_gtotw = -1
integer :: id_uhbt = -1, id_frhatu = -1, id_vhbt = -1, id_frhatv = -1
integer :: id_frhatu1 = -1, id_frhatv1 = -1
integer :: id_BTC_FA_u_EE = -1, id_BTC_FA_u_E0 = -1, id_BTC_FA_u_W0 = -1, id_BTC_FA_u_WW = -1
integer :: id_BTC_ubt_EE = -1, id_BTC_ubt_WW = -1
integer :: id_BTC_FA_v_NN = -1, id_BTC_FA_v_N0 = -1, id_BTC_FA_v_S0 = -1, id_BTC_FA_v_SS = -1
integer :: id_BTC_vbt_NN = -1, id_BTC_vbt_SS = -1
integer :: id_BTC_FA_u_rat0 = -1, id_BTC_FA_v_rat0 = -1, id_BTC_FA_h_rat0 = -1
integer :: id_uhbt0 = -1, id_vhbt0 = -1
!>@}
end type barotropic_CS
!> A desciption of the functional dependence of transport at a u-point
type, private :: local_BT_cont_u_type
real :: FA_u_EE !< The effective open face area for zonal barotropic transport
!! drawing from locations far to the east [H L ~> m2 or kg m-1].
real :: FA_u_E0 !< The effective open face area for zonal barotropic transport
!! drawing from nearby to the east [H L ~> m2 or kg m-1].
real :: FA_u_W0 !< The effective open face area for zonal barotropic transport
!! drawing from nearby to the west [H L ~> m2 or kg m-1].
real :: FA_u_WW !< The effective open face area for zonal barotropic transport
!! drawing from locations far to the west [H L ~> m2 or kg m-1].
real :: uBT_WW !< uBT_WW is the barotropic velocity [L T-1 ~> m s-1], or with INTEGRAL_BT_CONTINUITY
!! the time-integrated barotropic velocity [L ~> m], beyond which the marginal
!! open face area is FA_u_WW. uBT_WW must be non-negative.
real :: uBT_EE !< uBT_EE is a barotropic velocity [L T-1 ~> m s-1], or with INTEGRAL_BT_CONTINUITY
!! the time-integrated barotropic velocity [L ~> m], beyond which the marginal
!! open face area is FA_u_EE. uBT_EE must be non-positive.
real :: uh_crvW !< The curvature of face area with velocity for flow from the west [H T2 L-1 ~> s2 or kg s2 m-3]
!! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY.
real :: uh_crvE !< The curvature of face area with velocity for flow from the east [H T2 L-1 ~> s2 or kg s2 m-3]
!! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY.
real :: uh_WW !< The zonal transport when ubt=ubt_WW [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent
!! time-integrated transport with INTEGRAL_BT_CONTINUITY [H L2 ~> m3 or kg].
real :: uh_EE !< The zonal transport when ubt=ubt_EE [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent
!! time-integrated transport with INTEGRAL_BT_CONTINUITY [H L2 ~> m3 or kg].
end type local_BT_cont_u_type
!> A desciption of the functional dependence of transport at a v-point
type, private :: local_BT_cont_v_type
real :: FA_v_NN !< The effective open face area for meridional barotropic transport
!! drawing from locations far to the north [H L ~> m2 or kg m-1].
real :: FA_v_N0 !< The effective open face area for meridional barotropic transport
!! drawing from nearby to the north [H L ~> m2 or kg m-1].
real :: FA_v_S0 !< The effective open face area for meridional barotropic transport
!! drawing from nearby to the south [H L ~> m2 or kg m-1].
real :: FA_v_SS !< The effective open face area for meridional barotropic transport
!! drawing from locations far to the south [H L ~> m2 or kg m-1].
real :: vBT_SS !< vBT_SS is the barotropic velocity [L T-1 ~> m s-1], or with INTEGRAL_BT_CONTINUITY
!! the time-integrated barotropic velocity [L ~> m], beyond which the marginal
!! open face area is FA_v_SS. vBT_SS must be non-negative.
real :: vBT_NN !< vBT_NN is the barotropic velocity [L T-1 ~> m s-1], or with INTEGRAL_BT_CONTINUITY
!! the time-integrated barotropic velocity [L ~> m], beyond which the marginal
!! open face area is FA_v_NN. vBT_NN must be non-positive.
real :: vh_crvS !< The curvature of face area with velocity for flow from the south [H T2 L-1 ~> s2 or kg s2 m-3]
!! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY.
real :: vh_crvN !< The curvature of face area with velocity for flow from the north [H T2 L-1 ~> s2 or kg s2 m-3]
!! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY.
real :: vh_SS !< The meridional transport when vbt=vbt_SS [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent
!! time-integrated transport with INTEGRAL_BT_CONTINUITY [H L2 ~> m3 or kg].
real :: vh_NN !< The meridional transport when vbt=vbt_NN [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent
!! time-integrated transport with INTEGRAL_BT_CONTINUITY [H L2 ~> m3 or kg].
end type local_BT_cont_v_type
!> A container for passing around active tracer point memory limits
type, private :: memory_size_type
!>@{ Currently active memory limits
integer :: isdw, iedw, jsdw, jedw ! The memory limits of the wide halo arrays.
!>@}
end type memory_size_type
!>@{ CPU time clock IDs
integer :: id_clock_sync=-1, id_clock_calc=-1
integer :: id_clock_calc_pre=-1, id_clock_calc_post=-1
integer :: id_clock_pass_step=-1, id_clock_pass_pre=-1, id_clock_pass_post=-1
!>@}
!>@{ Enumeration values for various schemes
integer, parameter :: HARMONIC = 1
integer, parameter :: ARITHMETIC = 2
integer, parameter :: HYBRID = 3
integer, parameter :: FROM_BT_CONT = 4
integer, parameter :: HYBRID_BT_CONT = 5
character*(20), parameter :: HYBRID_STRING = "HYBRID"
character*(20), parameter :: HARMONIC_STRING = "HARMONIC"
character*(20), parameter :: ARITHMETIC_STRING = "ARITHMETIC"
character*(20), parameter :: BT_CONT_STRING = "FROM_BT_CONT"
!>@}
contains
!> This subroutine time steps the barotropic equations explicitly.
!! For gravity waves, anything between a forwards-backwards scheme
!! and a simulated backwards Euler scheme is used, with bebt between
!! 0.0 and 1.0 determining the scheme. In practice, bebt must be of
!! order 0.2 or greater. A forwards-backwards treatment of the
!! Coriolis terms is always used.
subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, &
eta_PF_in, U_Cor, V_Cor, accel_layer_u, accel_layer_v, &
eta_out, uhbtav, vhbtav, G, GV, US, CS, &
visc_rem_u, visc_rem_v, etaav, ADp, OBC, BT_cont, eta_PF_start, &
taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0)
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: U_in !< The initial (3-D) zonal
!! velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: V_in !< The initial (3-D) meridional
!! velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eta_in !< The initial barotropic free surface height
!! anomaly or column mass anomaly [H ~> m or kg m-2].
real, intent(in) :: dt !< The time increment to integrate over [T ~> s].
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: bc_accel_u !< The zonal baroclinic accelerations,
!! [L T-2 ~> m s-2].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: bc_accel_v !< The meridional baroclinic accelerations,
!! [L T-2 ~> m s-2].
type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: pbce !< The baroclinic pressure anomaly in each layer
!! due to free surface height anomalies
!! [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eta_PF_in !< The 2-D eta field (either SSH anomaly or
!! column mass anomaly) that was used to calculate the input
!! pressure gradient accelerations (or its final value if
!! eta_PF_start is provided [H ~> m or kg m-2].
!! Note: eta_in, pbce, and eta_PF_in must have up-to-date
!! values in the first point of their halos.
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: U_Cor !< The (3-D) zonal velocities used to
!! calculate the Coriolis terms in bc_accel_u [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: V_Cor !< The (3-D) meridional velocities used to
!! calculate the Coriolis terms in bc_accel_u [L T-1 ~> m s-1].
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: accel_layer_u !< The zonal acceleration of each layer due
!! to the barotropic calculation [L T-2 ~> m s-2].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: accel_layer_v !< The meridional acceleration of each layer
!! due to the barotropic calculation [L T-2 ~> m s-2].
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_out !< The final barotropic free surface
!! height anomaly or column mass anomaly [H ~> m or kg m-2].
real, dimension(SZIB_(G),SZJ_(G)), intent(out) :: uhbtav !< the barotropic zonal volume or mass
!! fluxes averaged through the barotropic steps
!! [H L2 T-1 ~> m3 or kg s-1].
real, dimension(SZI_(G),SZJB_(G)), intent(out) :: vhbtav !< the barotropic meridional volume or mass
!! fluxes averaged through the barotropic steps
!! [H L2 T-1 ~> m3 or kg s-1].
type(barotropic_CS), pointer :: CS !< The control structure returned by a
!! previous call to barotropic_init.
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: visc_rem_u !< Both the fraction of the momentum
!! originally in a layer that remains after a time-step of
!! viscosity, and the fraction of a time-step's worth of a
!! barotropic acceleration that a layer experiences after
!! viscosity is applied, in the zonal direction. Nondimensional
!! between 0 (at the bottom) and 1 (far above the bottom).
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: visc_rem_v !< Ditto for meridional direction [nondim].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: etaav !< The free surface height or column mass
!! averaged over the barotropic integration [H ~> m or kg m-2].
type(accel_diag_ptrs), optional, pointer :: ADp !< Acceleration diagnostic pointers
type(ocean_OBC_type), optional, pointer :: OBC !< The open boundary condition structure.
type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe
!! the effective open face areas as a function of barotropic
!! flow.
real, dimension(:,:), optional, pointer :: eta_PF_start !< The eta field consistent with the pressure
!! gradient at the start of the barotropic stepping
!! [H ~> m or kg m-2].
real, dimension(:,:), optional, pointer :: taux_bot !< The zonal bottom frictional stress from
!! ocean to the seafloor [R L Z T-2 ~> Pa].
real, dimension(:,:), optional, pointer :: tauy_bot !< The meridional bottom frictional stress
!! from ocean to the seafloor [R L Z T-2 ~> Pa].
real, dimension(:,:,:), optional, pointer :: uh0 !< The zonal layer transports at reference
!! velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(:,:,:), optional, pointer :: u_uh0 !< The velocities used to calculate
!! uh0 [L T-1 ~> m s-1]
real, dimension(:,:,:), optional, pointer :: vh0 !< The zonal layer transports at reference
!! velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(:,:,:), optional, pointer :: v_vh0 !< The velocities used to calculate
!! vh0 [L T-1 ~> m s-1]
! Local variables
real :: ubt_Cor(SZIB_(G),SZJ_(G)) ! The barotropic velocities that had been
real :: vbt_Cor(SZI_(G),SZJB_(G)) ! used to calculate the input Coriolis
! terms [L T-1 ~> m s-1].
real :: wt_u(SZIB_(G),SZJ_(G),SZK_(GV)) ! wt_u and wt_v are the
real :: wt_v(SZI_(G),SZJB_(G),SZK_(GV)) ! normalized weights to
! be used in calculating barotropic velocities, possibly with
! sums less than one due to viscous losses. Nondimensional.
real, dimension(SZIB_(G),SZJ_(G)) :: &
av_rem_u, & ! The weighted average of visc_rem_u, nondimensional.
tmp_u, & ! A temporary array at u points.
ubt_st, & ! The zonal barotropic velocity at the start of timestep [L T-1 ~> m s-1].
ubt_dt ! The zonal barotropic velocity tendency [L T-2 ~> m s-2].
real, dimension(SZI_(G),SZJB_(G)) :: &
av_rem_v, & ! The weighted average of visc_rem_v, nondimensional.
tmp_v, & ! A temporary array at v points.
vbt_st, & ! The meridional barotropic velocity at the start of timestep [L T-1 ~> m s-1].
vbt_dt ! The meridional barotropic velocity tendency [L T-2 ~> m s-2].
real, dimension(SZI_(G),SZJ_(G)) :: &
tmp_h, & ! A temporary array at h points.
e_anom ! The anomaly in the sea surface height or column mass
! averaged between the beginning and end of the time step,
! relative to eta_PF, with SAL effects included [H ~> m or kg m-2].
! These are always allocated with symmetric memory and wide halos.
real :: q(SZIBW_(CS),SZJBW_(CS)) ! A pseudo potential vorticity [T-1 Z-1 ~> s-1 m-1]
! or [T-1 H-1 ~> s-1 m-1 or m2 s-1 kg-1]
real, dimension(SZIBW_(CS),SZJW_(CS)) :: &
ubt, & ! The zonal barotropic velocity [L T-1 ~> m s-1].
bt_rem_u, & ! The fraction of the barotropic zonal velocity that remains
! after a time step, the remainder being lost to bottom drag.
! bt_rem_u is a nondimensional number between 0 and 1.
BT_force_u, & ! The vertical average of all of the u-accelerations that are
! not explicitly included in the barotropic equation [L T-2 ~> m s-2].
u_accel_bt, & ! The difference between the zonal acceleration from the
! barotropic calculation and BT_force_u [L T-2 ~> m s-2].
uhbt, & ! The zonal barotropic thickness fluxes [H L2 T-1 ~> m3 s-1 or kg s-1].
uhbt0, & ! The difference between the sum of the layer zonal thickness
! fluxes and the barotropic thickness flux using the same
! velocity [H L2 T-1 ~> m3 s-1 or kg s-1].
ubt_old, & ! The starting value of ubt in a barotropic step [L T-1 ~> m s-1].
ubt_first, & ! The starting value of ubt in a series of barotropic steps [L T-1 ~> m s-1].
ubt_sum, & ! The sum of ubt over the time steps [L T-1 ~> m s-1].
ubt_int, & ! The running time integral of ubt over the time steps [L ~> m].
uhbt_sum, & ! The sum of uhbt over the time steps [H L2 T-1 ~> m3 s-1 or kg s-1].
uhbt_int, & ! The running time integral of uhbt over the time steps [H L2 ~> m3].
ubt_wtd, & ! A weighted sum used to find the filtered final ubt [L T-1 ~> m s-1].
ubt_trans, & ! The latest value of ubt used for a transport [L T-1 ~> m s-1].
azon, bzon, & ! _zon & _mer are the values of the Coriolis force which
czon, dzon, & ! are applied to the neighboring values of vbtav & ubtav,
amer, bmer, & ! respectively to get the barotropic inertial rotation
cmer, dmer, & ! [T-1 ~> s-1].
Cor_u, & ! The zonal Coriolis acceleration [L T-2 ~> m s-2].
Cor_ref_u, & ! The zonal barotropic Coriolis acceleration due
! to the reference velocities [L T-2 ~> m s-2].
PFu, & ! The zonal pressure force acceleration [L T-2 ~> m s-2].
Rayleigh_u, & ! A Rayleigh drag timescale operating at u-points [T-1 ~> s-1].
PFu_bt_sum, & ! The summed zonal barotropic pressure gradient force [L T-2 ~> m s-2].
Coru_bt_sum, & ! The summed zonal barotropic Coriolis acceleration [L T-2 ~> m s-2].
DCor_u, & ! An averaged depth or total thickness at u points [Z ~> m] or [H ~> m or kg m-2].
Datu ! Basin depth at u-velocity grid points times the y-grid
! spacing [H L ~> m2 or kg m-1].
real, dimension(SZIW_(CS),SZJBW_(CS)) :: &
vbt, & ! The meridional barotropic velocity [L T-1 ~> m s-1].
bt_rem_v, & ! The fraction of the barotropic meridional velocity that
! remains after a time step, the rest being lost to bottom
! drag. bt_rem_v is a nondimensional number between 0 and 1.
BT_force_v, & ! The vertical average of all of the v-accelerations that are
! not explicitly included in the barotropic equation [L T-2 ~> m s-2].
v_accel_bt, & ! The difference between the meridional acceleration from the
! barotropic calculation and BT_force_v [L T-2 ~> m s-2].
vhbt, & ! The meridional barotropic thickness fluxes [H L2 T-1 ~> m3 s-1 or kg s-1].
vhbt0, & ! The difference between the sum of the layer meridional
! thickness fluxes and the barotropic thickness flux using
! the same velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
vbt_old, & ! The starting value of vbt in a barotropic step [L T-1 ~> m s-1].
vbt_first, & ! The starting value of ubt in a series of barotropic steps [L T-1 ~> m s-1].
vbt_sum, & ! The sum of vbt over the time steps [L T-1 ~> m s-1].
vbt_int, & ! The running time integral of vbt over the time steps [L ~> m].
vhbt_sum, & ! The sum of vhbt over the time steps [H L2 T-1 ~> m3 s-1 or kg s-1].
vhbt_int, & ! The running time integral of vhbt over the time steps [H L2 ~> m3].
vbt_wtd, & ! A weighted sum used to find the filtered final vbt [L T-1 ~> m s-1].
vbt_trans, & ! The latest value of vbt used for a transport [L T-1 ~> m s-1].
Cor_v, & ! The meridional Coriolis acceleration [L T-2 ~> m s-2].
Cor_ref_v, & ! The meridional barotropic Coriolis acceleration due
! to the reference velocities [L T-2 ~> m s-2].
PFv, & ! The meridional pressure force acceleration [L T-2 ~> m s-2].
Rayleigh_v, & ! A Rayleigh drag timescale operating at v-points [T-1 ~> s-1].
PFv_bt_sum, & ! The summed meridional barotropic pressure gradient force,
! [L T-2 ~> m s-2].
Corv_bt_sum, & ! The summed meridional barotropic Coriolis acceleration,
! [L T-2 ~> m s-2].
DCor_v, & ! An averaged depth or total thickness at v points [Z ~> m] or [H ~> m or kg m-2].
Datv ! Basin depth at v-velocity grid points times the x-grid
! spacing [H L ~> m2 or kg m-1].
real, target, dimension(SZIW_(CS),SZJW_(CS)) :: &
eta, & ! The barotropic free surface height anomaly or column mass
! anomaly [H ~> m or kg m-2]
eta_pred ! A predictor value of eta [H ~> m or kg m-2] like eta.
real, dimension(:,:), pointer :: &
eta_PF_BT ! A pointer to the eta array (either eta or eta_pred) that
! determines the barotropic pressure force [H ~> m or kg m-2]
real, dimension(SZIW_(CS),SZJW_(CS)) :: &
eta_sum, & ! eta summed across the timesteps [H ~> m or kg m-2].
eta_wtd, & ! A weighted estimate used to calculate eta_out [H ~> m or kg m-2].
eta_IC, & ! A local copy of the initial 2-D eta field (eta_in) [H ~> m or kg m-2]
eta_PF, & ! A local copy of the 2-D eta field (either SSH anomaly or
! column mass anomaly) that was used to calculate the input
! pressure gradient accelerations [H ~> m or kg m-2].
eta_PF_1, & ! The initial value of eta_PF, when interp_eta_PF is
! true [H ~> m or kg m-2].
d_eta_PF, & ! The change in eta_PF over the barotropic time stepping when
! interp_eta_PF is true [H ~> m or kg m-2].
gtot_E, & ! gtot_X is the effective total reduced gravity used to relate
gtot_W, & ! free surface height deviations to pressure forces (including
gtot_N, & ! GFS and baroclinic contributions) in the barotropic momentum
gtot_S, & ! equations half a grid-point in the X-direction (X is N, S, E, or W)
! from the thickness point [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
! (See Hallberg, J Comp Phys 1997 for a discussion.)
eta_src, & ! The source of eta per barotropic timestep [H ~> m or kg m-2].
dyn_coef_eta, & ! The coefficient relating the changes in eta to the
! dynamic surface pressure under rigid ice
! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
p_surf_dyn ! A dynamic surface pressure under rigid ice [L2 T-2 ~> m2 s-2].
type(local_BT_cont_u_type), dimension(SZIBW_(CS),SZJW_(CS)) :: &
BTCL_u ! A repackaged version of the u-point information in BT_cont.
type(local_BT_cont_v_type), dimension(SZIW_(CS),SZJBW_(CS)) :: &
BTCL_v ! A repackaged version of the v-point information in BT_cont.
! End of wide-sized variables.
real, dimension(SZIBW_(CS),SZJW_(CS)) :: &
ubt_prev, ubt_sum_prev, ubt_wtd_prev, & ! Previous velocities stored for OBCs [L T-1 ~> m s-1]
uhbt_prev, uhbt_sum_prev, & ! Previous transports stored for OBCs [L2 H T-1 ~> m3 s-1]
ubt_int_prev, & ! Previous value of time-integrated velocity stored for OBCs [L ~> m]
uhbt_int_prev ! Previous value of time-integrated transport stored for OBCs [L2 H ~> m3]
real, dimension(SZIW_(CS),SZJBW_(CS)) :: &
vbt_prev, vbt_sum_prev, vbt_wtd_prev, & ! Previous velocities stored for OBCs [L T-1 ~> m s-1]
vhbt_prev, vhbt_sum_prev, & ! Previous transports stored for OBCs [L2 H T-1 ~> m3 s-1]
vbt_int_prev, & ! Previous value of time-integrated velocity stored for OBCs [L ~> m]
vhbt_int_prev ! Previous value of time-integrated transport stored for OBCs [L2 H ~> m3]
real :: mass_to_Z ! The depth unit conversion divided by the mean density (Rho0) [Z m-1 R-1 ~> m3 kg-1].
real :: mass_accel_to_Z ! The inverse of the mean density (Rho0) [R-1 ~> m3 kg-1].
real :: visc_rem ! A work variable that may equal visc_rem_[uv]. Nondim.
real :: vel_prev ! The previous velocity [L T-1 ~> m s-1].
real :: dtbt ! The barotropic time step [T ~> s].
real :: bebt ! A copy of CS%bebt [nondim].
real :: be_proj ! The fractional amount by which velocities are projected
! when project_velocity is true. For now be_proj is set
! to equal bebt, as they have similar roles and meanings.
real :: Idt ! The inverse of dt [T-1 ~> s-1].
real :: det_de ! The partial derivative due to self-attraction and loading
! of the reference geopotential with the sea surface height.
! This is typically ~0.09 or less.
real :: dgeo_de ! The constant of proportionality between geopotential and
! sea surface height. It is a nondimensional number of
! order 1. For stability, this may be made larger
! than the physical problem would suggest.
real :: Instep ! The inverse of the number of barotropic time steps to take.
real :: wt_end ! The weighting of the final value of eta_PF [nondim]
integer :: nstep ! The number of barotropic time steps to take.
type(time_type) :: &
time_bt_start, & ! The starting time of the barotropic steps.
time_step_end, & ! The end time of a barotropic step.
time_end_in ! The end time for diagnostics when this routine started.
real :: time_int_in ! The diagnostics' time interval when this routine started.
real :: Htot_avg ! The average total thickness of the tracer columns adjacent to a
! velocity point [H ~> m or kg m-2]
logical :: do_hifreq_output ! If true, output occurs every barotropic step.
logical :: use_BT_cont, do_ave, find_etaav, find_PF, find_Cor
logical :: integral_BT_cont ! If true, update the barotropic continuity equation directly
! from the initial condition using the time-integrated barotropic velocity.
logical :: ice_is_rigid, nonblock_setup, interp_eta_PF
logical :: project_velocity, add_uh0
real :: dyn_coef_max ! The maximum stable value of dyn_coef_eta
! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
real :: ice_strength = 0.0 ! The effective strength of the ice [L2 Z-1 T-2 ~> m s-2].
real :: Idt_max2 ! The squared inverse of the local maximum stable
! barotropic time step [T-2 ~> s-2].
real :: H_min_dyn ! The minimum depth to use in limiting the size of the
! dynamic surface pressure for stability [H ~> m or kg m-2].
real :: H_eff_dx2 ! The effective total thickness divided by the grid spacing
! squared [H L-2 ~> m-1 or kg m-4].
real :: u_max_cor, v_max_cor ! The maximum corrective velocities [L T-1 ~> m s-1].
real :: uint_cor, vint_cor ! The maximum time-integrated corrective velocities [L ~> m].
real :: Htot ! The total thickness [H ~> m or kg m-2].
real :: eta_cor_max ! The maximum fluid that can be added as a correction to eta [H ~> m or kg m-2].
real :: accel_underflow ! An acceleration that is so small it should be zeroed out [L T-2 ~> m s-2].
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m or kg m-2].
real :: Idtbt ! The inverse of the barotropic time step [T-1 ~> s-1]
real, allocatable, dimension(:) :: wt_vel, wt_eta, wt_accel, wt_trans, wt_accel2
real :: sum_wt_vel, sum_wt_eta, sum_wt_accel, sum_wt_trans
real :: I_sum_wt_vel, I_sum_wt_eta, I_sum_wt_accel, I_sum_wt_trans
real :: dt_filt ! The half-width of the barotropic filter [T ~> s].
real :: trans_wt1, trans_wt2 ! The weights used to compute ubt_trans and vbt_trans
integer :: nfilter
logical :: apply_OBCs, apply_OBC_flather, apply_OBC_open
type(memory_size_type) :: MS
character(len=200) :: mesg
integer :: isv, iev, jsv, jev ! The valid array size at the end of a step.
integer :: stencil ! The stencil size of the algorithm, often 1 or 2.
integer :: isvf, ievf, jsvf, jevf, num_cycles
integer :: err_count ! A counter to limit the volume of error messages written to stdout.
integer :: i, j, k, n
integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
integer :: ioff, joff
integer :: l_seg
if (.not.associated(CS)) call MOM_error(FATAL, &
"btstep: Module MOM_barotropic must be initialized before it is used.")
if (.not.CS%split) return
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%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
MS%isdw = CS%isdw ; MS%iedw = CS%iedw ; MS%jsdw = CS%jsdw ; MS%jedw = CS%jedw
h_neglect = GV%H_subroundoff
err_count = 0
Idt = 1.0 / dt
accel_underflow = CS%vel_underflow * Idt
use_BT_cont = .false.
if (present(BT_cont)) use_BT_cont = (associated(BT_cont))
integral_BT_cont = use_BT_cont .and. CS%integral_BT_cont
interp_eta_PF = .false.
if (present(eta_PF_start)) interp_eta_PF = (associated(eta_PF_start))
project_velocity = CS%BT_project_velocity
! Figure out the fullest arrays that could be updated.
stencil = 1
if ((.not.use_BT_cont) .and. CS%Nonlinear_continuity .and. &
(CS%Nonlin_cont_update_period > 0)) stencil = 2
do_ave = query_averaging_enabled(CS%diag)
find_etaav = present(etaav)
find_PF = (do_ave .and. ((CS%id_PFu_bt > 0) .or. (CS%id_PFv_bt > 0)))
find_Cor = (do_ave .and. ((CS%id_Coru_bt > 0) .or. (CS%id_Corv_bt > 0)))
add_uh0 = .false.
if (present(uh0)) add_uh0 = associated(uh0)
if (add_uh0 .and. .not.(present(vh0) .and. present(u_uh0) .and. &
present(v_vh0))) call MOM_error(FATAL, &
"btstep: vh0, u_uh0, and v_vh0 must be present if uh0 is used.")
if (add_uh0 .and. .not.(associated(vh0) .and. associated(u_uh0) .and. &
associated(v_vh0))) call MOM_error(FATAL, &
"btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.")
! This can be changed to try to optimize the performance.
nonblock_setup = G%nonblocking_updates
if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)
apply_OBCs = .false. ; CS%BT_OBC%apply_u_OBCs = .false. ; CS%BT_OBC%apply_v_OBCs = .false.
apply_OBC_open = .false.
apply_OBC_flather = .false.
if (present(OBC)) then ; if (associated(OBC)) then
CS%BT_OBC%apply_u_OBCs = OBC%open_u_BCs_exist_globally .or. OBC%specified_u_BCs_exist_globally
CS%BT_OBC%apply_v_OBCs = OBC%open_v_BCs_exist_globally .or. OBC%specified_v_BCs_exist_globally
apply_OBC_flather = open_boundary_query(OBC, apply_Flather_OBC=.true.)
apply_OBC_open = open_boundary_query(OBC, apply_open_OBC=.true.)
apply_OBCs = open_boundary_query(OBC, apply_specified_OBC=.true.) .or. &
apply_OBC_flather .or. apply_OBC_open
if (apply_OBC_flather .and. .not.GV%Boussinesq) call MOM_error(FATAL, &
"btstep: Flather open boundary conditions have not yet been "// &
"implemented for a non-Boussinesq model.")
endif ; endif
num_cycles = 1
if (CS%use_wide_halos) &
num_cycles = min((is-CS%isdw) / stencil, (js-CS%jsdw) / stencil)
isvf = is - (num_cycles-1)*stencil ; ievf = ie + (num_cycles-1)*stencil
jsvf = js - (num_cycles-1)*stencil ; jevf = je + (num_cycles-1)*stencil
nstep = CEILING(dt/CS%dtbt - 0.0001)
if (is_root_PE() .and. (nstep /= CS%nstep_last)) then
write(mesg,'("btstep is using a dynamic barotropic timestep of ", ES12.6, &
& " seconds, max ", ES12.6, ".")') (US%T_to_s*dt/nstep), US%T_to_s*CS%dtbt_max
call MOM_mesg(mesg, 3)
endif
CS%nstep_last = nstep
! Set the actual barotropic time step.
Instep = 1.0 / real(nstep)
dtbt = dt * Instep
Idtbt = 1.0 / dtbt
bebt = CS%bebt
be_proj = CS%bebt
mass_accel_to_Z = 1.0 / GV%Rho0
mass_to_Z = US%m_to_Z / GV%Rho0
!--- setup the weight when computing vbt_trans and ubt_trans
if (project_velocity) then
trans_wt1 = (1.0 + be_proj); trans_wt2 = -be_proj
else
trans_wt1 = bebt ; trans_wt2 = (1.0-bebt)
endif
do_hifreq_output = .false.
if ((CS%id_ubt_hifreq > 0) .or. (CS%id_vbt_hifreq > 0) .or. &
(CS%id_eta_hifreq > 0) .or. (CS%id_eta_pred_hifreq > 0) .or. &
(CS%id_uhbt_hifreq > 0) .or. (CS%id_vhbt_hifreq > 0)) then
do_hifreq_output = query_averaging_enabled(CS%diag, time_int_in, time_end_in)
if (do_hifreq_output) &
time_bt_start = time_end_in - real_to_time(US%T_to_s*dt)
endif
!--- begin setup for group halo update
if (id_clock_pass_pre > 0) call cpu_clock_begin(id_clock_pass_pre)
if (.not. CS%linearized_BT_PV) then
call create_group_pass(CS%pass_q_DCor, q, CS%BT_Domain, To_All, position=CORNER)
call create_group_pass(CS%pass_q_DCor, DCor_u, DCor_v, CS%BT_Domain, &
To_All+Scalar_Pair)
endif
if ((Isq > is-1) .or. (Jsq > js-1)) &
call create_group_pass(CS%pass_tmp_uv, tmp_u, tmp_v, G%Domain)
call create_group_pass(CS%pass_gtot, gtot_E, gtot_N, CS%BT_Domain, &
To_All+Scalar_Pair, AGRID)
call create_group_pass(CS%pass_gtot, gtot_W, gtot_S, CS%BT_Domain, &
To_All+Scalar_Pair, AGRID)
if (CS%dynamic_psurf) &
call create_group_pass(CS%pass_eta_bt_rem, dyn_coef_eta, CS%BT_Domain)
if (interp_eta_PF) then
call create_group_pass(CS%pass_eta_bt_rem, eta_PF_1, CS%BT_Domain)
call create_group_pass(CS%pass_eta_bt_rem, d_eta_PF, CS%BT_Domain)
else
call create_group_pass(CS%pass_eta_bt_rem, eta_PF, CS%BT_Domain)
endif
if (integral_BT_cont) &
call create_group_pass(CS%pass_eta_bt_rem, eta_IC, CS%BT_Domain)
call create_group_pass(CS%pass_eta_bt_rem, eta_src, CS%BT_Domain)
! The following halo updates are not needed without wide halos. RWH
! We do need them after all.
! if (ievf > ie) then
call create_group_pass(CS%pass_eta_bt_rem, bt_rem_u, bt_rem_v, &
CS%BT_Domain, To_All+Scalar_Pair)
if (CS%linear_wave_drag) &
call create_group_pass(CS%pass_eta_bt_rem, Rayleigh_u, Rayleigh_v, &
CS%BT_Domain, To_All+Scalar_Pair)
! endif
! The following halo update is not needed without wide halos. RWH
if (((G%isd > CS%isdw) .or. (G%jsd > CS%jsdw)) .or. (Isq <= is-1) .or. (Jsq <= js-1)) &
call create_group_pass(CS%pass_force_hbt0_Cor_ref, BT_force_u, BT_force_v, CS%BT_Domain)
if (add_uh0) call create_group_pass(CS%pass_force_hbt0_Cor_ref, uhbt0, vhbt0, CS%BT_Domain)
call create_group_pass(CS%pass_force_hbt0_Cor_ref, Cor_ref_u, Cor_ref_v, CS%BT_Domain)
if (.not. use_BT_cont) then
call create_group_pass(CS%pass_Dat_uv, Datu, Datv, CS%BT_Domain, To_All+Scalar_Pair)
endif
call create_group_pass(CS%pass_eta_ubt, eta, CS%BT_Domain)
call create_group_pass(CS%pass_eta_ubt, ubt, vbt, CS%BT_Domain)
if (integral_BT_cont) then
call create_group_pass(CS%pass_eta_ubt, ubt_int, vbt_int, CS%BT_Domain)
! This is only needed with integral_BT_cont, OBCs and multiple barotropic steps between halo updates.
if (apply_OBC_open) &
call create_group_pass(CS%pass_eta_ubt, uhbt_int, vhbt_int, CS%BT_Domain)
endif
call create_group_pass(CS%pass_ubt_Cor, ubt_Cor, vbt_Cor, G%Domain)
! These passes occur at the end of the routine, as data is being readied to
! share with the main part of the MOM6 code.
if (find_etaav) then
call create_group_pass(CS%pass_etaav, etaav, G%Domain)
endif
call create_group_pass(CS%pass_e_anom, e_anom, G%Domain)
call create_group_pass(CS%pass_ubta_uhbta, CS%ubtav, CS%vbtav, G%Domain)
call create_group_pass(CS%pass_ubta_uhbta, uhbtav, vhbtav, G%Domain)
if (id_clock_pass_pre > 0) call cpu_clock_end(id_clock_pass_pre)
!--- end setup for group halo update
! Calculate the constant coefficients for the Coriolis force terms in the
! barotropic momentum equations. This has to be done quite early to start
! the halo update that needs to be completed before the next calculations.
if (CS%linearized_BT_PV) then
!$OMP parallel do default(shared)
do J=jsvf-2,jevf+1 ; do I=isvf-2,ievf+1
q(I,J) = CS%q_D(I,j)
enddo ; enddo
!$OMP parallel do default(shared)
do j=jsvf-1,jevf+1 ; do I=isvf-2,ievf+1
DCor_u(I,j) = CS%D_u_Cor(I,j)
enddo ; enddo
!$OMP parallel do default(shared)
do J=jsvf-2,jevf+1 ; do i=isvf-1,ievf+1
DCor_v(i,J) = CS%D_v_Cor(i,J)
enddo ; enddo
else
q(:,:) = 0.0 ; DCor_u(:,:) = 0.0 ; DCor_v(:,:) = 0.0
if (GV%Boussinesq) then
!$OMP parallel do default(shared)
do j=js,je ; do I=is-1,ie
DCor_u(I,j) = 0.5 * (max(GV%Z_to_H*G%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + &
max(GV%Z_to_H*G%bathyT(i,j) + eta_in(i,j), 0.0) )
enddo ; enddo
!$OMP parallel do default(shared)
do J=js-1,je ; do i=is,ie
DCor_v(i,J) = 0.5 * (max(GV%Z_to_H*G%bathyT(i,j+1) + eta_in(i+1,j), 0.0) + &
max(GV%Z_to_H*G%bathyT(i,j) + eta_in(i,j), 0.0) )
enddo ; enddo
!$OMP parallel do default(shared)
do J=js-1,je ; do I=is-1,ie
q(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * &
((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / &
(max((G%areaT(i,j) * max(GV%Z_to_H*G%bathyT(i,j) + eta_in(i,j), 0.0) + &
G%areaT(i+1,j+1) * max(GV%Z_to_H*G%bathyT(i+1,j+1) + eta_in(i+1,j+1), 0.0)) + &
(G%areaT(i+1,j) * max(GV%Z_to_H*G%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + &
G%areaT(i,j+1) * max(GV%Z_to_H*G%bathyT(i,j+1) + eta_in(i,j+1), 0.0)), h_neglect) )
enddo ; enddo
else
!$OMP parallel do default(shared)
do j=js,je ; do I=is-1,ie
DCor_u(I,j) = 0.5 * (eta_in(i+1,j) + eta_in(i,j))
enddo ; enddo
!$OMP parallel do default(shared)
do J=js-1,je ; do i=is,ie
DCor_v(i,J) = 0.5 * (eta_in(i,j+1) + eta_in(i,j))
enddo ; enddo
!$OMP parallel do default(shared)
do J=js-1,je ; do I=is-1,ie
q(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * &
((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) / &
(max((G%areaT(i,j) * eta_in(i,j) + G%areaT(i+1,j+1) * eta_in(i+1,j+1)) + &
(G%areaT(i+1,j) * eta_in(i+1,j) + G%areaT(i,j+1) * eta_in(i,j+1)), h_neglect) )
enddo ; enddo
endif
! With very wide halos, q and D need to be calculated on the available data
! domain and then updated onto the full computational domain.
! These calculations can be done almost immediately, but the halo updates
! must be done before the [abcd]mer and [abcd]zon are calculated.
if (id_clock_calc_pre > 0) call cpu_clock_end(id_clock_calc_pre)
if (nonblock_setup) then
call start_group_pass(CS%pass_q_DCor, CS%BT_Domain, clock=id_clock_pass_pre)
else
call do_group_pass(CS%pass_q_DCor, CS%BT_Domain, clock=id_clock_pass_pre)
endif
if (id_clock_calc_pre > 0) call cpu_clock_begin(id_clock_calc_pre)
endif
! Zero out various wide-halo arrays.
!$OMP parallel do default(shared)
do j=CS%jsdw,CS%jedw ; do i=CS%isdw,CS%iedw
gtot_E(i,j) = 0.0 ; gtot_W(i,j) = 0.0
gtot_N(i,j) = 0.0 ; gtot_S(i,j) = 0.0
eta(i,j) = 0.0
eta_PF(i,j) = 0.0
if (interp_eta_PF) then
eta_PF_1(i,j) = 0.0 ; d_eta_PF(i,j) = 0.0
endif
if (integral_BT_cont) then
eta_IC(i,j) = 0.0
endif
p_surf_dyn(i,j) = 0.0
if (CS%dynamic_psurf) dyn_coef_eta(i,j) = 0.0
enddo ; enddo
! The halo regions of various arrays need to be initialized to
! non-NaNs in case the neighboring domains are not part of the ocean.
! Otherwise a halo update later on fills in the correct values.
!$OMP parallel do default(shared)
do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw
Cor_ref_u(I,j) = 0.0 ; BT_force_u(I,j) = 0.0 ; ubt(I,j) = 0.0
Datu(I,j) = 0.0 ; bt_rem_u(I,j) = 0.0 ; uhbt0(I,j) = 0.0
enddo ; enddo
!$OMP parallel do default(shared)
do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw
Cor_ref_v(i,J) = 0.0 ; BT_force_v(i,J) = 0.0 ; vbt(i,J) = 0.0
Datv(i,J) = 0.0 ; bt_rem_v(i,J) = 0.0 ; vhbt0(i,J) = 0.0
enddo ; enddo
if (CS%linear_wave_drag) then
!$OMP parallel do default(shared)
do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw
Rayleigh_u(I,j) = 0.0
enddo ; enddo
!$OMP parallel do default(shared)
do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw
Rayleigh_v(i,J) = 0.0
enddo ; enddo
endif
! Copy input arrays into their wide-halo counterparts.
if (interp_eta_PF) then
!$OMP parallel do default(shared)
do j=G%jsd,G%jed ; do i=G%isd,G%ied ! Was "do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1" but doing so breaks OBC. Not sure why?
eta(i,j) = eta_in(i,j)
eta_PF_1(i,j) = eta_PF_start(i,j)
d_eta_PF(i,j) = eta_PF_in(i,j) - eta_PF_start(i,j)
enddo ; enddo
else
!$OMP parallel do default(shared)
do j=G%jsd,G%jed ; do i=G%isd,G%ied !: Was "do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1" but doing so breaks OBC. Not sure why?
eta(i,j) = eta_in(i,j)
eta_PF(i,j) = eta_PF_in(i,j)
enddo ; enddo
endif
if (integral_BT_cont) then
!$OMP parallel do default(shared)
do j=G%jsd,G%jed ; do i=G%isd,G%ied
eta_IC(i,j) = eta_in(i,j)
enddo ; enddo
endif
!$OMP parallel do default(shared) private(visc_rem)
do k=1,nz ; do j=js,je ; do I=is-1,ie
! rem needs greater than visc_rem_u and 1-Instep/visc_rem_u.
! The 0.5 below is just for safety.
if (visc_rem_u(I,j,k) <= 0.0) then ; visc_rem = 0.0