-
Notifications
You must be signed in to change notification settings - Fork 96
/
ocean_tracer.F90
4998 lines (4206 loc) · 220 KB
/
ocean_tracer.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
module ocean_tracer_mod
#define COMP isc:iec,jsc:jec
!
!<CONTACT EMAIL="[email protected]"> Stephen Griffies
!</CONTACT>
!
!<CONTACT EMAIL="[email protected]"> Matt Harrison
!</CONTACT>
!
!<CONTACT EMAIL="[email protected]"> Richard D. Slater (initialization)
!</CONTACT>
!
!<OVERVIEW>
! This module time steps the tracer fields.
!</OVERVIEW>
!
!<DESCRIPTION>
! This module time steps the tracer fields.
! Initialization for the tracer packages is done as well.
!</DESCRIPTION>
!
! <INFO>
!
! <REFERENCE>
! R.C. Pacanowski and S.M. Griffies
! The MOM3 Manual (1999)
! </REFERENCE>
!
! <REFERENCE>
! S.M. Griffies, M.J. Harrison, R.C. Pacanowski, and A. Rosati
! A Technical Guide to MOM4 (2004)
! </REFERENCE>
!
! <REFERENCE>
! S.M. Griffies, Fundamentals of ocean climate models (2004)
! </REFERENCE>
!
! </INFO>
!
!<NAMELIST NAME="ocean_tracer_nml">
!
! <DATA NAME="zero_tendency" TYPE="logical">
! If true, then will freeze the tracer fields.
! </DATA>
!
! <DATA NAME="zero_tracer_source" TYPE="logical">
! To remove the T_prog%source contribution to tracer
! evolution. For debugging purposes. Default
! zero_tracer_source=.false.
! </DATA>
!
! <DATA NAME="limit_age_tracer" TYPE="logical">
! Limit the values of age tracer to be less than
! total run time and greater than zero.
! Default limit_age_tracer=.false.
! </DATA>
!
! <DATA NAME="age_tracer_max_init" TYPE="real" UNITS="years">
! Initial maximum age tracer. This nml provides the ability to
! start an integration with an age tracer that is not initialized
! to zero, say if we took an age tracer from another spin-up.
! Default age_tracer_max_init=0.0.
! </DATA>
!
! <DATA NAME="remap_depth_to_s_init" TYPE="logical">
! For remapping initial tracer distributions, generally determined
! according to depth vertical coordinates using the mom preprocessing
! schemes, onto s-coordinates. This method is of use for initializing
! terrain following coordinate simulations with mom.
! </DATA>
!
! <DATA NAME="frazil_heating_before_vphysics" TYPE="logical">
! For computing frazil heating before the implicit vertical physics
! (which includes boundary fluxes), and before vertical convection.
! This is the order that CM2.0 and CM2.1 performed their calculations
! of frazil. It is arguable that one should NOT do frazil until the
! end of a time step, after vertical physics and after surface
! boundary fluxes.
! Default frazil_heating_before_vphysics=.false.
! </DATA>
!
! <DATA NAME="frazil_heating_after_vphysics" TYPE="logical">
! For computing frazil heating after the implicit vertical physics
! (which includes boundary fluxes), and after vertical convection.
! This is the recommended method.
! Default frazil_heating_after_vphysics=.false.
! </DATA>
!
! <DATA NAME="tmask_limit_ts_same" TYPE="logical">
! tmask_limit is derived separately for the tracers. However,
! it may be appropriate to have the mask be the same for temp
! and salinity, in which case the neutral physics fluxes are
! self-consistent. But for some cases, such as when running with
! linear eos, may not wish to have the temp and salinity coupled
! when computing the mask.
! </DATA>
!
! <DATA NAME="compute_tmask_limit_on" TYPE="logical">
! For updating the tmaks_limit array. This calculation is
! recommended for the following physics and advection schemes:
! 1/ quicker advection
! 2/ neutral physics
! 3/ submesoscale closure.
! The default is compute_tmask_limit_on=.true., but if none
! of the above schemes is used, then some time savings can be
! realized by setting compute_tmask_limit_on=.false.
! </DATA>
!
! <DATA NAME="debug_this_module" TYPE="logical">
! For debugging the tracer module
! </DATA>
! <DATA NAME="write_a_restart" TYPE="logical">
! Set true to write a restart. False setting only for rare
! cases where wish to benchmark model without measuring the cost
! of writing restarts and associated chksums.
! Default is write_a_restart=.true.
! </DATA>
!
! <DATA NAME="interpolate_tprog_to_pbott" TYPE="logical">
! To linear interpolate the initial conditions for prognostic
! tracers to the partial bottom cells. Default
! interpolate_tprog_to_pbott=.true.
! </DATA>
!
! <DATA NAME="interpolate_tdiag_to_pbott" TYPE="logical">
! To linear interpolate the initial conditions for diagnostic
! tracers to the partial bottom cells. Default
! interpolate_tdiag_to_pbott=.false.
! </DATA>
!
! <DATA NAME="inflow_nboundary" TYPE="logical">
! For adding an inflow transport from the northern boundary
! which brings in temp and salinity according to inflow data
! files. Default is inflow_nboundary=.false.
! </DATA>
!
! <DATA NAME="ocean_tpm_debug" TYPE="logical">
! For debugging ocean tracer package manager.
! </DATA>
!
! <DATA NAME="use_tempsalt_check_range" TYPE="logical">
! To call a check to see that temperature and salinity
! are within their pre-selected range.
! Default use_tempsalt_check_range=.false. since this
! check may incur some cost that users should be aware of.
! </DATA>
!</NAMELIST>
!
use constants_mod, only: epsln, kelvin
use diag_manager_mod, only: register_diag_field, send_data
use field_manager_mod, only: fm_string_len, fm_type_name_len
use field_manager_mod, only: fm_get_length
use field_manager_mod, only: fm_dump_list, fm_change_list, fm_loop_over_list
use fms_mod, only: read_data, file_exist, field_exist
use fms_mod, only: open_namelist_file, check_nml_error, close_file
use fms_mod, only: FATAL, WARNING, NOTE, stdout, stdlog
use fms_io_mod, only: register_restart_field, save_restart, restore_state
use fms_io_mod, only: restart_file_type, reset_field_pointer, reset_field_name
use fms_io_mod, only: field_size
use fm_util_mod, only: fm_util_get_real, fm_util_get_integer
use fm_util_mod, only: fm_util_get_logical, fm_util_get_string
use fm_util_mod, only: fm_util_get_string_array
use fm_util_mod, only: fm_util_check_for_bad_fields, fm_util_set_value
use mpp_domains_mod, only: mpp_update_domains
use mpp_domains_mod, only: mpp_global_sum, NON_BITWISE_EXACT_SUM
use mpp_io_mod, only: mpp_open, fieldtype
use mpp_io_mod, only: MPP_NETCDF, MPP_OVERWR, MPP_ASCII, MPP_RDONLY, MPP_SINGLE, MPP_MULTI
use mpp_mod, only: input_nml_file, mpp_error, mpp_pe, mpp_root_pe, mpp_broadcast, ALL_PES, mpp_max, mpp_min, mpp_chksum
use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end
use mpp_mod, only: CLOCK_COMPONENT, CLOCK_SUBCOMPONENT, CLOCK_MODULE, CLOCK_ROUTINE
use platform_mod, only: i8_kind
use time_manager_mod, only: time_type, set_time, increment_time, operator( + )
use transport_matrix_mod, only: do_transport_matrix, transport_matrix_store_explicit
use ocean_blob_mod, only: ocean_blob_implicit, adjust_L_thickness
use ocean_convect_mod, only: convection
use ocean_domains_mod, only: get_local_indices
use ocean_density_mod, only: neutral_density, update_ocean_density_salinity
use ocean_frazil_mod, only: compute_frazil_heating
use ocean_bih_tracer_mod, only: bih_tracer
#ifdef USE_OCEAN_BGC
use ocean_generic_mod, only: ocean_generic_get_field, ocean_generic_get_field_pointer
use ocean_generic_mod, only: ocean_generic_set_pointer
#endif
use ocean_lap_tracer_mod, only: lap_tracer
use ocean_obc_mod, only: ocean_obc_tracer, ocean_obc_update_boundary, ocean_obc_tracer_init
use ocean_parameters_mod, only: ADVECT_UPWIND, ADVECT_2ND_ORDER, ADVECT_4TH_ORDER, ADVECT_6TH_ORDER
use ocean_parameters_mod, only: ADVECT_QUICKER, ADVECT_QUICKMOM3, ADVECT_MDFL_SUP_B, ADVECT_MDFL_SWEBY
use ocean_parameters_mod, only: ADVECT_PSOM, ADVECT_MDPPM, ADVECT_DST_LINEAR, ADVECT_MDMDT_TEST
use ocean_parameters_mod, only: ADVECT_MDPPM_TEST, ADVECT_MDFL_SWEBY_TEST, ADVECT_DST_LINEAR_TEST
use ocean_parameters_mod, only: missing_value, sec_in_yr_r, rho0, rho0r, cp_ocean
use ocean_parameters_mod, only: TWO_LEVEL, THREE_LEVEL
use ocean_parameters_mod, only: CONSERVATIVE_TEMP, POTENTIAL_TEMP
use ocean_parameters_mod, only: QUASI_HORIZONTAL, TERRAIN_FOLLOWING
use ocean_parameters_mod, only: GEOPOTENTIAL
use ocean_parameters_mod, only: DEPTH_BASED
use ocean_passive_mod, only: passive_tracer_init, update_tracer_passive
use ocean_shortwave_mod, only: ocean_irradiance_init
use ocean_tempsalt_mod, only: contemp_from_pottemp, pottemp_from_contemp, tempsalt_check_range
use ocean_tpm_mod, only: ocean_tpm_init
use ocean_tpm_util_mod, only: otpm_set_tracer_package
use ocean_tracer_advect_mod, only: horz_advect_tracer, vert_advect_tracer
use ocean_tracer_diag_mod, only: send_tracer_variance
use ocean_tracer_util_mod, only: rebin_onto_rho, diagnose_mass_of_layer
use ocean_tracer_util_mod, only: tracer_prog_chksum, tracer_diag_chksum, tracer_min_max
use ocean_thickness_mod, only: update_E_thickness
use ocean_types_mod, only: ocean_grid_type, ocean_domain_type, ocean_thickness_type
use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_diag_tracer_type
use ocean_types_mod, only: ocean_adv_vel_type, ocean_external_mode_type
use ocean_types_mod, only: ocean_public_type, ocean_density_type, ocean_options_type
use ocean_types_mod, only: ocean_lagrangian_type, ocean_velocity_type, blob_diag_type
use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, diagnose_sum, write_chksum_3d
use ocean_tracer_util_mod, only: diagnose_3d_rho
use ocean_vert_mix_mod, only: vert_diffuse, vert_diffuse_implicit
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk6, wrk1_2d
implicit none
private
logical :: prog_module_initialized = .false.
logical :: diag_module_initialized = .false.
character(len=256) :: version='CVS $Id: ocean_tracer.F90,v 20.0 2013/12/14 00:17:20 fms Exp $'
character(len=256) :: tagname='Tag $Name: tikal $'
character(len=48), parameter :: mod_name = 'ocean_tracer_mod'
integer :: num_tracers =0
integer :: num_prog_tracers =0
integer :: num_diag_tracers =0
integer :: num_family_tracers=0
#include <ocean_memory.h>
type(ocean_grid_type), pointer :: Grd =>NULL()
type(ocean_domain_type), pointer :: Dom =>NULL()
! for time steps and implicit vertical mixing
integer :: tendency = 0
real :: dtime = 0.0
real :: dtimer = 0.0
real :: dtts = 0.0
real :: dtuv = 0.0
real :: cp_oceanr
integer :: index_temp =-1
integer :: index_salt =-1
integer :: index_frazil =-1
integer :: index_temp_sq =-1
integer :: index_salt_sq =-1
integer :: index_diag_temp =-1
integer :: index_added_heat =-1
integer :: index_redist_heat =-1
! for obc
logical :: have_obc=.false.
! for global normalization
real :: cellarea_r
! for vertical coordinate class
integer :: vert_coordinate_class
! for setting the temperature variable
integer :: prog_temp_variable=0
! for CMIP units (K rather than degrees C if cmip_version<=5)
real :: cmip_offset = 0.0
character(len=32) :: temp_units='degrees C'
! for possible inflow tracer settings
real, dimension(:,:,:), allocatable :: mask_inflow
real, dimension(:,:,:), allocatable :: salt_inflow
real, dimension(:,:,:), allocatable :: temp_inflow
! identification numbers for mpp clocks
integer :: id_clock_bih_tracer
integer :: id_clock_lap_tracer
integer :: id_clock_vert_diffuse
integer :: id_clock_vert_diffuse_implicit
integer :: id_clock_tracer_advect
integer :: id_clock_frazil
integer :: id_clock_convection
integer :: id_clock_blob
integer :: id_clock_blob_implicit
integer :: id_clock_adjust_L_thickness
integer :: id_clock_adjust_E_thickness
! for restart
integer :: num_tracer_restart = 0
integer, allocatable :: id_tracer_restart(:)
integer, allocatable :: id_tracer_type(:)
character(len=64), allocatable :: tracer_restart_file(:)
type(restart_file_type), allocatable :: Tracer_restart(:)
!work array on neutral density space
integer :: neutralrho_nk
real, dimension(:,:,:), allocatable :: nrho_work
real, dimension(:,:,:), allocatable :: nrho_work2
real, dimension(:,:,:), allocatable :: nrho_mass
! internally set for computing watermass diagnstics
logical :: compute_watermass_diag = .false.
! for diagnostics
logical :: used
integer, allocatable, dimension(:) :: id_eta_smooth
integer, allocatable, dimension(:) :: id_eta_smooth_on_nrho
integer, allocatable, dimension(:) :: id_pbot_smooth
integer, allocatable, dimension(:) :: id_prog
integer, allocatable, dimension(:) :: id_progT
integer, allocatable, dimension(:) :: id_prog_explicit
integer, allocatable, dimension(:) :: id_prog_rhodzt
integer, allocatable, dimension(:) :: id_prog_rhodztT
integer, allocatable, dimension(:) :: id_prog_int_rhodz
integer, allocatable, dimension(:) :: id_prog_int_rhodzT
integer, allocatable, dimension(:) :: id_prog_on_depth
integer, allocatable, dimension(:) :: id_tendency_conc
integer, allocatable, dimension(:) :: id_tendency_concL
integer, allocatable, dimension(:) :: id_tendency_concT
integer, allocatable, dimension(:) :: id_tendency
integer, allocatable, dimension(:) :: id_tendency_on_nrho
integer, allocatable, dimension(:) :: id_tendencyL
integer, allocatable, dimension(:) :: id_tendencyT
integer, allocatable, dimension(:) :: id_tendency_expl
integer, allocatable, dimension(:) :: id_tendency_explL
integer, allocatable, dimension(:) :: id_tendency_explT
integer, allocatable, dimension(:) :: id_surf_tracer
integer, allocatable, dimension(:) :: id_surf_tracerT
integer, allocatable, dimension(:) :: id_surf_tracer_sq
integer, allocatable, dimension(:) :: id_diag_surf_tracer
integer, allocatable, dimension(:) :: id_diag_surf_tracer_sq
integer, allocatable, dimension(:) :: id_bott_tracer
integer, allocatable, dimension(:) :: id_bott_tracerT
integer, allocatable, dimension(:) :: id_diag
integer, allocatable, dimension(:) :: id_diag_total
integer, allocatable, dimension(:) :: id_tmask_limit
integer :: id_neut_rho_tendency
integer :: id_pot_rho_tendency
integer :: id_neut_rho_tendency_on_nrho
integer :: id_wdian_rho_tendency
integer :: id_wdian_rho_tendency_on_nrho
integer :: id_tform_rho_tendency
integer :: id_tform_rho_tendency_on_nrho
integer :: id_neut_temp_tendency
integer :: id_neut_temp_tendency_on_nrho
integer :: id_wdian_temp_tendency
integer :: id_wdian_temp_tendency_on_nrho
integer :: id_tform_temp_tendency
integer :: id_tform_temp_tendency_on_nrho
integer :: id_neut_salt_tendency
integer :: id_neut_salt_tendency_on_nrho
integer :: id_wdian_salt_tendency
integer :: id_wdian_salt_tendency_on_nrho
integer :: id_tform_salt_tendency
integer :: id_tform_salt_tendency_on_nrho
integer :: id_neut_rho_smooth
integer :: id_pot_rho_smooth
integer :: id_neut_rho_smooth_on_nrho
integer :: id_wdian_rho_smooth
integer :: id_wdian_rho_smooth_on_nrho
integer :: id_tform_rho_smooth
integer :: id_tform_rho_smooth_on_nrho
integer :: id_eta_tend_smooth
integer :: id_eta_tend_smooth_glob
integer :: id_neut_temp_smooth
integer :: id_neut_temp_smooth_on_nrho
integer :: id_wdian_temp_smooth
integer :: id_wdian_temp_smooth_on_nrho
integer :: id_tform_temp_smooth
integer :: id_tform_temp_smooth_on_nrho
integer :: id_neut_salt_smooth
integer :: id_neut_salt_smooth_on_nrho
integer :: id_wdian_salt_smooth
integer :: id_wdian_salt_smooth_on_nrho
integer :: id_tform_salt_smooth
integer :: id_tform_salt_smooth_on_nrho
integer :: id_neut_rho_pme
integer :: id_pot_rho_pme
integer :: id_wdian_rho_pme
integer :: id_tform_rho_pme
integer :: id_neut_rho_pme_on_nrho
integer :: id_wdian_rho_pme_on_nrho
integer :: id_tform_rho_pme_on_nrho
integer :: id_neut_rho_pbl_pme_kn
integer :: id_neut_rho_pbl_pme_kn_on_nrho
integer :: id_wdian_rho_pbl_pme_kn
integer :: id_wdian_rho_pbl_pme_kn_on_nrho
integer :: id_tform_rho_pbl_pme_kn
integer :: id_tform_rho_pbl_pme_kn_on_nrho
integer :: id_neut_temp_pbl_pme_kn
integer :: id_neut_temp_pbl_pme_kn_on_nrho
integer :: idwdian_temp_pbl_pme_kn
integer :: idwdian_temp_pbl_pme_kn_on_nrho
integer :: idtform_temp_pbl_pme_kn
integer :: idtform_temp_pbl_pme_kn_on_nrho
integer :: id_neut_salt_pbl_pme_kn
integer :: id_neut_salt_pbl_pme_kn_on_nrho
integer :: idwdian_salt_pbl_pme_kn
integer :: idwdian_salt_pbl_pme_kn_on_nrho
integer :: idtform_salt_pbl_pme_kn
integer :: idtform_salt_pbl_pme_kn_on_nrho
integer :: id_neut_rho_pbl_pme_pr
integer :: id_neut_rho_pbl_pme_pr_on_nrho
integer :: id_wdian_rho_pbl_pme_pr
integer :: id_wdian_rho_pbl_pme_pr_on_nrho
integer :: id_tform_rho_pbl_pme_pr
integer :: id_tform_rho_pbl_pme_pr_on_nrho
integer :: id_neut_temp_pbl_pme_pr
integer :: id_neut_temp_pbl_pme_pr_on_nrho
integer :: idwdian_temp_pbl_pme_pr
integer :: idwdian_temp_pbl_pme_pr_on_nrho
integer :: idtform_temp_pbl_pme_pr
integer :: idtform_temp_pbl_pme_pr_on_nrho
integer :: id_neut_salt_pbl_pme_pr
integer :: id_neut_salt_pbl_pme_pr_on_nrho
integer :: idwdian_salt_pbl_pme_pr
integer :: idwdian_salt_pbl_pme_pr_on_nrho
integer :: idtform_salt_pbl_pme_pr
integer :: idtform_salt_pbl_pme_pr_on_nrho
integer :: id_frazil_on_nrho
integer :: id_neut_rho_frazil
integer :: id_pot_rho_frazil
integer :: id_wdian_rho_frazil
integer :: id_tform_rho_frazil
integer :: id_neut_rho_frazil_on_nrho
integer :: id_wdian_rho_frazil_on_nrho
integer :: id_tform_rho_frazil_on_nrho
integer :: id_eta_tend_frazil
integer :: id_eta_tend_frazil_glob
integer :: id_mass_t_on_nrho
integer :: id_mass_t_tendency_on_nrho
integer :: id_mass_nrho_layer
integer :: id_mass_nrho_tendency_layer
! for ascii output
integer :: unit=6
public update_ocean_tracer
public ocean_prog_tracer_init
public ocean_diag_tracer_init
public ocean_tracer_diagnostics_init
public ocean_tracer_end
public compute_tmask_limit
public ocean_tracer_restart
private update_advection_only
private remap_s_to_depth
private remap_depth_to_s
private inflow_nboundary_init
private watermass_diag
private send_tracer_diagnostics
!---------------nml settings---------------
logical :: zero_tendency = .false.
logical :: zero_tracer_source = .false.
logical :: debug_this_module = .false.
logical :: tmask_limit_ts_same = .true.
logical :: write_a_restart = .true.
logical :: remap_depth_to_s_init = .false.
logical :: inflow_nboundary = .false.
logical :: interpolate_tprog_to_pbott = .true.
logical :: interpolate_tdiag_to_pbott = .false.
logical :: limit_age_tracer = .false.
logical :: frazil_heating_before_vphysics = .false.
logical :: frazil_heating_after_vphysics = .false.
logical :: compute_tmask_limit_on = .true.
logical :: use_tempsalt_check_range = .false.
real :: age_tracer_max_init = 0.0
! for tracer package manager
logical :: ocean_tpm_debug = .false.
namelist /ocean_tracer_nml/ debug_this_module, zero_tendency, zero_tracer_source, write_a_restart, &
ocean_tpm_debug, tmask_limit_ts_same, remap_depth_to_s_init, &
inflow_nboundary, interpolate_tprog_to_pbott, interpolate_tdiag_to_pbott,&
limit_age_tracer, age_tracer_max_init, &
frazil_heating_before_vphysics, frazil_heating_after_vphysics, &
compute_tmask_limit_on, use_tempsalt_check_range
contains
!#######################################################################
! <FUNCTION NAME="ocean_prog_tracer_init">
!
! <DESCRIPTION>
! Initialization code for prognostic tracers, returning a pointer to
! the T_prog array.
! </DESCRIPTION>
!
function ocean_prog_tracer_init (Grid, Thickness, Ocean_options, Domain, Time, Time_steps, &
num_prog, vert_coordinate_type, obc, cmip_units, cmip_version, use_blobs, debug) &
result (T_prog)
type(ocean_grid_type), intent(in), target :: Grid
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_options_type), intent(inout) :: Ocean_options
type(ocean_domain_type), intent(in), target :: Domain
type(ocean_time_type), intent(in) :: Time
type(ocean_time_steps_type), intent(in) :: Time_steps
integer, intent(out) :: num_prog
integer, intent(in) :: vert_coordinate_type
logical, intent(in) :: obc
logical, intent(in) :: cmip_units
integer, intent(in) :: cmip_version
logical, intent(in) :: use_blobs
logical, intent(in), optional :: debug
! return value
type(ocean_prog_tracer_type), dimension(:), pointer :: T_prog
integer :: i, j, k, n, kb, l
integer :: ierr, num_diag
integer :: tau, taum1, taup1
integer :: ioun, io_status
integer, dimension(4) :: siz
integer :: frazil_heating_order=0
real :: fact
character(len=32) :: name
character(len=128) :: filename
logical :: initialize_as_a_passive_tracer=.false.
character(len=33) :: prog_name
character(len=138) :: prog_longname
character(len=48), parameter :: sub_name = 'ocean_prog_tracer_init'
character(len=256), parameter :: error_header = '==>Error from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
character(len=256), parameter :: warn_header = '==>Warning from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
character(len=256), parameter :: note_header = '==>Note from ' // trim(mod_name) // &
'(' // trim(sub_name) // '): '
! variables for tracer package
integer :: ind
real, dimension(2) :: range_array
character(len=64) :: caller_str
character(len=fm_string_len) :: string_fm
character(len=fm_type_name_len) :: typ
character(len=fm_string_len), pointer, dimension(:) :: good_list
integer :: stdoutunit, stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if (prog_module_initialized) then
call mpp_error(FATAL, trim(error_header) // ' Prognostic tracers already initialized')
endif
nullify(T_prog)
write( stdlogunit,'(/a/)') trim(version)
have_obc = obc
dtts = Time_steps%dtts
dtuv = Time_steps%dtuv
cp_oceanr = 1.0/cp_ocean
! provide for namelist over-ride
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_tracer_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_tracer_nml')
#else
ioun = open_namelist_file()
read (ioun, ocean_tracer_nml,iostat=io_status)
ierr = check_nml_error(io_status,'ocean_tracer_nml')
call close_file (ioun)
#endif
write (stdoutunit,'(/)')
write (stdoutunit, ocean_tracer_nml)
write (stdlogunit, ocean_tracer_nml)
if (PRESENT(debug) .and. .not. debug_this_module) then
debug_this_module = debug
endif
if(debug_this_module) then
write(stdoutunit,'(a)') '==>Note: running with debug_this_module=.true. so will print lots of checksums.'
endif
if(cmip_units .and. cmip_version <= 5) then
cmip_offset = kelvin
temp_units = 'K'
else
cmip_offset = 0.0
temp_units = 'degrees C'
endif
if(zero_tendency) then
call mpp_error(NOTE, trim(note_header) // ' zero_tendency=true so will not time step tracer fields.')
Ocean_options%tracer_tendency = 'Did NOT time step prognostic tracer fields.'
else
Ocean_options%tracer_tendency = 'Time stepped the prognostic tracer fields.'
endif
if(zero_tracer_source) then
call mpp_error(NOTE, &
trim(note_header) // ' zero_tracer_source=true so remove T_prog%source from evolution.')
endif
if(.not. write_a_restart) then
write(stdoutunit,'(a)') '==>Note: running ocean_tracer with write_a_restart=.false.'
write(stdoutunit,'(a)') ' Will NOT write restart file, and so cannot restart the run.'
endif
if(frazil_heating_before_vphysics) then
frazil_heating_order=frazil_heating_order+1
write(stdoutunit,'(a)') '==>Note: frazil heating called before vertical physics and before boundary fluxes.'
write(stdoutunit,'(a)') ' This method is retained for legacy purposes: it is NOT recommended for new runs. '
write(stdoutunit,'(a)') ' '
endif
if(frazil_heating_after_vphysics) then
frazil_heating_order=frazil_heating_order+1
write(stdoutunit,'(a)') '==>Note: frazil heating called after vertical physics and after boundary fluxes.'
write(stdoutunit,'(a)') ' This is the recommended method. '
write(stdoutunit,'(a)') ' '
endif
if(frazil_heating_order>1) then
write(stdoutunit,'(a)') '==>Error from ocean_tracer_mod: choose just one temporal order for frazil heating.'
write(stdoutunit,'(a)') ' '
call mpp_error(FATAL, &
trim(note_header) // ' can only choose one temporal order for frazil heating.')
endif
if(frazil_heating_order==0) then
write(stdoutunit,'(a)') '==>Error from ocean_tracer_mod: MUST specify order for frazil heating in ocean_tracer.'
write(stdoutunit,'(a)') ' '
call mpp_error(FATAL, &
trim(note_header) // ' Need to specify an order for calling frazil heating from ocean_tracer_mod.')
endif
! some time step information
if (dtts /= dtuv) then
write (stdoutunit,'(/a)') trim(warn_header) // ' Asynchronous timesteps (dtts > dtuv) imply inaccurate transients.'
write (stdoutunit,'(a)') ' and total tracer (i.e. heat content) is not conserved.'
write (stdoutunit,'(a,f5.2,a/)') ' dtts =',dtts/dtuv,' times larger than dtuv.'
else
call mpp_error(NOTE, trim(note_header) // ' Synchronous timesteps have been specified (dtts = dtuv).')
endif
tendency = Time_steps%tendency
dtime = Time_steps%dtime_t
dtimer = 1.0/(dtime+epsln)
tau = Time%tau
taum1 = Time%taum1
taup1 = Time%taup1
! initialize clock ids
id_clock_bih_tracer = mpp_clock_id('(Ocean tracer: bih tracer) ' ,grain=CLOCK_MODULE)
id_clock_lap_tracer = mpp_clock_id('(Ocean tracer: lap tracer) ' ,grain=CLOCK_MODULE)
id_clock_vert_diffuse = mpp_clock_id('(Ocean tracer: vert diffuse) ' ,grain=CLOCK_MODULE)
id_clock_vert_diffuse_implicit = mpp_clock_id('(Ocean tracer: vert diffuse impl) ',grain=CLOCK_MODULE)
id_clock_tracer_advect = mpp_clock_id('(Ocean tracer: advection) ' ,grain=CLOCK_MODULE)
id_clock_frazil = mpp_clock_id('(Ocean tracer: frazil) ' ,grain=CLOCK_MODULE)
id_clock_convection = mpp_clock_id('(Ocean tracer: convection) ' ,grain=CLOCK_MODULE)
id_clock_blob = mpp_clock_id('(Ocean tracer: Lagrangian blobs) ' ,grain=CLOCK_MODULE)
id_clock_blob_implicit = mpp_clock_id('(Ocean tracer: implicit blobs)' ,grain=CLOCK_ROUTINE)
id_clock_adjust_L_thickness = mpp_clock_id('(Ocean tracer: adjust the L thickness)',grain=CLOCK_ROUTINE)
id_clock_adjust_E_thickness = mpp_clock_id('(Ocean tracer: adjust the E thickness)',grain=CLOCK_ROUTINE)
! perform requested debugging for ocean tracer package manager
if (ocean_tpm_debug) then
write (stdoutunit,*) ' '
write (stdoutunit,*) 'Dumping field tree at start of ocean_prog_tracer_init'
call write_timestamp(Time%model_time)
if (.not. fm_dump_list('/', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping tracer tree')
endif
endif
! call routine to initialize the required tracer package
if (otpm_set_tracer_package('required', caller=trim(mod_name)//'('//trim(sub_name)//')') .le. 0) then
call mpp_error(FATAL, trim(error_header) // ' Could not set required packages list')
endif
! perform requested debugging for ocean tracer package manager
if (ocean_tpm_debug) then !{
write (stdoutunit,*) ' '
write (stdoutunit,*) 'Dumping /ocean_mod field tree before ocean_tpm_init'
call write_timestamp(Time%model_time)
if (.not. fm_dump_list('/ocean_mod', recursive = .true.)) then !{
call mpp_error(FATAL, trim(error_header) // ' Problem dumping tracer tree')
endif !}
endif !}
! call the initialization routine for the ocean tracer packages
call ocean_tpm_init(Domain, Grid, Time, Time_steps, &
Ocean_options, debug)
! perform requested debugging for ocean tracer package manager
if (ocean_tpm_debug) then
write (stdoutunit,*) ' '
write (stdoutunit,*) 'Dumping /ocean_mod field tree after ocean_tpm_init'
call write_timestamp(Time%model_time)
if (.not. fm_dump_list('/ocean_mod', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping tracer tree')
endif
endif
! check for any errors in the number of fields in the tracer_packages list
good_list => fm_util_get_string_array('/ocean_mod/GOOD/good_tracer_packages', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
if (associated(good_list)) then
call fm_util_check_for_bad_fields('/ocean_mod/tracer_packages', good_list, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
deallocate(good_list)
else
call mpp_error(FATAL,trim(error_header) // ' Empty "good_tracer_packages" list')
endif
! check for any errors in the number of fields in the prog_tracers list
good_list => fm_util_get_string_array('/ocean_mod/GOOD/good_prog_tracers', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
if (associated(good_list)) then
call fm_util_check_for_bad_fields('/ocean_mod/prog_tracers', good_list, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
deallocate(good_list)
else
call mpp_error(FATAL,trim(error_header) // ' Empty "good_prog_tracers" list')
endif
! check for any errors in the number of fields in the namelists list
good_list => fm_util_get_string_array('/ocean_mod/GOOD/good_namelists', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
if (associated(good_list)) then
call fm_util_check_for_bad_fields('/ocean_mod/namelists', good_list, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
deallocate(good_list)
!else
!call mpp_error(FATAL,trim(error_header) // ' Empty "good_namelists" list')
endif
! get the number of tracers
! for now, we will not try to dynamically allocate the tracer arrays here
num_prog_tracers = fm_get_length('/ocean_mod/prog_tracers')
num_diag = fm_get_length('/ocean_mod/diag_tracers')
num_prog=num_prog_tracers
! allocate arrays based on the number of prognostic tracers
allocate( T_prog (num_prog_tracers) )
allocate( id_eta_smooth (num_prog_tracers) )
allocate( id_eta_smooth_on_nrho(num_prog_tracers) )
allocate( id_pbot_smooth (num_prog_tracers) )
allocate( id_prog (num_prog_tracers) )
allocate( id_prog_explicit (num_prog_tracers) )
allocate( id_prog_rhodzt (num_prog_tracers) )
allocate( id_prog_int_rhodz (num_prog_tracers) )
allocate( id_prog_on_depth (num_prog_tracers) )
allocate( id_tendency_conc (num_prog_tracers) )
allocate( id_tendency (num_prog_tracers) )
allocate( id_tendency_on_nrho(num_prog_tracers) )
allocate( id_tendency_expl (num_prog_tracers) )
allocate( id_surf_tracer (num_prog_tracers) )
allocate( id_surf_tracer_sq (num_prog_tracers) )
allocate( id_bott_tracer (num_prog_tracers) )
if (use_blobs) then
allocate( id_progT (num_prog_tracers) )
allocate( id_prog_rhodztT (num_prog_tracers) )
allocate( id_prog_int_rhodzT(num_prog_tracers) )
allocate( id_surf_tracerT (num_prog_tracers) )
allocate( id_bott_tracerT (num_prog_tracers) )
allocate( id_tendency_concL (num_prog_tracers) )
allocate( id_tendency_concT (num_prog_tracers) )
allocate( id_tendencyL (num_prog_tracers) )
allocate( id_tendencyT (num_prog_tracers) )
allocate( id_tendency_explL (num_prog_tracers) )
allocate( id_tendency_explT (num_prog_tracers) )
endif
allocate( Tracer_restart (num_prog+num_diag) )
allocate( id_tracer_restart (num_prog+num_diag) )
allocate( tracer_restart_file (num_prog+num_diag) )
allocate( id_tracer_type (num_prog+num_diag) )
id_eta_smooth(:) = -1
id_eta_smooth_on_nrho(:)= -1
id_pbot_smooth(:) = -1
id_prog(:) = -1
id_prog_explicit(:) = -1
id_prog_rhodzt(:) = -1
id_prog_int_rhodz(:) = -1
id_prog_on_depth(:) = -1
id_tendency_conc(:) = -1
id_tendency(:) = -1
id_tendency_on_nrho(:)= -1
id_tendency_expl(:) = -1
id_surf_tracer(:) = -1
id_surf_tracer_sq(:) = -1
id_bott_tracer(:) = -1
if (use_blobs) then
id_progT(:) = -1
id_prog_rhodztT(:) = -1
id_prog_int_rhodzT(:)= -1
id_surf_tracerT(:) = -1
id_bott_tracerT(:) = -1
id_tendency_concL(:) = -1
id_tendency_concT(:) = -1
id_tendencyL(:) = -1
id_tendencyT(:) = -1
id_tendency_explL(:) = -1
id_tendency_explT(:) = -1
endif
! set logical to determine when to mpp_update tracers
do n=1,num_prog_tracers-1
T_prog(n)%complete=.false.
enddo
T_prog(num_prog_tracers)%complete=.true.
! dump the lists for the tracer packages
write (stdoutunit,*)
write (stdoutunit,*) 'Dumping tracer_packages tracer tree'
if (.not. fm_dump_list('/ocean_mod/tracer_packages', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping tracer_packages tracer tree')
endif
write (stdoutunit,*)
write (stdoutunit,*) 'Dumping prog_tracers tracer tree'
if (.not. fm_dump_list('/ocean_mod/prog_tracers', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping prog_tracers tracer tree')
endif
write (stdoutunit,*)
write (stdoutunit,*) 'Dumping namelists tracer tree'
if (.not. fm_dump_list('/ocean_mod/namelists', recursive = .true.)) then
call mpp_error(FATAL, trim(error_header) // ' Problem dumping namelists tracer tree')
endif
!### finished with initializing t_prog arrays
!### finished with call to tracer startup routine
! set local array indices
Grd => Grid
Dom => Domain
#ifndef MOM_STATIC_ARRAYS
call get_local_indices(Dom, isd, ied, jsd, jed, isc, iec, jsc, jec)
nk=Grd%nk
#endif
do n=1,num_prog_tracers
#ifndef MOM_STATIC_ARRAYS
allocate( T_prog(n)%field(isd:ied,jsd:jed,nk,3))
allocate( T_prog(n)%th_tendency(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%tendency(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%source(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%eta_smooth(isd:ied,jsd:jed))
allocate( T_prog(n)%pbot_smooth(isd:ied,jsd:jed))
allocate( T_prog(n)%wrk1(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%tmask_limit(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%K33_implicit(isd:ied,jsd:jed,nk))
if (use_blobs) then
allocate( T_prog(n)%fieldT(isd:ied,jsd:jed,nk))
allocate( T_prog(n)%sum_blob(isd:ied,jsd:jed,nk,3))
allocate( T_prog(n)%tend_blob(isd:ied,jsd:jed,nk))
endif
#endif
T_prog(n)%field(:,:,:,:) = 0.0
T_prog(n)%th_tendency(:,:,:) = 0.0
T_prog(n)%tendency(:,:,:) = 0.0
T_prog(n)%source(:,:,:) = 0.0
T_prog(n)%eta_smooth(:,:) = 0.0
T_prog(n)%pbot_smooth(:,:) = 0.0
T_prog(n)%wrk1(:,:,:) = 0.0
T_prog(n)%tmask_limit(:,:,:) = 0.0
T_prog(n)%K33_implicit(:,:,:) = 0.0
T_prog(n)%neutral_physics_limit = .false.
if (use_blobs) then
!note that fieldT is initialised in ocean_blob_init
T_prog(n)%fieldT(:,:,:) = 0.0
T_prog(n)%sum_blob(:,:,:,:) = 0.0
T_prog(n)%tend_blob(:,:,:) = 0.0
endif
enddo
if(inflow_nboundary) then
call inflow_nboundary_init
endif
! for global surface area normalization
cellarea_r = 1.0/(epsln + Grd%tcellsurf)
! fill the field table entries for the prognostic tracers
n = 0
do while (fm_loop_over_list('/ocean_mod/prog_tracers', name, typ, ind)) !{
if (typ .ne. 'list') then !{
call mpp_error(FATAL, trim(error_header) // ' ' // trim(name) // ' is not a list')
else !}{
n = n + 1 ! increment the array index
if (n .ne. ind) then !{
write (stdoutunit,*) trim(warn_header), ' Tracer index, ', ind, &
' does not match array index, ', n, ' for ', trim(name)
endif !}
! save the name
T_prog(n)%name = name
if (.not. fm_change_list('/ocean_mod/prog_tracers/' // trim(name))) then !{
call mpp_error(FATAL, trim(error_header) // ' Problem changing to ' // trim(name))
endif !}
caller_str = 'ocean_tracer_mod(ocean_prog_tracer_init)'
! save the units
T_prog(n)%units = fm_util_get_string('units', caller = caller_str, scalar = .true.)
! save the type
T_prog(n)%type = fm_util_get_string('type', caller = caller_str, scalar = .true.)
! save the longname
T_prog(n)%longname = fm_util_get_string('longname', caller = caller_str, scalar = .true.)
! save the conversion
T_prog(n)%conversion = fm_util_get_real('conversion', caller = caller_str, scalar = .true.)
! save the offset
T_prog(n)%offset = fm_util_get_real('offset', caller = caller_str, scalar = .true.)
! get the min and max of the tracer
T_prog(n)%min_tracer = fm_util_get_real('min_tracer', caller = caller_str, scalar = .true.)
T_prog(n)%max_tracer = fm_util_get_real('max_tracer', caller = caller_str, scalar = .true.)
! get the min and max of the range for analysis
T_prog(n)%min_range = fm_util_get_real('min_range', caller = caller_str, scalar = .true.)
T_prog(n)%max_range = fm_util_get_real('max_range', caller = caller_str, scalar = .true.)
! get the flux unit
T_prog(n)%flux_units = fm_util_get_string('flux_units', caller = caller_str, scalar = .true.)
! get the min and max of the flux range for analysis
T_prog(n)%min_flux_range = fm_util_get_real('min_flux_range', caller = caller_str, scalar = .true.)
T_prog(n)%max_flux_range = fm_util_get_real('max_flux_range', caller = caller_str, scalar = .true.)
! save the restart file
T_prog(n)%restart_file = fm_util_get_string('restart_file', caller = caller_str, scalar = .true.)
! save flag for whether the tracer must have a value in the restart file
T_prog(n)%const_init_tracer = fm_util_get_logical('const_init_tracer', caller = caller_str, scalar = .true.)
! save value to globally initialize this tracer (optional)
T_prog(n)%const_init_value = fm_util_get_real('const_init_value', caller = caller_str, scalar = .true.)
! get the horizontal-advection-scheme
string_fm = fm_util_get_string('horizontal-advection-scheme', caller = caller_str, scalar = .true.)
select case (trim(string_fm))
case ('upwind')
T_prog(n)%horz_advect_scheme = ADVECT_UPWIND
case ('2nd_order')
T_prog(n)%horz_advect_scheme = ADVECT_2ND_ORDER
case ('4th_order')
T_prog(n)%horz_advect_scheme = ADVECT_4TH_ORDER
case ('quicker')
T_prog(n)%horz_advect_scheme = ADVECT_QUICKER
case ('quickMOM3')
T_prog(n)%horz_advect_scheme = ADVECT_QUICKMOM3
case ('6th_order')
T_prog(n)%horz_advect_scheme = ADVECT_6TH_ORDER
case ('mdfl_sup_b')
T_prog(n)%horz_advect_scheme = ADVECT_MDFL_SUP_B
case ('mdfl_sweby')
T_prog(n)%horz_advect_scheme = ADVECT_MDFL_SWEBY
case ('mdfl_sweby_test')
T_prog(n)%horz_advect_scheme = ADVECT_MDFL_SWEBY_TEST
case ('dst_linear')