forked from mom-ocean/MOM6
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMOM.F90
4075 lines (3617 loc) · 204 KB
/
MOM.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
!> The central module of the MOM6 ocean model
module MOM
! This file is part of MOM6. See LICENSE.md for the license.
! Infrastructure modules
use MOM_array_transform, only : rotate_array, rotate_vector
use MOM_debugging, only : MOM_debugging_init, hchksum, uvchksum
use MOM_debugging, only : check_redundant
use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum
use MOM_checksum_packages, only : MOM_accel_chksum, MOM_surface_chksum
use MOM_coms, only : num_PEs
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT
use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
use MOM_diag_mediator, only : diag_mediator_init, enable_averaging, enable_averages
use MOM_diag_mediator, only : diag_mediator_infrastructure_init
use MOM_diag_mediator, only : diag_set_state_ptrs, diag_update_remap_grids
use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
use MOM_diag_mediator, only : register_diag_field, register_cell_measure
use MOM_diag_mediator, only : set_axes_info, diag_ctrl, diag_masks_set
use MOM_diag_mediator, only : set_masks_for_axes
use MOM_diag_mediator, only : diag_grid_storage, diag_grid_storage_init
use MOM_diag_mediator, only : diag_save_grids, diag_restore_grids
use MOM_diag_mediator, only : diag_copy_storage_to_diag, diag_copy_diag_to_storage
use MOM_domains, only : MOM_domains_init
use MOM_domains, only : sum_across_PEs, pass_var, pass_vector
use MOM_domains, only : clone_MOM_domain, deallocate_MOM_domain
use MOM_domains, only : To_North, To_East, To_South, To_West
use MOM_domains, only : To_All, Omit_corners, CGRID_NE, SCALAR_PAIR
use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
use MOM_domains, only : start_group_pass, complete_group_pass, Omit_Corners
use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : read_param, get_param, log_version, param_file_type
use MOM_forcing_type, only : forcing, mech_forcing
use MOM_forcing_type, only : MOM_forcing_chksum, MOM_mech_forcing_chksum
use MOM_get_input, only : Get_MOM_Input, directories
use MOM_io, only : MOM_io_init, vardesc, var_desc
use MOM_io, only : slasher, file_exists, MOM_read_data
use MOM_obsolete_params, only : find_obsolete_params
use MOM_restart, only : register_restart_field, register_restart_pair
use MOM_restart, only : query_initialized, save_restart, restart_registry_lock
use MOM_restart, only : restart_init, is_new_run, determine_is_new_run, MOM_restart_CS
use MOM_spatial_means, only : global_mass_integral
use MOM_time_manager, only : time_type, real_to_time, time_type_to_real, operator(+)
use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/)
use MOM_time_manager, only : operator(>=), operator(==), increment_date
use MOM_time_manager, only : set_time
use MOM_unit_tests, only : unit_tests
! MOM core modules
use MOM_ALE, only : ALE_init, ALE_end, ALE_main, ALE_CS, adjustGridForIntegrity
use MOM_ALE, only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile
use MOM_ALE, only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags
use MOM_ALE_sponge, only : rotate_ALE_sponge, update_ALE_sponge_field
use MOM_barotropic, only : Barotropic_CS
use MOM_boundary_update, only : call_OBC_register, OBC_register_end, update_OBC_CS
use MOM_coord_initialization, only : MOM_initialize_coord
use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member
use MOM_diabatic_driver, only : adiabatic, adiabatic_driver_init, diabatic_driver_end
use MOM_diagnostics, only : calculate_diagnostic_fields, MOM_diagnostics_init
use MOM_diagnostics, only : register_transport_diags, post_transport_diagnostics
use MOM_diagnostics, only : register_surface_diags, write_static_fields
use MOM_diagnostics, only : post_surface_dyn_diags, post_surface_thermo_diags
use MOM_diagnostics, only : diagnostics_CS, surface_diag_IDs, transport_diag_IDs
use MOM_diagnostics, only : MOM_diagnostics_end
use MOM_dynamics_unsplit, only : step_MOM_dyn_unsplit, register_restarts_dyn_unsplit
use MOM_dynamics_unsplit, only : initialize_dyn_unsplit, end_dyn_unsplit
use MOM_dynamics_unsplit, only : MOM_dyn_unsplit_CS
use MOM_dynamics_split_RK2, only : step_MOM_dyn_split_RK2, register_restarts_dyn_split_RK2
use MOM_dynamics_split_RK2, only : initialize_dyn_split_RK2, end_dyn_split_RK2
use MOM_dynamics_split_RK2, only : MOM_dyn_split_RK2_CS
use MOM_dynamics_unsplit_RK2, only : step_MOM_dyn_unsplit_RK2, register_restarts_dyn_unsplit_RK2
use MOM_dynamics_unsplit_RK2, only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2
use MOM_dynamics_unsplit_RK2, only : MOM_dyn_unsplit_RK2_CS
use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
use MOM_EOS, only : EOS_init, calculate_density, calculate_TFreeze, EOS_domain
use MOM_fixed_initialization, only : MOM_initialize_fixed
use MOM_forcing_type, only : allocate_forcing_type, allocate_mech_forcing
use MOM_forcing_type, only : deallocate_mech_forcing, deallocate_forcing_type
use MOM_forcing_type, only : rotate_forcing, rotate_mech_forcing
use MOM_grid, only : ocean_grid_type, MOM_grid_init, MOM_grid_end
use MOM_grid, only : set_first_direction, rescale_grid_bathymetry
use MOM_hor_index, only : hor_index_type, hor_index_init
use MOM_hor_index, only : rotate_hor_index
use MOM_interface_heights, only : find_eta
use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init, VarMix_end
use MOM_lateral_mixing_coeffs, only : calc_resoln_function, calc_depth_function, VarMix_CS
use MOM_MEKE, only : MEKE_alloc_register_restart, step_forward_MEKE
use MOM_MEKE, only : MEKE_CS, MEKE_init, MEKE_end
use MOM_MEKE_types, only : MEKE_type
use MOM_mixed_layer_restrat, only : mixedlayer_restrat, mixedlayer_restrat_init, mixedlayer_restrat_CS
use MOM_mixed_layer_restrat, only : mixedlayer_restrat_register_restarts
use MOM_obsolete_diagnostics, only : register_obsolete_diagnostics
use MOM_open_boundary, only : ocean_OBC_type, OBC_registry_type
use MOM_open_boundary, only : register_temp_salt_segments
use MOM_open_boundary, only : open_boundary_register_restarts
use MOM_open_boundary, only : update_segment_tracer_reservoirs
use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init
use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML
use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS
use MOM_set_visc, only : set_visc_init, set_visc_end
use MOM_shared_initialization, only : write_ocean_geometry_file
use MOM_sponge, only : init_sponge_diags, sponge_CS
use MOM_state_initialization, only : MOM_initialize_state
use MOM_sum_output, only : write_energy, accumulate_net_input
use MOM_sum_output, only : MOM_sum_output_init, MOM_sum_output_end
use MOM_sum_output, only : sum_output_CS
use MOM_ALE_sponge, only : init_ALE_sponge_diags, ALE_sponge_CS
use MOM_thickness_diffuse, only : thickness_diffuse, thickness_diffuse_init
use MOM_thickness_diffuse, only : thickness_diffuse_end, thickness_diffuse_CS
use MOM_tracer_advect, only : advect_tracer, tracer_advect_init
use MOM_tracer_advect, only : tracer_advect_end, tracer_advect_CS
use MOM_tracer_hor_diff, only : tracer_hordiff, tracer_hor_diff_init
use MOM_tracer_hor_diff, only : tracer_hor_diff_end, tracer_hor_diff_CS
use MOM_tracer_registry, only : tracer_registry_type, register_tracer, tracer_registry_init
use MOM_tracer_registry, only : register_tracer_diagnostics, post_tracer_diagnostics_at_sync
use MOM_tracer_registry, only : post_tracer_transport_diagnostics
use MOM_tracer_registry, only : preALE_tracer_diagnostics, postALE_tracer_diagnostics
use MOM_tracer_registry, only : lock_tracer_registry, tracer_registry_end
use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_CS
use MOM_tracer_flow_control, only : tracer_flow_control_init, call_tracer_surface_state
use MOM_tracer_flow_control, only : tracer_flow_control_end
use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid
use MOM_transcribe_grid, only : rotate_dyngrid
use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init
use MOM_unit_scaling, only : unit_scaling_end, fix_restart_unit_scaling
use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state
use MOM_variables, only : thermo_var_ptrs, vertvisc_type
use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state
use MOM_variables, only : rotate_surface_state
use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd
use MOM_verticalGrid, only : fix_restart_scaling
use MOM_verticalGrid, only : get_thickness_units, get_flux_units, get_tr_flux_units
use MOM_wave_interface, only : wave_parameters_CS, waves_end
use MOM_wave_interface, only : Update_Stokes_Drift
! ODA modules
use MOM_oda_driver_mod, only : ODA_CS, oda, init_oda, oda_end
use MOM_oda_driver_mod, only : set_prior_tracer, set_analysis_time, apply_oda_tracer_increments
use MOM_oda_incupd, only : oda_incupd_CS, init_oda_incupd_diags
! Offline modules
use MOM_offline_main, only : offline_transport_CS, offline_transport_init, update_offline_fields
use MOM_offline_main, only : insert_offline_main, extract_offline_main, post_offline_convergence_diags
use MOM_offline_main, only : register_diags_offline_transport, offline_advection_ale
use MOM_offline_main, only : offline_redistribute_residual, offline_diabatic_ale
use MOM_offline_main, only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean
use MOM_offline_main, only : offline_advection_layer, offline_transport_end
use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline
use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query, initialize_ice_shelf
use MOM_particles_mod, only : particles, particles_init, particles_run, particles_save_restart, particles_end
implicit none ; private
#include <MOM_memory.h>
! 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.
!> A structure with diagnostic IDs of the state variables
type MOM_diag_IDs
!>@{ 3-d state field diagnostic IDs
integer :: id_u = -1, id_v = -1, id_h = -1
!>@}
!> 2-d state field diagnotic ID
integer :: id_ssh_inst = -1
end type MOM_diag_IDs
!> Control structure for the MOM module, including the variables that describe
!! the state of the ocean.
type, public :: MOM_control_struct ; private
real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: &
h, & !< layer thickness [H ~> m or kg m-2]
T, & !< potential temperature [degC]
S !< salinity [ppt]
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
u, & !< zonal velocity component [L T-1 ~> m s-1]
uh, & !< uh = u * h * dy at u grid points [H L2 T-1 ~> m3 s-1 or kg s-1]
uhtr !< accumulated zonal thickness fluxes to advect tracers [H L2 ~> m3 or kg]
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
v, & !< meridional velocity [L T-1 ~> m s-1]
vh, & !< vh = v * h * dx at v grid points [H L2 T-1 ~> m3 s-1 or kg s-1]
vhtr !< accumulated meridional thickness fluxes to advect tracers [H L2 ~> m3 or kg]
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: ssh_rint
!< A running time integral of the sea surface height [T m ~> s m].
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: ave_ssh_ibc
!< time-averaged (over a forcing time step) sea surface height
!! with a correction for the inverse barometer [m]
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: eta_av_bc
!< free surface height or column mass time averaged over the last
!! baroclinic dynamics time step [H ~> m or kg m-2]
real, dimension(:,:), pointer :: &
Hml => NULL() !< active mixed layer depth [Z ~> m]
real :: time_in_cycle !< The running time of the current time-stepping cycle
!! in calls that step the dynamics, and also the length of
!! the time integral of ssh_rint [T ~> s].
real :: time_in_thermo_cycle !< The running time of the current time-stepping
!! cycle in calls that step the thermodynamics [T ~> s].
type(ocean_grid_type) :: G_in !< Input grid metric
type(ocean_grid_type), pointer :: G => NULL() !< Model grid metric
logical :: rotate_index = .false. !< True if index map is rotated
type(verticalGrid_type), pointer :: &
GV => NULL() !< structure containing vertical grid info
type(unit_scale_type), pointer :: &
US => NULL() !< structure containing various unit conversion factors
type(thermo_var_ptrs) :: tv !< structure containing pointers to available thermodynamic fields
real :: t_dyn_rel_adv !< The time of the dynamics relative to tracer advection and lateral mixing
!! [T ~> s], or equivalently the elapsed time since advectively updating the
!! tracers. t_dyn_rel_adv is invariably positive and may span multiple coupling timesteps.
integer :: n_dyn_steps_in_adv !< The number of dynamics time steps that contributed to uhtr
!! and vhtr since the last time tracer advection occured.
real :: t_dyn_rel_thermo !< The time of the dynamics relative to diabatic processes and remapping
!! [T ~> s]. t_dyn_rel_thermo can be negative or positive depending on whether
!! the diabatic processes are applied before or after the dynamics and may span
!! multiple coupling timesteps.
real :: t_dyn_rel_diag !< The time of the diagnostics relative to diabatic processes and remapping
!! [T ~> s]. t_dyn_rel_diag is always positive, since the diagnostics must lag.
logical :: preadv_h_stored = .false. !< If true, the thicknesses from before the advective cycle
!! have been stored for use in diagnostics.
type(diag_ctrl) :: diag !< structure to regulate diagnostic output timing
type(vertvisc_type) :: visc !< structure containing vertical viscosities,
!! bottom drag viscosities, and related fields
type(MEKE_type), pointer :: MEKE => NULL() !< structure containing fields
!! related to the Mesoscale Eddy Kinetic Energy
logical :: adiabatic !< If true, there are no diapycnal mass fluxes, and no calls
!! to routines to calculate or apply diapycnal fluxes.
logical :: diabatic_first !< If true, apply diabatic and thermodynamic processes before time
!! stepping the dynamics.
logical :: use_ALE_algorithm !< If true, use the ALE algorithm rather than layered
!! isopycnal/stacked shallow water mode. This logical is set by calling the
!! function useRegridding() from the MOM_regridding module.
logical :: alternate_first_direction !< If true, alternate whether the x- or y-direction
!! updates occur first in directionally split parts of the calculation.
real :: first_dir_restart = -1.0 !< A real copy of G%first_direction for use in restart files
logical :: offline_tracer_mode = .false.
!< If true, step_offline() is called instead of step_MOM().
!! This is intended for running MOM6 in offline tracer mode
type(time_type), pointer :: Time !< pointer to the ocean clock
real :: dt !< (baroclinic) dynamics time step [T ~> s]
real :: dt_therm !< thermodynamics time step [T ~> s]
logical :: thermo_spans_coupling !< If true, thermodynamic and tracer time
!! steps can span multiple coupled time steps.
integer :: nstep_tot = 0 !< The total number of dynamic timesteps tcaaken
!! so far in this run segment
logical :: count_calls = .false. !< If true, count the calls to step_MOM, rather than the
!! number of dynamics steps in nstep_tot
logical :: debug !< If true, write verbose checksums for debugging purposes.
integer :: ntrunc !< number u,v truncations since last call to write_energy
integer :: cont_stencil !< The stencil for thickness from the continuity solver.
! These elements are used to control the dynamics updates.
logical :: do_dynamics !< If false, does not call step_MOM_dyn_*. This is an
!! undocumented run-time flag that is fragile.
logical :: split !< If true, use the split time stepping scheme.
logical :: use_RK2 !< If true, use RK2 instead of RK3 in unsplit mode
!! (i.e., no split between barotropic and baroclinic).
logical :: thickness_diffuse !< If true, diffuse interface height w/ a diffusivity KHTH.
logical :: thickness_diffuse_first !< If true, diffuse thickness before dynamics.
logical :: mixedlayer_restrat !< If true, use submesoscale mixed layer restratifying scheme.
logical :: useMEKE !< If true, call the MEKE parameterization.
logical :: useWaves !< If true, update Stokes drift
logical :: use_p_surf_in_EOS !< If true, always include the surface pressure contributions
!! in equation of state calculations.
logical :: use_diabatic_time_bug !< If true, uses the wrong calendar time for diabatic processes,
!! as was done in MOM6 versions prior to February 2018.
real :: dtbt_reset_period !< The time interval between dynamic recalculation of the
!! barotropic time step [s]. If this is negative dtbt is never
!! calculated, and if it is 0, dtbt is calculated every step.
type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period.
type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated.
real, dimension(:,:), pointer :: frac_shelf_h => NULL() !< fraction of total area occupied
!! by ice shelf [nondim]
real, dimension(:,:,:), pointer :: &
h_pre_dyn => NULL(), & !< The thickness before the transports [H ~> m or kg m-2].
T_pre_dyn => NULL(), & !< Temperature before the transports [degC].
S_pre_dyn => NULL() !< Salinity before the transports [ppt].
type(accel_diag_ptrs) :: ADp !< structure containing pointers to accelerations,
!! for derived diagnostics (e.g., energy budgets)
type(cont_diag_ptrs) :: CDp !< structure containing pointers to continuity equation
!! terms, for derived diagnostics (e.g., energy budgets)
real, dimension(:,:,:), pointer :: &
u_prev => NULL(), & !< previous value of u stored for diagnostics [L T-1 ~> m s-1]
v_prev => NULL() !< previous value of v stored for diagnostics [L T-1 ~> m s-1]
logical :: interp_p_surf !< If true, linearly interpolate surface pressure
!! over the coupling time step, using specified value
!! at the end of the coupling step. False by default.
logical :: p_surf_prev_set !< If true, p_surf_prev has been properly set from
!! a previous time-step or the ocean restart file.
!! This is only valid when interp_p_surf is true.
real, dimension(:,:), pointer :: &
p_surf_prev => NULL(), & !< surface pressure [R L2 T-2 ~> Pa] at end previous call to step_MOM
p_surf_begin => NULL(), & !< surface pressure [R L2 T-2 ~> Pa] at start of step_MOM_dyn_...
p_surf_end => NULL() !< surface pressure [R L2 T-2 ~> Pa] at end of step_MOM_dyn_...
! Variables needed to reach between start and finish phases of initialization
logical :: write_IC !< If true, then the initial conditions will be written to file
character(len=120) :: IC_file !< A file into which the initial conditions are
!! written in a new run if SAVE_INITIAL_CONDS is true.
logical :: calc_rho_for_sea_lev !< If true, calculate rho to convert pressure to sea level
! These elements are used to control the calculation and error checking of the surface state
real :: Hmix !< Diagnostic mixed layer thickness over which to
!! average surface tracer properties when a bulk
!! mixed layer is not used [Z ~> m], or a negative value
!! if a bulk mixed layer is being used.
real :: HFrz !< If HFrz > 0, the nominal depth over which melt potential is
!! computed [Z ~> m]. The actual depth over which melt potential is
!! computed is min(HFrz, OBLD), where OBLD is the boundary layer depth.
!! If HFrz <= 0 (default), melt potential will not be computed.
real :: Hmix_UV !< Depth scale over which to average surface flow to
!! feedback to the coupler/driver [Z ~> m] when
!! bulk mixed layer is not used, or a negative value
!! if a bulk mixed layer is being used.
logical :: check_bad_sfc_vals !< If true, scan surface state for ridiculous values.
real :: bad_val_ssh_max !< Maximum SSH before triggering bad value message [Z ~> m]
real :: bad_val_sst_max !< Maximum SST before triggering bad value message [degC]
real :: bad_val_sst_min !< Minimum SST before triggering bad value message [degC]
real :: bad_val_sss_max !< Maximum SSS before triggering bad value message [ppt]
real :: bad_val_col_thick !< Minimum column thickness before triggering bad value message [Z ~> m]
logical :: answers_2018 !< If true, use expressions for the surface properties that recover
!! the answers from the end of 2018. Otherwise, use more appropriate
!! expressions that differ at roundoff for non-Boussinsq cases.
logical :: use_particles !< Turns on the particles package
character(len=10) :: particle_type !< Particle types include: surface(default), profiling and sail drone.
type(MOM_diag_IDs) :: IDs !< Handles used for diagnostics.
type(transport_diag_IDs) :: transport_IDs !< Handles used for transport diagnostics.
type(surface_diag_IDs) :: sfc_IDs !< Handles used for surface diagnostics.
type(diag_grid_storage) :: diag_pre_sync !< The grid (thicknesses) before remapping
type(diag_grid_storage) :: diag_pre_dyn !< The grid (thicknesses) before dynamics
! The remainder of this type provides pointers to child module control structures.
type(MOM_dyn_unsplit_CS), pointer :: dyn_unsplit_CSp => NULL()
!< Pointer to the control structure used for the unsplit dynamics
type(MOM_dyn_unsplit_RK2_CS), pointer :: dyn_unsplit_RK2_CSp => NULL()
!< Pointer to the control structure used for the unsplit RK2 dynamics
type(MOM_dyn_split_RK2_CS), pointer :: dyn_split_RK2_CSp => NULL()
!< Pointer to the control structure used for the mode-split RK2 dynamics
type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp => NULL()
!< Pointer to the control structure used for the isopycnal height diffusive transport.
!! This is also common referred to as Gent-McWilliams diffusion
type(mixedlayer_restrat_CS), pointer :: mixedlayer_restrat_CSp => NULL()
!< Pointer to the control structure used for the mixed layer restratification
type(set_visc_CS), pointer :: set_visc_CSp => NULL()
!< Pointer to the control structure used to set viscosities
type(diabatic_CS), pointer :: diabatic_CSp => NULL()
!< Pointer to the control structure for the diabatic driver
type(MEKE_CS), pointer :: MEKE_CSp => NULL()
!< Pointer to the control structure for the MEKE updates
type(VarMix_CS), pointer :: VarMix => NULL()
!< Pointer to the control structure for the variable mixing module
type(Barotropic_CS), pointer :: Barotropic_CSp => NULL()
!< Pointer to the control structure for the barotropic module
type(tracer_registry_type), pointer :: tracer_Reg => NULL()
!< Pointer to the MOM tracer registry
type(tracer_advect_CS), pointer :: tracer_adv_CSp => NULL()
!< Pointer to the MOM tracer advection control structure
type(tracer_hor_diff_CS), pointer :: tracer_diff_CSp => NULL()
!< Pointer to the MOM along-isopycnal tracer diffusion control structure
type(tracer_flow_control_CS), pointer :: tracer_flow_CSp => NULL()
!< Pointer to the control structure that orchestrates the calling of tracer packages
! Although update_OBC_CS is not used directly outside of initialization, other modules
! set pointers to this type, so it should be kept for the duration of the run.
type(update_OBC_CS), pointer :: update_OBC_CSp => NULL()
!< Pointer to the control structure for updating open boundary condition properties
type(ocean_OBC_type), pointer :: OBC => NULL()
!< Pointer to the MOM open boundary condition type
type(sponge_CS), pointer :: sponge_CSp => NULL()
!< Pointer to the layered-mode sponge control structure
type(ALE_sponge_CS), pointer :: ALE_sponge_CSp => NULL()
!< Pointer to the oda incremental update control structure
type(oda_incupd_CS), pointer :: oda_incupd_CSp => NULL()
!< Pointer to the ALE-mode sponge control structure
type(ALE_CS), pointer :: ALE_CSp => NULL()
!< Pointer to the Arbitrary Lagrangian Eulerian (ALE) vertical coordinate control structure
! Pointers to control structures used for diagnostics
type(sum_output_CS), pointer :: sum_output_CSp => NULL()
!< Pointer to the globally summed output control structure
type(diagnostics_CS), pointer :: diagnostics_CSp => NULL()
!< Pointer to the MOM diagnostics control structure
type(offline_transport_CS), pointer :: offline_CSp => NULL()
!< Pointer to the offline tracer transport control structure
logical :: ensemble_ocean !< if true, this run is part of a
!! larger ensemble for the purpose of data assimilation
!! or statistical analysis.
type(ODA_CS), pointer :: odaCS => NULL() !< a pointer to the control structure for handling
!! ensemble model state vectors and data assimilation
!! increments and priors
type(particles), pointer :: particles => NULL() !<Lagrangian particles
end type MOM_control_struct
public initialize_MOM, finish_MOM_initialization, MOM_end
public step_MOM, step_offline
public extract_surface_state, get_ocean_stocks
public get_MOM_state_elements, MOM_state_is_synchronized
public allocate_surface_state, deallocate_surface_state
!>@{ CPU time clock IDs
integer :: id_clock_ocean
integer :: id_clock_dynamics
integer :: id_clock_thermo
integer :: id_clock_tracer
integer :: id_clock_diabatic
integer :: id_clock_adiabatic
integer :: id_clock_continuity ! also in dynamics s/r
integer :: id_clock_thick_diff
integer :: id_clock_BBL_visc
integer :: id_clock_ml_restrat
integer :: id_clock_diagnostics
integer :: id_clock_Z_diag
integer :: id_clock_init
integer :: id_clock_MOM_init
integer :: id_clock_pass ! also in dynamics d/r
integer :: id_clock_pass_init ! also in dynamics d/r
integer :: id_clock_ALE
integer :: id_clock_other
integer :: id_clock_offline_tracer
integer :: id_clock_unit_tests
!>@}
contains
!> This subroutine orchestrates the time stepping of MOM. The adiabatic
!! dynamics are stepped by calls to one of the step_MOM_dyn_...routines.
!! The action of lateral processes on tracers occur in calls to
!! advect_tracer and tracer_hordiff. Vertical mixing and possibly remapping
!! occur inside of diabatic.
subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS, &
Waves, do_dynamics, do_thermodynamics, start_cycle, &
end_cycle, cycle_length, reset_therm)
type(mech_forcing), target, intent(inout) :: forces_in !< A structure with the driving mechanical forces
type(forcing), target, intent(inout) :: fluxes_in !< A structure with pointers to themodynamic,
!! tracer and mass exchange forcing fields
type(surface), target, intent(inout) :: sfc_state !< surface ocean state
type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type
real, intent(in) :: time_int_in !< time interval covered by this run segment [s].
type(MOM_control_struct), intent(inout), target :: CS !< control structure from initialize_MOM
type(Wave_parameters_CS), &
optional, pointer :: Waves !< An optional pointer to a wave property CS
logical, optional, intent(in) :: do_dynamics !< Present and false, do not do updates due
!! to the dynamics.
logical, optional, intent(in) :: do_thermodynamics !< Present and false, do not do updates due
!! to the thermodynamics or remapping.
logical, optional, intent(in) :: start_cycle !< This indicates whether this call is to be
!! treated as the first call to step_MOM in a
!! time-stepping cycle; missing is like true.
logical, optional, intent(in) :: end_cycle !< This indicates whether this call is to be
!! treated as the last call to step_MOM in a
!! time-stepping cycle; missing is like true.
real, optional, intent(in) :: cycle_length !< The amount of time in a coupled time
!! stepping cycle [s].
logical, optional, intent(in) :: reset_therm !< This indicates whether the running sums of
!! thermodynamic quantities should be reset.
!! If missing, this is like start_cycle.
! local variables
type(ocean_grid_type), pointer :: G => NULL() ! pointer to a structure containing
! metrics and related information
type(ocean_grid_type), pointer :: G_in => NULL() ! Input grid metric
type(verticalGrid_type), pointer :: GV => NULL() ! Pointer to the vertical grid structure
type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing
! various unit conversion factors
integer :: ntstep ! time steps between tracer updates or diabatic forcing
integer :: n_max ! number of steps to take in this call
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
real :: time_interval ! time interval covered by this run segment [T ~> s].
real :: dt ! baroclinic time step [T ~> s]
real :: dtdia ! time step for diabatic processes [T ~> s]
real :: dt_therm ! a limited and quantized version of CS%dt_therm [T ~> s]
real :: dt_therm_here ! a further limited value of dt_therm [T ~> s]
real :: wt_end, wt_beg
real :: bbl_time_int ! The amount of time over which the calculated BBL
! properties will apply, for use in diagnostics, or 0
! if it is not to be calculated anew [T ~> s].
real :: rel_time = 0.0 ! relative time since start of this call [T ~> s].
logical :: calc_dtbt ! Indicates whether the dynamically adjusted
! barotropic time step needs to be updated.
logical :: do_advection ! If true, it is time to advect tracers.
logical :: do_calc_bbl ! If true, calculate the boundary layer properties.
logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans
! multiple dynamic timesteps.
logical :: do_dyn ! If true, dynamics are updated with this call.
logical :: do_thermo ! If true, thermodynamics and remapping may be applied with this call.
logical :: cycle_start ! If true, do calculations that are only done at the start of
! a stepping cycle (whatever that may mean).
logical :: cycle_end ! If true, do calculations and diagnostics that are only done at
! the end of a stepping cycle (whatever that may mean).
logical :: therm_reset ! If true, reset running sums of thermodynamic quantities.
real :: cycle_time ! The length of the coupled time-stepping cycle [T ~> s].
real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
ssh ! sea surface height, which may be based on eta_av [m]
real, dimension(:,:,:), pointer :: &
u => NULL(), & ! u : zonal velocity component [L T-1 ~> m s-1]
v => NULL(), & ! v : meridional velocity component [L T-1 ~> m s-1]
h => NULL() ! h : layer thickness [H ~> m or kg m-2]
real, dimension(:,:), pointer :: &
p_surf => NULL() ! A pointer to the ocean surface pressure [R L2 T-2 ~> Pa].
real :: I_wt_ssh ! The inverse of the time weights [T-1 ~> s-1]
type(time_type) :: Time_local, end_time_thermo, Time_temp
type(group_pass_type) :: pass_tau_ustar_psurf
logical :: showCallTree
! External forcing fields on the model index map
type(mech_forcing), pointer :: forces ! Mechanical forcing
type(forcing), pointer :: fluxes ! Boundary fluxes
type(surface), pointer :: sfc_state_diag ! Surface boundary fields
integer :: turns ! Number of quarter turns from input to model indexing
G => CS%G ; G_in => CS%G_in ; GV => CS%GV ; US => CS%US
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
u => CS%u ; v => CS%v ; h => CS%h
time_interval = US%s_to_T*time_int_in
do_dyn = .true. ; if (present(do_dynamics)) do_dyn = do_dynamics
do_thermo = .true. ; if (present(do_thermodynamics)) do_thermo = do_thermodynamics
if (.not.(do_dyn .or. do_thermo)) call MOM_error(FATAL,"Step_MOM: "//&
"Both do_dynamics and do_thermodynamics are false, which makes no sense.")
cycle_start = .true. ; if (present(start_cycle)) cycle_start = start_cycle
cycle_end = .true. ; if (present(end_cycle)) cycle_end = end_cycle
cycle_time = time_interval ; if (present(cycle_length)) cycle_time = US%s_to_T*cycle_length
therm_reset = cycle_start ; if (present(reset_therm)) therm_reset = reset_therm
call cpu_clock_begin(id_clock_ocean)
call cpu_clock_begin(id_clock_other)
if (CS%debug) then
call MOM_state_chksum("Beginning of step_MOM ", u, v, h, CS%uh, CS%vh, G, GV, US)
endif
showCallTree = callTree_showQuery()
if (showCallTree) call callTree_enter("step_MOM(), MOM.F90")
! Rotate the forces from G_in to G
if (CS%rotate_index) then
turns = G%HI%turns
allocate(forces)
call allocate_mech_forcing(forces_in, G, forces)
call rotate_mech_forcing(forces_in, turns, forces)
allocate(fluxes)
call allocate_forcing_type(fluxes_in, G, fluxes)
call rotate_forcing(fluxes_in, fluxes, turns)
else
forces => forces_in
fluxes => fluxes_in
endif
! First determine the time step that is consistent with this call and an
! integer fraction of time_interval.
if (do_dyn) then
n_max = 1
if (time_interval > CS%dt) n_max = ceiling(time_interval/CS%dt - 0.001)
dt = time_interval / real(n_max)
thermo_does_span_coupling = (CS%thermo_spans_coupling .and. &
(CS%dt_therm > 1.5*cycle_time))
if (thermo_does_span_coupling) then
! Set dt_therm to be an integer multiple of the coupling time step.
dt_therm = cycle_time * floor(CS%dt_therm / cycle_time + 0.001)
ntstep = floor(dt_therm/dt + 0.001)
elseif (.not.do_thermo) then
dt_therm = CS%dt_therm
if (present(cycle_length)) dt_therm = min(CS%dt_therm, US%s_to_T*cycle_length)
! ntstep is not used.
else
ntstep = MAX(1, MIN(n_max, floor(CS%dt_therm/dt + 0.001)))
dt_therm = dt*ntstep
endif
if (associated(forces%p_surf)) p_surf => forces%p_surf
if (.not.associated(forces%p_surf)) CS%interp_p_surf = .false.
CS%tv%p_surf => NULL()
if (CS%use_p_surf_in_EOS .and. associated(forces%p_surf)) CS%tv%p_surf => forces%p_surf
!---------- Initiate group halo pass of the forcing fields
call cpu_clock_begin(id_clock_pass)
if (.not.associated(forces%taux) .or. .not.associated(forces%tauy)) &
call MOM_error(FATAL,'step_MOM:forces%taux,tauy not associated')
call create_group_pass(pass_tau_ustar_psurf, forces%taux, forces%tauy, G%Domain)
if (associated(forces%ustar)) &
call create_group_pass(pass_tau_ustar_psurf, forces%ustar, G%Domain)
if (associated(forces%p_surf)) &
call create_group_pass(pass_tau_ustar_psurf, forces%p_surf, G%Domain)
if (G%nonblocking_updates) then
call start_group_pass(pass_tau_ustar_psurf, G%Domain)
else
call do_group_pass(pass_tau_ustar_psurf, G%Domain)
endif
call cpu_clock_end(id_clock_pass)
else
! This step only updates the thermodynamics so setting timesteps is simpler.
n_max = 1
if ((time_interval > CS%dt_therm) .and. (CS%dt_therm > 0.0)) &
n_max = ceiling(time_interval/CS%dt_therm - 0.001)
dt = time_interval / real(n_max)
dt_therm = dt ; ntstep = 1
if (associated(fluxes%p_surf)) p_surf => fluxes%p_surf
CS%tv%p_surf => NULL()
if (associated(fluxes%p_surf)) then
if (CS%use_p_surf_in_EOS) CS%tv%p_surf => fluxes%p_surf
endif
if (CS%UseWaves) call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass)
endif
if (therm_reset) then
CS%time_in_thermo_cycle = 0.0
if (associated(CS%tv%frazil)) CS%tv%frazil(:,:) = 0.0
if (associated(CS%tv%salt_deficit)) CS%tv%salt_deficit(:,:) = 0.0
if (associated(CS%tv%TempxPmE)) CS%tv%TempxPmE(:,:) = 0.0
if (associated(CS%tv%internal_heat)) CS%tv%internal_heat(:,:) = 0.0
endif
if (cycle_start) then
CS%time_in_cycle = 0.0
do j=js,je ; do i=is,ie ; CS%ssh_rint(i,j) = 0.0 ; enddo ; enddo
if (associated(CS%VarMix)) then
call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag)
call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix)
call calc_depth_function(G, CS%VarMix)
call disable_averaging(CS%diag)
endif
endif
if (do_dyn) then
if (G%nonblocking_updates) &
call complete_group_pass(pass_tau_ustar_psurf, G%Domain, clock=id_clock_pass)
if (CS%interp_p_surf) then
if (.not.associated(CS%p_surf_end)) allocate(CS%p_surf_end(isd:ied,jsd:jed))
if (.not.associated(CS%p_surf_begin)) allocate(CS%p_surf_begin(isd:ied,jsd:jed))
if (.not.CS%p_surf_prev_set) then
do j=jsd,jed ; do i=isd,ied
CS%p_surf_prev(i,j) = forces%p_surf(i,j)
enddo ; enddo
CS%p_surf_prev_set = .true.
endif
else
CS%p_surf_end => forces%p_surf
endif
if (CS%UseWaves) then
! Update wave information, which is presently kept static over each call to step_mom
call enable_averages(time_interval, Time_start + real_to_time(US%T_to_s*time_interval), CS%diag)
call Update_Stokes_Drift(G, GV, US, Waves, h, forces%ustar)
call disable_averaging(CS%diag)
endif
else ! not do_dyn.
if (CS%UseWaves) & ! Diagnostics are not enabled in this call.
call Update_Stokes_Drift(G, GV, US, Waves, h, fluxes%ustar)
endif
if (CS%debug) then
if (cycle_start) &
call MOM_state_chksum("Before steps ", u, v, h, CS%uh, CS%vh, G, GV, US)
if (cycle_start) call check_redundant("Before steps ", u, v, G)
if (do_dyn) call MOM_mech_forcing_chksum("Before steps", forces, G, US, haloshift=0)
if (do_dyn) call check_redundant("Before steps ", forces%taux, forces%tauy, G)
endif
call cpu_clock_end(id_clock_other)
rel_time = 0.0
do n=1,n_max
if (CS%use_diabatic_time_bug) then
! This wrong form of update was used until Feb 2018, recovered with CS%use_diabatic_time_bug=T.
CS%Time = Time_start + real_to_time(US%T_to_s*int(floor(rel_time+0.5*dt+0.5)))
rel_time = rel_time + dt
else
rel_time = rel_time + dt ! The relative time at the end of the step.
! Set the universally visible time to the middle of the time step.
CS%Time = Time_start + real_to_time(US%T_to_s*(rel_time - 0.5*dt))
endif
! Set the local time to the end of the time step.
Time_local = Time_start + real_to_time(US%T_to_s*rel_time)
if (showCallTree) call callTree_enter("DT cycles (step_MOM) n=",n)
! Update the vertically extensive diagnostic grids so that they are
! referenced to the beginning timestep
call diag_update_remap_grids(CS%diag, update_intensive = .false., update_extensive = .true. )
!===========================================================================
! This is the first place where the diabatic processes and remapping could occur.
if (CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0) .and. do_thermo) then ! do thermodynamics.
if (.not.do_dyn) then
dtdia = dt
elseif (thermo_does_span_coupling) then
dtdia = dt_therm
if ((fluxes%dt_buoy_accum > 0.0) .and. (dtdia > time_interval) .and. &
(abs(fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia)) then
call MOM_error(FATAL, "step_MOM: Mismatch between long thermodynamic "//&
"timestep and time over which buoyancy fluxes have been accumulated.")
endif
call MOM_error(FATAL, "MOM is not yet set up to have restarts that work "//&
"with THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
else
dtdia = dt*min(ntstep,n_max-(n-1))
endif
end_time_thermo = Time_local
if (dtdia > dt .and. .not. CS%use_diabatic_time_bug) then
! If necessary, temporarily reset CS%Time to the center of the period covered
! by the call to step_MOM_thermo, noting that they begin at the same time.
! This step was missing prior to Feb 2018, and is skipped with CS%use_diabatic_time_bug=T.
CS%Time = CS%Time + real_to_time(0.5*US%T_to_s*(dtdia-dt))
endif
if (dtdia > dt .or. CS%use_diabatic_time_bug) then
! The end-time of the diagnostic interval needs to be set ahead if there
! are multiple dynamic time steps worth of thermodynamics applied here.
! This line was not conditional prior to Feb 2018, recovered with CS%use_diabatic_time_bug=T.
end_time_thermo = Time_local + real_to_time(US%T_to_s*(dtdia-dt))
endif
! Apply diabatic forcing, do mixing, and regrid.
call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, &
end_time_thermo, .true., Waves=Waves)
CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia
! The diabatic processes are now ahead of the dynamics by dtdia.
CS%t_dyn_rel_thermo = -dtdia
if (showCallTree) call callTree_waypoint("finished diabatic_first (step_MOM)")
if (dtdia > dt .and. .not. CS%use_diabatic_time_bug) & ! Reset CS%Time to its previous value.
! This step was missing prior to Feb 2018, recovered with CS%use_diabatic_time_bug=T.
CS%Time = Time_start + real_to_time(US%T_to_s*(rel_time - 0.5*dt))
endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))"
if (do_dyn) then
! Store pre-dynamics thicknesses for proper diagnostic remapping for transports or
! advective tendencies. If there are more than one dynamics steps per advective
! step (i.e DT_THERM > DT), this needs to be stored at the first dynamics call.
if (.not.CS%preadv_h_stored .and. (CS%t_dyn_rel_adv == 0.)) then
call diag_copy_diag_to_storage(CS%diag_pre_dyn, h, CS%diag)
CS%preadv_h_stored = .true.
endif
! The pre-dynamics velocities might be stored for debugging truncations.
if (associated(CS%u_prev) .and. associated(CS%v_prev)) then
do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB
CS%u_prev(I,j,k) = u(I,j,k)
enddo ; enddo ; enddo
do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied
CS%v_prev(I,j,k) = v(i,J,k)
enddo ; enddo ; enddo
endif
dt_therm_here = dt_therm
if (do_thermo .and. do_dyn .and. .not.thermo_does_span_coupling) &
dt_therm_here = dt*min(ntstep, n_max-n+1)
! Indicate whether the bottom boundary layer properties need to be
! recalculated, and if so for how long an interval they are valid.
bbl_time_int = 0.0
if (do_thermo) then
if ((CS%t_dyn_rel_adv == 0.0) .or. (n==1)) &
bbl_time_int = max(dt, min(dt_therm - CS%t_dyn_rel_adv, dt*(1+n_max-n)) )
else
if ((CS%t_dyn_rel_adv == 0.0) .or. ((n==1) .and. cycle_start)) &
bbl_time_int = min(dt_therm, cycle_time)
endif
if (CS%interp_p_surf) then
wt_end = real(n) / real(n_max)
wt_beg = real(n-1) / real(n_max)
do j=jsd,jed ; do i=isd,ied
CS%p_surf_end(i,j) = wt_end * forces%p_surf(i,j) + &
(1.0-wt_end) * CS%p_surf_prev(i,j)
CS%p_surf_begin(i,j) = wt_beg * forces%p_surf(i,j) + &
(1.0-wt_beg) * CS%p_surf_prev(i,j)
enddo ; enddo
endif
call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, &
dt_therm_here, bbl_time_int, CS, &
Time_local, Waves=Waves)
!===========================================================================
! This is the start of the tracer advection part of the algorithm.
if (thermo_does_span_coupling .or. .not.do_thermo) then
do_advection = (CS%t_dyn_rel_adv + 0.5*dt > dt_therm)
else
do_advection = ((MOD(n,ntstep) == 0) .or. (n==n_max))
endif
if (do_advection) then ! Do advective transport and lateral tracer mixing.
call step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local)
if (CS%diabatic_first .and. abs(CS%t_dyn_rel_thermo) > 1e-6*dt) call MOM_error(FATAL, &
"step_MOM: Mismatch between the dynamics and diabatic times "//&
"with DIABATIC_FIRST.")
endif
endif ! end of (do_dyn)
!===========================================================================
! This is the second place where the diabatic processes and remapping could occur.
if ((CS%t_dyn_rel_adv==0.0) .and. do_thermo .and. (.not.CS%diabatic_first)) then
dtdia = CS%t_dyn_rel_thermo
! If the MOM6 dynamic and thermodynamic time stepping is being orchestrated
! by the coupler, the value of diabatic_first does not matter.
if ((CS%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) dtdia = dt
if (CS%thermo_spans_coupling .and. (CS%dt_therm > 1.5*cycle_time) .and. &
(abs(dt_therm - dtdia) > 1e-6*dt_therm)) then
call MOM_error(FATAL, "step_MOM: Mismatch between dt_therm and dtdia "//&
"before call to diabatic.")
endif
! If necessary, temporarily reset CS%Time to the center of the period covered
! by the call to step_MOM_thermo, noting that they end at the same time.
! This step was missing prior to Feb 2018, and is skipped with CS%use_diabatic_time_bug=T.
if (dtdia > dt .and. .not. CS%use_diabatic_time_bug) &
CS%Time = CS%Time - real_to_time(0.5*US%T_to_s*(dtdia-dt))
! Apply diabatic forcing, do mixing, and regrid.
call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, &
Time_local, .false., Waves=Waves)
CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia
if ((CS%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) then
! The diabatic processes are now ahead of the dynamics by dtdia.
CS%t_dyn_rel_thermo = -dtdia
else ! The diabatic processes and the dynamics are synchronized.
CS%t_dyn_rel_thermo = 0.0
endif
! Reset CS%Time to its previous value.
! This step was missing prior to Feb 2018, and is skipped with CS%use_diabatic_time_bug=T.
if (dtdia > dt .and. .not. CS%use_diabatic_time_bug) &
CS%Time = Time_start + real_to_time(US%T_to_s*(rel_time - 0.5*dt))
endif
if (do_dyn) then
call cpu_clock_begin(id_clock_dynamics)
! Determining the time-average sea surface height is part of the algorithm.
! This may be eta_av if Boussinesq, or need to be diagnosed if not.
CS%time_in_cycle = CS%time_in_cycle + dt
call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, eta_to_m=1.0, dZref=G%Z_ref)
do j=js,je ; do i=is,ie
CS%ssh_rint(i,j) = CS%ssh_rint(i,j) + dt*ssh(i,j)
enddo ; enddo
if (CS%IDs%id_ssh_inst > 0) call post_data(CS%IDs%id_ssh_inst, ssh, CS%diag)
call cpu_clock_end(id_clock_dynamics)
endif
!===========================================================================
! Calculate diagnostics at the end of the time step if the state is self-consistent.
if (MOM_state_is_synchronized(CS)) then
!### Perhaps this should be if (CS%t_dyn_rel_thermo == 0.0)
call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
! Diagnostics that require the complete state to be up-to-date can be calculated.
call enable_averages(CS%t_dyn_rel_diag, Time_local, CS%diag)
call calculate_diagnostic_fields(u, v, h, CS%uh, CS%vh, CS%tv, CS%ADp, &
CS%CDp, p_surf, CS%t_dyn_rel_diag, CS%diag_pre_sync,&
G, GV, US, CS%diagnostics_CSp)
call post_tracer_diagnostics_at_sync(CS%Tracer_reg, h, CS%diag_pre_sync, CS%diag, G, GV, CS%t_dyn_rel_diag)
call diag_copy_diag_to_storage(CS%diag_pre_sync, h, CS%diag)
if (showCallTree) call callTree_waypoint("finished calculate_diagnostic_fields (step_MOM)")
call disable_averaging(CS%diag)
CS%t_dyn_rel_diag = 0.0
call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
endif
if (do_dyn .and. .not.CS%count_calls) CS%nstep_tot = CS%nstep_tot + 1
if (showCallTree) call callTree_leave("DT cycles (step_MOM)")
enddo ! complete the n loop
if (CS%count_calls .and. cycle_start) CS%nstep_tot = CS%nstep_tot + 1
call cpu_clock_begin(id_clock_other)
if (CS%time_in_cycle > 0.0) then
I_wt_ssh = 1.0/CS%time_in_cycle
do j=js,je ; do i=is,ie
ssh(i,j) = CS%ssh_rint(i,j)*I_wt_ssh
CS%ave_ssh_ibc(i,j) = ssh(i,j)
enddo ; enddo
if (do_dyn) then
call adjust_ssh_for_p_atm(CS%tv, G, GV, US, CS%ave_ssh_ibc, forces%p_surf_SSH, &
CS%calc_rho_for_sea_lev)
elseif (do_thermo) then
call adjust_ssh_for_p_atm(CS%tv, G, GV, US, CS%ave_ssh_ibc, fluxes%p_surf_SSH, &
CS%calc_rho_for_sea_lev)
endif
endif
if (do_dyn .and. CS%interp_p_surf) then ; do j=jsd,jed ; do i=isd,ied
CS%p_surf_prev(i,j) = forces%p_surf(i,j)
enddo ; enddo ; endif
if (CS%ensemble_ocean) then
! store ensemble vector in odaCS
call set_prior_tracer(CS%Time, G, GV, CS%h, CS%tv, CS%odaCS)
! call DA interface
call oda(CS%Time,CS%odaCS)
! update the time for the next analysis step if needed
call set_analysis_time(CS%Time,CS%odaCS)
endif
if (showCallTree) call callTree_waypoint("calling extract_surface_state (step_MOM)")
! NOTE: sfc_state uses input indexing, since it is also used by drivers.
call extract_surface_state(CS, sfc_state)
! Do diagnostics that only occur at the end of a complete forcing step.
if (cycle_end) then
if (CS%rotate_index) then
allocate(sfc_state_diag)
call rotate_surface_state(sfc_state, G_in, sfc_state_diag, G, turns)
else
sfc_state_diag => sfc_state
endif
call cpu_clock_begin(id_clock_diagnostics)
if (CS%time_in_cycle > 0.0) then
call enable_averages(CS%time_in_cycle, Time_local, CS%diag)
call post_surface_dyn_diags(CS%sfc_IDs, G, CS%diag, sfc_state_diag, ssh)
endif
if (CS%time_in_thermo_cycle > 0.0) then
call enable_averages(CS%time_in_thermo_cycle, Time_local, CS%diag)
call post_surface_thermo_diags(CS%sfc_IDs, G, GV, US, CS%diag, CS%time_in_thermo_cycle, &
sfc_state_diag, CS%tv, ssh, CS%ave_ssh_ibc)
endif
call disable_averaging(CS%diag)
call cpu_clock_end(id_clock_diagnostics)
endif
! Accumulate the surface fluxes for assessing conservation
if (do_thermo .and. fluxes%fluxes_used) &
call accumulate_net_input(fluxes, sfc_state, CS%tv, fluxes%dt_buoy_accum, &
G, US, CS%sum_output_CSp)
if (MOM_state_is_synchronized(CS)) &
call write_energy(CS%u, CS%v, CS%h, CS%tv, Time_local, CS%nstep_tot, &
G, GV, US, CS%sum_output_CSp, CS%tracer_flow_CSp, &
dt_forcing=real_to_time(US%T_to_s*time_interval) )
call cpu_clock_end(id_clock_other)
! De-rotate fluxes and copy back to the input, since they can be changed.
if (CS%rotate_index) then
call rotate_forcing(fluxes, fluxes_in, -turns)
call rotate_mech_forcing(forces, -turns, forces_in)
call deallocate_mech_forcing(forces)
deallocate(forces)
call deallocate_forcing_type(fluxes)
deallocate(fluxes)
endif
if (showCallTree) call callTree_leave("step_MOM()")
call cpu_clock_end(id_clock_ocean)
end subroutine step_MOM
!> Time step the ocean dynamics, including the momentum and continuity equations
subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
bbl_time_int, CS, Time_local, Waves)
type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
real, dimension(:,:), pointer :: p_surf_begin !< A pointer (perhaps NULL) to the surface
!! pressure at the beginning of this dynamic
!! step, intent in [R L2 T-2 ~> Pa].
real, dimension(:,:), pointer :: p_surf_end !< A pointer (perhaps NULL) to the surface
!! pressure at the end of this dynamic step,
!! intent in [R L2 T-2 ~> Pa].
real, intent(in) :: dt !< time interval covered by this call [T ~> s].
real, intent(in) :: dt_thermo !< time interval covered by any updates that may
!! span multiple dynamics steps [T ~> s].
real, intent(in) :: bbl_time_int !< time interval over which updates to the
!! bottom boundary layer properties will apply [T ~> s],
!! or zero not to update the properties.
type(MOM_control_struct), intent(inout), target :: CS !< control structure from initialize_MOM
type(time_type), intent(in) :: Time_local !< End time of a segment, as a time type
type(wave_parameters_CS), &