-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathMOM_generic_tracer.F90
1160 lines (986 loc) · 59 KB
/
MOM_generic_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
!> Drives the generic version of tracers TOPAZ and CFC and other GFDL BGC components
module MOM_generic_tracer
! This file is part of MOM6. See LICENSE.md for the license.
#include <MOM_memory.h>
! The following macro is usually defined in <fms_platform.h> but since MOM6 should not directly
! include files from FMS we replicate the macro lines here:
#ifdef NO_F2000
#define _ALLOCATED associated
#else
#define _ALLOCATED allocated
#endif
! ### These imports should not reach into FMS directly ###
use field_manager_mod, only: fm_string_len
use generic_tracer, only: generic_tracer_register, generic_tracer_get_diag_list
use generic_tracer, only: generic_tracer_init, generic_tracer_source, generic_tracer_register_diag
use generic_tracer, only: generic_tracer_coupler_get, generic_tracer_coupler_set
use generic_tracer, only: generic_tracer_end, generic_tracer_get_list, do_generic_tracer
use generic_tracer, only: generic_tracer_update_from_bottom,generic_tracer_vertdiff_G
use generic_tracer, only: generic_tracer_coupler_accumulate
use g_tracer_utils, only: g_tracer_get_name,g_tracer_set_values,g_tracer_set_common,g_tracer_get_common
use g_tracer_utils, only: g_tracer_get_next,g_tracer_type,g_tracer_is_prog,g_tracer_flux_init
use g_tracer_utils, only: g_tracer_send_diag,g_tracer_get_values
use g_tracer_utils, only: g_tracer_get_pointer,g_tracer_get_alias,g_tracer_set_csdiag
use g_tracer_utils, only: g_tracer_get_obc_segment_props
use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS
use MOM_ALE_sponge, only : ALE_sponge_CS, initialize_ALE_sponge
use MOM_coms, only : EFP_type, max_across_PEs, min_across_PEs, PE_here
use MOM_diagnose_mld, only : diagnoseMLDbyDensityDifference, diagnoseMLDbyEnergy
use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use MOM_diag_mediator, only : diag_ctrl, get_diag_time_end
use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_forcing_type, only : forcing, optics_type
use MOM_grid, only : ocean_grid_type
use MOM_hor_index, only : hor_index_type
use MOM_interface_heights, only : thickness_to_dz
use MOM_io, only : file_exists, MOM_read_data, slasher
use MOM_open_boundary, only : ocean_OBC_type
use MOM_open_boundary, only : register_obgc_segments, fill_obgc_segments
use MOM_open_boundary, only : set_obgc_segments_props
use MOM_restart, only : register_restart_field, query_initialized, set_initialized, MOM_restart_CS
use MOM_spatial_means, only : global_area_mean, global_mass_int_EFP, array_global_min_max
use MOM_sponge, only : set_up_sponge_field, sponge_CS
use MOM_time_manager, only : time_type, set_time
use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_Z_init, only : tracer_Z_init
use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : surface, thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
implicit none ; private
!> A state hidden in module data that is very much not allowed in MOM6
! ### This needs to be fixed
logical :: g_registered = .false.
public register_MOM_generic_tracer, initialize_MOM_generic_tracer
public MOM_generic_tracer_column_physics, MOM_generic_tracer_surface_state
public end_MOM_generic_tracer, MOM_generic_tracer_get
public MOM_generic_tracer_stock
public MOM_generic_flux_init
public MOM_generic_tracer_min_max
public MOM_generic_tracer_fluxes_accumulate
public register_MOM_generic_tracer_segments
!> Control structure for generic tracers
type, public :: MOM_generic_tracer_CS ; private
character(len = 200) :: IC_file !< The file in which the generic tracer initial values can
!! be found, or an empty string for internal initialization.
logical :: Z_IC_file !< If true, the generic_tracer IC_file is in Z-space. The default is false.
real :: tracer_IC_val = 0.0 !< The initial value assigned to tracers, in
!! concentration units [conc]
real :: tracer_land_val = -1.0 !< The values of tracers used where land is masked out, in
!! concentration units [conc]
logical :: tracers_may_reinit !< If true, tracers may go through the
!! initialization code if they are not found in the restart files.
logical :: mld_pha_calc = .False. !< If true, use a fixed value for photoacclimation MLD
real :: mld_pha_val = 0.0 !< The value of fixed photoacclimation MLD
logical :: mld_pha_use_delta_rho = .False. !< If true, use a density diference to find the MLD
real :: mld_pha_href = 0.0 !< The reference depth for density difference based MLD
real :: mld_pha_drho = 0.03 !< The density thershold for a density difference based MLD
logical :: mld_pha_use_delta_eng = .False. !< If true, use an energy diference to find the MLD
real :: mld_pha_deng = 25.0 !< The energy threshold for an energy d ifference based MLD
type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
!! regulate the timing of diagnostic output.
type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Restart control structure
type(ocean_OBC_type), pointer :: OBC => NULL() !<open boundary condition type
!> Pointer to the first element of the linked list of generic tracers.
type(g_tracer_type), pointer :: g_tracer_list => NULL()
end type MOM_generic_tracer_CS
contains
!> Initializes the generic tracer packages and adds their tracers to the list
!! Adds the tracers in the list of generic tracers to the set of MOM tracers (i.e., MOM-register them)
!! Register these tracers for restart
function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
type(hor_index_type), intent(in) :: HI !< Horizontal index ranges
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module
type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer
!! advection and diffusion module.
type(MOM_restart_CS), target, intent(inout) :: restart_CS !< MOM restart control struct
! Local variables
logical :: register_MOM_generic_tracer
logical :: obc_has
! This include declares and sets the variable "version".
# include "version_variable.h"
character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer'
character(len=200) :: inputdir ! The directory where NetCDF input files are.
! These can be overridden later in via the field manager?
integer :: ntau, axes(3)
type(g_tracer_type), pointer :: g_tracer, g_tracer_next
character(len=fm_string_len) :: g_tracer_name, longname,units
character(len=fm_string_len) :: obc_src_file_name, obc_src_field_name
real :: lfac_in ! Multiplicative factor used in setting the tracer-specific inverse length
! scales associated with inflowing tracer reservoirs at OBCs [nondim]
real :: lfac_out ! Multiplicative factor used in setting the tracer-specific inverse length
! scales associated with outflowing tracer reservoirs at OBCs [nondim]
real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc]
real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc]
real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim]
integer, dimension(SZI_(HI),SZJ_(HI)) :: grid_kmt ! A 2-d array of nk
register_MOM_generic_tracer = .false.
if (associated(CS)) then
call MOM_error(FATAL, "register_MOM_generic_tracer called with an "// &
"associated control structure.")
endif
allocate(CS)
!Register all the generic tracers used and create the list of them.
!This can be called by ALL PE's. No array fields allocated.
if (.not. g_registered) then
call generic_tracer_register()
g_registered = .true.
endif
! Read all relevant parameters and write them to the model log.
call log_version(param_file, sub_name, version, "")
call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE", CS%IC_file, &
"The file in which the generic tracer initial values can "//&
"be found, or an empty string for internal initialization.", &
default=" ")
if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then
! Add the directory if CS%IC_file is not already a complete path.
call get_param(param_file, sub_name, "INPUTDIR", inputdir, default=".")
CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file)
call log_param(param_file, sub_name, "INPUTDIR/GENERIC_TRACER_IC_FILE", CS%IC_file)
endif
call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE_IS_Z", CS%Z_IC_file, &
"If true, GENERIC_TRACER_IC_FILE is in depth space, not "//&
"layer space.",default=.false.)
call get_param(param_file, sub_name, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, &
"If true, tracers may go through the initialization code "//&
"if they are not found in the restart files. Otherwise "//&
"it is a fatal error if tracers are not found in the "//&
"restart files of a restarted run.", default=.false.)
CS%restart_CSp => restart_CS
ntau=1 ! MOM needs the fields at only one time step
! At this point G%mask2dT and CS%diag%axesTL are not allocated.
! postpone diag_registeration to initialize_MOM_generic_tracer
!Fields cannot be diag registered as they are allocated and have to registered later.
grid_tmask(:,:,:) = 0.0
grid_kmt(:,:) = 0
axes(:) = -1
!
! Initialize all generic tracers
!
call generic_tracer_init(HI%isc,HI%iec,HI%jsc,HI%jec,HI%isd,HI%ied,HI%jsd,HI%jed,&
GV%ke,ntau,axes,grid_tmask,grid_kmt,set_time(0,0))
!
! MOM-register the generic tracers
!
!Get the tracer list
call generic_tracer_get_list(CS%g_tracer_list)
if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//&
": No tracer in the list.")
! For each tracer name get its T_prog index and get its fields
g_tracer=>CS%g_tracer_list
do
call g_tracer_get_alias(g_tracer,g_tracer_name)
call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field)
call g_tracer_get_values(g_tracer,g_tracer_name,'longname', longname)
call g_tracer_get_values(g_tracer,g_tracer_name,'units',units )
!!nnz: MOM field is 3D. Does this affect performance? Need it be override field?
tr_ptr => tr_field(:,:,:,1)
! Register prognostic tracer for horizontal advection, diffusion, and restarts.
if (g_tracer_is_prog(g_tracer)) then
call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, &
name=g_tracer_name, longname=longname, units=units, &
registry_diags=.true., & !### CHANGE TO TRUE?
restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit)
else
call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, &
restart_CS, longname=longname, units=units)
endif
!traverse the linked list till hit NULL
call g_tracer_get_next(g_tracer, g_tracer_next)
if (.NOT. associated(g_tracer_next)) exit
g_tracer=>g_tracer_next
enddo
register_MOM_generic_tracer = .true.
end function register_MOM_generic_tracer
!> Register OBC segments for generic tracers
subroutine register_MOM_generic_tracer_segments(CS, GV, OBC, tr_Reg, param_file)
type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether,
!! where, and what open boundary conditions are used.
type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer
!! advection and diffusion module.
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
! Local variables
logical :: obc_has
! This include declares and sets the variable "version".
# include "version_variable.h"
character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer_segments'
type(g_tracer_type), pointer :: g_tracer,g_tracer_next
character(len=fm_string_len) :: g_tracer_name
character(len=fm_string_len) :: obc_src_file_name, obc_src_field_name
real :: lfac_in ! Multiplicative factor used in setting the tracer-specific inverse length
! scales associated with inflowing tracer reservoirs at OBCs [nondim]
real :: lfac_out ! Multiplicative factor used in setting the tracer-specific inverse length
! scales associated with outflowing tracer reservoirs at OBCs [nondim]
if (.NOT. associated(OBC)) return
!Get the tracer list
call generic_tracer_get_list(CS%g_tracer_list)
if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//&
": No tracer in the list.")
g_tracer=>CS%g_tracer_list
do
call g_tracer_get_alias(g_tracer,g_tracer_name)
if (g_tracer_is_prog(g_tracer)) then
call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ,&
obc_src_file_name,obc_src_field_name,lfac_in,lfac_out)
if (obc_has) then
call set_obgc_segments_props(OBC,g_tracer_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out)
call register_obgc_segments(GV, OBC, tr_Reg, param_file, g_tracer_name)
endif
endif
!traverse the linked list till hit NULL
call g_tracer_get_next(g_tracer, g_tracer_next)
if (.NOT. associated(g_tracer_next)) exit
g_tracer=>g_tracer_next
enddo
end subroutine register_MOM_generic_tracer_segments
!> Initialize phase II: Initialize required variables for generic tracers
!! There are some steps of initialization that cannot be done in register_MOM_generic_tracer
!! This is the place and time to do them:
!! Set the grid mask and initial time for all generic tracers.
!! Diag_register them.
!! Z_diag_register them.
!!
!! This subroutine initializes the NTR tracer fields in tr(:,:,:,:)
!! and it sets up the tracer output.
subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, tv, param_file, diag, OBC, &
CS, sponge_CSp, ALE_sponge_CSp)
logical, intent(in) :: restart !< .true. if the fields have already been
!! read from a restart file.
type(time_type), target, intent(in) :: day !< Time of the start of the run.
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic
!! variables
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
type(diag_ctrl), target, intent(in) :: diag !< Regulates diagnostic output.
type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether,
!! where, and what open boundary conditions are used.
type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module.
type(sponge_CS), pointer :: sponge_CSp !< Pointer to the control structure for the sponges.
type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< Pointer to the control structure for the
!! ALE sponges.
character(len=128), parameter :: sub_name = 'initialize_MOM_generic_tracer'
logical :: OK,obc_has
logical :: do_use_GT_sponge
integer :: i, j, k, isc, iec, jsc, jec, nk
type(g_tracer_type), pointer :: g_tracer,g_tracer_next
character(len=fm_string_len) :: g_tracer_name
real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc]
real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Layer vertical extent [Z ~> m]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim]
integer, dimension(SZI_(G),SZJ_(G)) :: grid_kmt ! A 2-d array of nk
!! 2010/02/04 Add code to re-initialize Generic Tracers if needed during a model simulation
!! By default, restart cpio should not contain a Generic Tracer IC file and step below will be skipped.
!! Ideally, the generic tracer IC file should have the tracers on Z levels.
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke
CS%diag=>diag
!Get the tracer list
if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//&
": No tracer in the list.")
!For each tracer name get its fields
g_tracer=>CS%g_tracer_list
call thickness_to_dz(h, tv, dz, G, GV, US)
do
if (INDEX(CS%IC_file, '_NULL_') /= 0) then
call MOM_error(WARNING, "The name of the IC_file "//trim(CS%IC_file)//&
" indicates no MOM initialization was asked for the generic tracers."//&
"Bypassing the MOM initialization of ALL generic tracers!")
exit
endif
call g_tracer_get_alias(g_tracer,g_tracer_name)
call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field)
tr_ptr => tr_field(:,:,:,1)
if (.not.restart .or. (CS%tracers_may_reinit .and. &
.not.query_initialized(tr_ptr, g_tracer_name, CS%restart_CSp))) then
if (g_tracer%requires_src_info ) then
call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
"initializing generic tracer "//trim(g_tracer_name)//&
" using MOM_initialize_tracer_from_Z ")
call MOM_initialize_tracer_from_Z(dz, tr_ptr, G, GV, US, param_file, &
src_file=g_tracer%src_file, src_var_nam=g_tracer%src_var_name, &
src_var_unit_conversion=g_tracer%src_var_unit_conversion, &
src_var_record=g_tracer%src_var_record, src_var_gridspec=g_tracer%src_var_gridspec, &
h_in_Z_units=.true.)
!Check/apply the bounds for each g_tracer
do k=1,nk ; do j=jsc,jec ; do i=isc,iec
if (tr_ptr(i,j,k) /= CS%tracer_land_val) then
if (tr_ptr(i,j,k) < g_tracer%src_var_valid_min) tr_ptr(i,j,k) = g_tracer%src_var_valid_min
!Jasmin does not want to apply the maximum for now
!if (tr_ptr(i,j,k) > g_tracer%src_var_valid_max) tr_ptr(i,j,k) = g_tracer%src_var_valid_max
endif
enddo ; enddo ; enddo
!jgj: Reset CASED to 0 below K=1
if ( (trim(g_tracer_name) == 'cased') .or. (trim(g_tracer_name) == 'ca13csed') ) then
do k=2,nk ; do j=jsc,jec ; do i=isc,iec
if (tr_ptr(i,j,k) /= CS%tracer_land_val) then
tr_ptr(i,j,k) = 0.0
endif
enddo ; enddo ; enddo
endif
elseif(.not. g_tracer%requires_restart) then
!Do nothing for this tracer, it is initialized by the tracer package
call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
"skip initialization of generic tracer "//trim(g_tracer_name))
else !Do it old way if the tracer is not registered to start from a specific source file.
!This path should be deprecated if all generic tracers are required to start from specified sources.
if (len_trim(CS%IC_file) > 0) then
! Read the tracer concentrations from a netcdf file.
if (.not.file_exists(CS%IC_file)) call MOM_error(FATAL, &
"initialize_MOM_Generic_tracer: Unable to open "//CS%IC_file)
if (CS%Z_IC_file) then
OK = tracer_Z_init(tr_ptr, h, CS%IC_file, g_tracer_name, G, GV, US)
if (.not.OK) then
OK = tracer_Z_init(tr_ptr, h, CS%IC_file, trim(g_tracer_name), G, GV, US)
if (.not.OK) call MOM_error(FATAL,"initialize_MOM_Generic_tracer: "//&
"Unable to read "//trim(g_tracer_name)//" from "//&
trim(CS%IC_file)//".")
endif
call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
"initialized generic tracer "//trim(g_tracer_name)//&
" using Generic Tracer File on Z: "//CS%IC_file)
else
! native grid
call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
"Using Generic Tracer IC file on native grid "//trim(CS%IC_file)//&
" for tracer "//trim(g_tracer_name))
call MOM_read_data(CS%IC_file, trim(g_tracer_name), tr_ptr, G%Domain)
endif
else
call MOM_error(FATAL,"initialize_MOM_generic_tracer: "//&
"check Generic Tracer IC filename "//trim(CS%IC_file)//&
" for tracer "//trim(g_tracer_name))
endif
endif
call set_initialized(tr_ptr, g_tracer_name, CS%restart_CSp)
endif
call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has )
if(obc_has .and. g_tracer_is_prog(g_tracer) .and. .not.restart) &
call fill_obgc_segments(G, GV, OBC, tr_ptr, g_tracer_name)
!traverse the linked list till hit NULL
call g_tracer_get_next(g_tracer, g_tracer_next)
if (.NOT. associated(g_tracer_next)) exit
g_tracer=>g_tracer_next
enddo
!! end section to re-initialize generic tracers
call get_param(param_file, "MOM", "PHA_MLD_CALC", CS%mld_pha_calc, &
"If false, use a fixed value for the photoacclimation mixed layer depth within the "//&
"generic tracer update. This MLD is only used for photoacclimation. This variable should "//&
"be set to true if using COBALTv3 for the BGC.", default=.false.)
if (CS%mld_pha_calc) then
call get_param(param_file, "MOM", "PHA_MLD_USE_DELTA_RHO", CS%mld_pha_use_delta_rho, &
"If true, use a density difference to find the photoacclimation mixed layer depth "//&
"within the generic tracer update. This MLD is only used for photoacclimation.", default=.false.)
call get_param(param_file, "MOM", "PHA_MLD_USE_DELTA_ENG", CS%mld_pha_use_delta_eng, &
"If true, use an energy difference to find the photoacclimation mixed layer depth "//&
"with the generic tracer update. This MLD is only used for photoacclimation.", default=.false.)
if (CS%mld_pha_use_delta_rho .and. CS%mld_pha_use_delta_eng) then
call MOM_error(FATAL, "PHA_MLD_CALC is set to true and PHA_MLD_USE_DELTA_RHO and PHA_MLD_USE_DELTA_ENG "// &
"are both true. Choose only one option for the calculated photoacclimation MLD!")
elseif ((.not.CS%mld_pha_use_delta_rho) .and. (.not.CS%mld_pha_use_delta_eng)) then
call MOM_error(FATAL, "PHA_MLD_CALC is set to true but PHA_MLD_USE_DELTA_RHO and PHA_MLD_USE_DELTA_ENG "// &
"are both false. Choose an option for the calculated photoacclimation MLD!")
endif
if (CS%mld_pha_use_delta_rho) then
call get_param(param_file, "MOM", "PHA_MLD_HREF", CS%mld_pha_href, &
"The reference depth for a density difference based photoacclimation MLD [m].", &
units='m', default=0.0, scale=US%m_to_Z, do_not_log=.not.CS%mld_pha_use_delta_rho)
call get_param(param_file, "MOM", "PHA_MLD_DRHO", CS%mld_pha_drho, &
"The density difference for a density difference based photoacclimation MLD [kg m-3].", &
units='kg/m3', default=0.03, scale=US%kg_m3_to_R, do_not_log=.not.CS%mld_pha_use_delta_rho)
elseif (CS%mld_pha_use_delta_eng) then
call get_param(param_file, "MOM", "PHA_MLD_DENG", CS%mld_pha_deng, &
"The energy for an energy difference based photoacclimation MLD.", default=25.0, &
units='J/m2',scale=US%W_m2_to_RZ3_T3*US%s_to_T, do_not_log=.not.CS%mld_pha_use_delta_eng)
endif
else
call get_param(param_file, "MOM", "PHA_MLD_VAL", CS%mld_pha_val, &
"The depth of photoacclimation if fixed depth is used [m].", &
units='m', default=0.0, scale=US%m_to_Z)
endif
!Now we can reset the grid mask, axes and time to their true values
!Note that grid_tmask must be set correctly on the data domain boundary
!so that coast mask can be deduced from it.
grid_tmask(:,:,:) = 0.0
grid_kmt(:,:) = 0
do j = G%jsd, G%jed ; do i = G%isd, G%ied
if (G%mask2dT(i,j) > 0.0) then
grid_tmask(i,j,:) = 1.0
grid_kmt(i,j) = GV%ke ! Tell the code that a layer thicker than 1m is the bottom layer.
endif
enddo ; enddo
call g_tracer_set_common(G%isc,G%iec,G%jsc,G%jec,G%isd,G%ied,G%jsd,G%jed,&
GV%ke,1,CS%diag%axesTL%handles,grid_tmask,grid_kmt,day)
call get_param(param_file, "initialize_sponges_file", "DO_SPONGE_GENERIC_TRACER", do_use_gt_sponge, &
"If true, then some generic tracers may be nudged.", default=.false.)
if (do_use_GT_sponge) call g_tracer_initialize_sponges(G, GV, US, CS, param_file, sponge_CSp, ALE_sponge_CSp, day)
! Register generic tracer modules diagnostics
#ifdef _USE_MOM6_DIAG
call g_tracer_set_csdiag(CS%diag)
#endif
call generic_tracer_register_diag()
#ifdef _USE_MOM6_DIAG
call g_tracer_set_csdiag(CS%diag)
#endif
end subroutine initialize_MOM_generic_tracer
subroutine g_tracer_initialize_sponges(G, GV, US, CS, param_file, Layer_CSp, ALE_CSp, Time)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G)) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module.
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
type(sponge_CS), pointer :: Layer_CSp !< A pointer that is set to point to the control
!! structure for this module (in layered mode).
type(ALE_sponge_CS), pointer :: ALE_CSp !< A pointer that is set to point to the control
!! structure for this module (in ALE mode).
type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in
!! overrides any value set for Time.
! Local variables
real, allocatable, dimension(:,:,:) :: eta ! The target interface heights [Z ~> m].
real, allocatable, dimension(:,:,:) :: dz ! The target interface thicknesses in height units [Z ~> m]
real, allocatable, dimension(:,:,:) :: h ! The target interface thicknesses [H ~> m or kg m-2].
real, allocatable, dimension(:,:,:) :: tmp_GT ! A temporary array for reading sponge target temperatures
! on the vertical grid of the input file [C ~> degC]
real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1]
integer :: i, j, k, is, ie, js, je, nz
integer :: isd, ied, jsd, jed
integer, dimension(4) :: siz
integer :: nz_data ! The size of the sponge source grid
character(len=40) :: tmp_var, Idamp_var, eta_var
character(len=40) :: mdl = "initialize_sponges_file"
character(len=200) :: damping_file, state_file, state_uv_file ! Strings for filenames
character(len=200) :: filename, inputdir ! Strings for file/path and path.
character(len=200) :: g_tracer_name
type(g_tracer_type), pointer :: g_tracer,g_tracer_next
logical :: do_sponge_gt
real, dimension(:,:,:), pointer :: g_tracer_ptr
logical :: use_ALE ! True if ALE is being used, False if in layered mode
logical :: time_space_interp_sponge ! If true use sponge data that need to be interpolated in both
! the horizontal dimension and in time prior to vertical remapping.
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
Idamp(:,:) = 0.0
call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
call get_param(param_file, mdl, "SPONGE_GT_DAMPING_FILE", damping_file, &
"The name of the file with the sponge damping rates.") !, &
call get_param(param_file, mdl, "SPONGE_IDAMP_VAR", Idamp_var, &
"The name of the inverse damping rate variable in "//&
"SPONGE_DAMPING_FILE.", default="Idamp")
call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.)
call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, &
"If True, perform on-the-fly regridding in lat-lon-time of sponge restoring data.", &
default=.false.)
! Read in sponge damping rate for tracers
filename = trim(inputdir)//trim(damping_file)
call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename)
if (.not.file_exists(filename, G%Domain)) &
call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
call MOM_read_data(filename, Idamp_var, Idamp(:,:), G%Domain, scale=US%T_to_s)
! Now register all of the fields which are damped in the sponge.
! By default, momentum is advected vertically within the sponge, but
! momentum is typically not damped within the sponge.
g_tracer=>CS%g_tracer_list
do ! check for each generic tracer if it is nudged
call g_tracer_get_alias(g_tracer,g_tracer_name)
call get_param(param_file, mdl, "DO_SPONGE_GT_"//trim(g_tracer_name), do_sponge_gt, &
"If true, then generic tracer "//trim(g_tracer_name)//" is nudged.", &
default=.false.)
if (do_sponge_gt) then
call get_param(param_file, mdl, "SPONGE_GT_"//trim(g_tracer_name)//"_FILE", state_file, &
"The name of the file with the state to damp the generic ", &
"tracer "//trim(g_tracer_name)//" toward.", default="sponge_"//trim(g_tracer_name)//".nc.")
call get_param(param_file, mdl, "SPONGE_GT_"//trim(g_tracer_name)//"_VAR", tmp_var, &
"The name of the variable to use in the sponge_GT file for generic ", &
"tracer "//trim(g_tracer_name)//".", default=trim(g_tracer_name))
call get_param(param_file, mdl, "SPONGE_GT_"//trim(g_tracer_name)//"_ETA_VAR", eta_var, &
"The name of the interface height variable in "//&
"SPONGE_GT_"//trim(g_tracer_name)//"_FILE.", default="ETA")
filename = trim(inputdir)//trim(state_file)
call log_param(param_file, mdl, "INPUTDIR/SPONGE_GT"//trim(g_tracer_name)//"_FILE", filename)
if (.not.file_exists(filename, G%Domain)) &
call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
! get the pointer for this tracer
call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',g_tracer_ptr)
if (use_ALE) then ! ALE mode
if (.not. time_space_interp_sponge) then
!call field_size(filename,eta_var,siz,no_domain=.true.)
if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) &
call MOM_error(FATAL,"initialize_sponge_file: Array size mismatch for sponge data.")
nz_data = siz(3)-1
allocate(eta(isd:ied,jsd:jed,nz_data+1))
allocate(dz(isd:ied,jsd:jed,nz_data))
call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z)
do j=js,je ; do i=is,ie
eta(i,j,nz_data+1) = -depth_tot(i,j)
enddo ; enddo
do k=nz_data,1,-1 ; do j=js,je ; do i=is,ie
if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) &
eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z
enddo ; enddo ; enddo
do k=1,nz_data ; do j=js,je ; do i=is,ie
dz(i,j,k) = eta(i,j,k)-eta(i,j,k+1)
enddo; enddo ; enddo
deallocate(eta)
allocate(tmp_GT(isd:ied,jsd:jed,nz_data))
call MOM_read_data(filename, tmp_var, tmp_GT(:,:,:), G%Domain) !, scale=US%degC_to_C)
call set_up_ALE_sponge_field(tmp_GT, G, GV, g_tracer_ptr, &
ALE_CSp, trim(g_tracer_name))
deallocate(tmp_GT)
deallocate(dz)
else
call set_up_ALE_sponge_field(filename, tmp_var, Time, G, GV, US, g_tracer_ptr, &
ALE_CSp, trim(g_tracer_name))
endif
endif
endif
call g_tracer_get_next(g_tracer, g_tracer_next)
if (.NOT. associated(g_tracer_next)) exit
g_tracer=>g_tracer_next
enddo
end subroutine g_tracer_initialize_sponges
!> Column physics for generic tracers.
!! Get the coupler values for generic tracers that exchange with atmosphere
!! Update generic tracer concentration fields from sources and sinks.
!! Vertically diffuse generic tracer concentration fields.
!! Update generic tracers from bottom and their bottom reservoir.
!!
!! This subroutine applies diapycnal diffusion and any other column
!! tracer physics or chemistry to the tracers from this file.
!! CFCs are relatively simple, as they are passive tracers. with only a surface
!! flux as a source.
subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, US, CS, tv, optics, &
evap_CFL_limit, minimum_forcing_depth)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: ea !< The amount of fluid entrained from the layer
!! above during this call [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: eb !< The amount of fluid entrained from the layer
!! below during this call [H ~> m or kg m-2].
type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
!! and tracer forcing fields.
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Hml !< Mixed layer depth [Z ~> m]
real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module.
type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
type(optics_type), intent(in) :: optics !< The structure containing optical properties.
real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can
!! be fluxed out of the top layer in a timestep [nondim]
! Stored previously in diabatic CS.
real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which fluxes
!! can be applied [H ~> m or kg m-2]
! Stored previously in diabatic CS.
! The arguments to this subroutine are redundant in that
! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
! Local variables
character(len=128), parameter :: sub_name = 'MOM_generic_tracer_column_physics'
type(g_tracer_type), pointer :: g_tracer, g_tracer_next
character(len=fm_string_len) :: g_tracer_name
real, dimension(:,:), pointer :: stf_array ! The surface flux of the tracer [conc kg m-2 s-1]
real, dimension(:,:), pointer :: trunoff_array ! The tracer concentration in the river runoff [conc]
real, dimension(:,:), pointer :: runoff_tracer_flux_array ! The runoff tracer flux [conc kg m-2 s-1]
real :: surface_field(SZI_(G),SZJ_(G)) ! The surface value of some field, here only used for salinity [S ~> ppt]
real :: dz_ml(SZI_(G),SZJ_(G)) ! The mixed layer depth in the MKS units used for generic tracers [m]
real :: sosga ! The global mean surface salinity [ppt]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: rho_dzt ! Layer mass per unit area [kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt ! Layer vertical extents [m]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! A work array of thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G)) :: mld_pha ! The mixed layer depth calculated for photoacclimation
! that is used in COBALTv3
integer :: i, j, k, isc, iec, jsc, jec, nk
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke
!Get the tracer list
if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL,&
trim(sub_name)//": No tracer in the list.")
#ifdef _USE_MOM6_DIAG
call g_tracer_set_csdiag(CS%diag)
#endif
!
!Extract the tracer surface fields from coupler and update tracer fields from sources
!
!call generic_tracer_coupler_get(fluxes%tr_fluxes)
!Niki: This is moved out to ocean_model_MOM.F90 because if dt_therm>dt_cpld we need to average
! the fluxes without coming into this subroutine.
! MOM5 has to modified to conform.
!
!Add contribution of river to surface flux
!
g_tracer=>CS%g_tracer_list
do
if (_ALLOCATED(g_tracer%trunoff) .and. (.NOT. g_tracer%runoff_added_to_stf)) then
call g_tracer_get_alias(g_tracer,g_tracer_name)
call g_tracer_get_pointer(g_tracer,g_tracer_name,'stf', stf_array)
call g_tracer_get_pointer(g_tracer,g_tracer_name,'trunoff',trunoff_array)
call g_tracer_get_pointer(g_tracer,g_tracer_name,'runoff_tracer_flux',runoff_tracer_flux_array)
!nnz: Why is fluxes%river = 0?
runoff_tracer_flux_array(:,:) = trunoff_array(:,:) * &
US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:)
stf_array = stf_array + runoff_tracer_flux_array
g_tracer%runoff_added_to_stf = .true.
endif
!traverse the linked list till hit NULL
call g_tracer_get_next(g_tracer, g_tracer_next)
if (.NOT. associated(g_tracer_next)) exit
g_tracer => g_tracer_next
enddo
!
!Prepare input arrays for source update
!
rho_dzt(:,:,:) = GV%H_to_kg_m2 * GV%Angstrom_H
do k=1,nk ; do j=jsc,jec ; do i=isc,iec
rho_dzt(i,j,k) = GV%H_to_kg_m2 * h_old(i,j,k)
enddo ; enddo ; enddo
dzt(:,:,:) = 1.0
call thickness_to_dz(h_old, tv, dzt, G, GV, US)
do k=1,nk ; do j=jsc,jec ; do i=isc,iec
dzt(i,j,k) = US%Z_to_m * dzt(i,j,k)
enddo ; enddo ; enddo
dz_ml(:,:) = 0.0
do j=jsc,jec ; do i=isc,iec
surface_field(i,j) = tv%S(i,j,1)
dz_ml(i,j) = US%Z_to_m * Hml(i,j)
enddo ; enddo
sosga = global_area_mean(surface_field, G, unscale=US%S_to_ppt)
mld_pha(:,:) = 0.0
if (.not.CS%mld_pha_calc) then
mld_pha(:,:) = CS%mld_pha_val
else
if (CS%mld_pha_use_delta_rho) then
call diagnoseMLDbyDensityDifference(-1, h_old, tv, CS%mld_pha_drho, G, GV, US, CS%diag, &
CS%mld_pha_href, id_ref_z=-1, id_ref_rho=-1, MLD_out=mld_pha)
elseif (CS%mld_pha_use_delta_eng) then
call diagnoseMLDbyEnergy((/-1, -1, -1/), h_old, tv, G, GV, US, (/CS%mld_pha_deng, &
CS%mld_pha_deng, CS%mld_pha_deng/), CS%diag, MLD_out=mld_pha)
endif
endif
!
!Calculate tendencies (i.e., field changes at dt) from the sources / sinks
!
if ((G%US%L_to_m == 1.0) .and. (G%US%s_to_T == 1.0) .and. (G%US%Z_to_m == 1.0) .and. &
(G%US%Q_to_J_kg == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0) .and. &
(US%C_to_degC == 1.0) .and. (US%S_to_ppt == 1.0)) then
! Avoid unnecessary copies when no unit conversion is needed.
call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
G%areaT, get_diag_time_end(CS%diag), &
optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, &
internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga, geolat=G%geolatT, &
photo_acc_dpth=mld_pha)
else
! tv%internal_heat is a null pointer unless DO_GEOTHERMAL = True,
! so we have to check and only do the scaling if it is associated.
if(associated(tv%internal_heat)) then
call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), &
optics%nbands, optics%max_wavelength_band, &
sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), &
opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), &
internal_heat=G%US%RZ_to_kg_m2*US%C_to_degC*tv%internal_heat(:,:), &
frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga, geolat=G%geolatT, &
photo_acc_dpth=mld_pha*US%Z_to_m)
else
call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), &
optics%nbands, optics%max_wavelength_band, &
sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), &
opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), &
frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga, geolat=G%geolatT, &
photo_acc_dpth=mld_pha*US%Z_to_m)
endif
endif
! This uses applyTracerBoundaryFluxesInOut to handle the change in tracer due to freshwater fluxes
! usually in ALE mode
if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
g_tracer=>CS%g_tracer_list
do
if (g_tracer_is_prog(g_tracer)) then
do k=1,nk ;do j=jsc,jec ; do i=isc,iec
h_work(i,j,k) = h_old(i,j,k)
if (g_tracer%diag_id_boundary_forcing_tend .gt. 0) then
g_tracer%boundary_forcing_tend(i,j,k) = g_tracer%field(i,j,k,1)
endif
enddo ; enddo ; enddo
call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), dt, &
fluxes, h_work, evap_CFL_limit, minimum_forcing_depth)
if (g_tracer%diag_id_boundary_forcing_tend .gt. 0) then
do k=1,nk ;do j=jsc,jec ; do i=isc,iec
g_tracer%boundary_forcing_tend(i,j,k)=G%mask2dT(i,j)*(g_tracer%field(i,j,k,1) &
- g_tracer%boundary_forcing_tend(i,j,k))/dt
enddo ; enddo ; enddo
endif
endif
!traverse the linked list till hit NULL
call g_tracer_get_next(g_tracer, g_tracer_next)
if (.NOT. associated(g_tracer_next)) exit
g_tracer=>g_tracer_next
enddo
endif
!
!Update Tr(n)%field from explicit vertical diffusion
!
! Use a tridiagonal solver to determine the concentrations after the
! surface source is applied and diapycnal advection and diffusion occurs.
if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
! Last arg is tau which is always 1 for MOM6
call generic_tracer_vertdiff_G(h_work, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1)
else
! Last arg is tau which is always 1 for MOM6
call generic_tracer_vertdiff_G(h_old, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1)
endif
! Update bottom fields after vertical processes
! Second arg is tau which is always 1 for MOM6
call generic_tracer_update_from_bottom(US%T_to_s*dt, 1, get_diag_time_end(CS%diag))
!Output diagnostics via diag_manager for all generic tracers and their fluxes
call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1)
#ifdef _USE_MOM6_DIAG
call g_tracer_set_csdiag(CS%diag)
#endif
end subroutine MOM_generic_tracer_column_physics
!> This subroutine calculates mass-weighted integral on the PE either
!! of all available tracer concentrations, or of a tracer that is
!! being requested specifically, returning the number of stocks it has
!! calculated. If the stock_index is present, only the stock corresponding
!! to that coded index is returned.
function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
type(EFP_type), dimension(:), intent(out) :: stocks !< The mass-weighted integrated amount of each
!! tracer, in kg times concentration units [kg conc]
type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module.
character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated.
character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated.
integer, optional, intent(in) :: stock_index !< The coded index of a specific stock
!! being sought.
integer :: MOM_generic_tracer_stock !< Return value, the
!! number of stocks calculated here.
! Local variables
type(g_tracer_type), pointer :: g_tracer, g_tracer_next
real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc]
real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc]
character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock'
integer :: m
MOM_generic_tracer_stock = 0
if (.not.associated(CS)) return
if (present(stock_index)) then ; if (stock_index > 0) then
! Check whether this stock is available from this routine.
! No stocks from this routine are being checked yet. Return 0.
return
endif ; endif
if (.NOT. associated(CS%g_tracer_list)) return ! No stocks.
m=1 ; g_tracer=>CS%g_tracer_list
do
call g_tracer_get_alias(g_tracer,names(m))
call g_tracer_get_values(g_tracer,names(m),'units',units(m))
units(m) = trim(units(m))//" kg"
call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field)
tr_ptr => tr_field(:,:,:,1)
stocks(m) = global_mass_int_EFP(h, G, GV, tr_ptr, on_PE_only=.true.)
!traverse the linked list till hit NULL
call g_tracer_get_next(g_tracer, g_tracer_next)
if (.NOT. associated(g_tracer_next)) exit
g_tracer=>g_tracer_next
m = m+1
enddo
MOM_generic_tracer_stock = m
end function MOM_generic_tracer_stock
!> This subroutine finds the global min and max of either of all available
!! tracer concentrations, or of a tracer that is being requested specifically,
!! returning the number of tracers it has evaluated.
!! It also optionally returns the locations of the extrema.
function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, G, CS, names, units, &
xgmin, ygmin, zgmin, xgmax, ygmax, zgmax)
integer, intent(in) :: ind_start !< The index of the tracer to start with
logical, dimension(:), intent(out) :: got_minmax !< Indicates whether the global min and
!! max are found for each tracer
real, dimension(:), intent(out) :: gmin !< Global minimum of each tracer [conc]
real, dimension(:), intent(out) :: gmax !< Global maximum of each tracer [conc]
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module.
character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated.
character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated.
real, dimension(:), optional, intent(out) :: xgmin !< The x-position of the global minimum in the
!! units of G%geoLonT, often [degrees_E] or [km] or [m]
real, dimension(:), optional, intent(out) :: ygmin !< The y-position of the global minimum in the
!! units of G%geoLatT, often [degrees_N] or [km] or [m]
real, dimension(:), optional, intent(out) :: zgmin !< The z-position of the global minimum [layer]
real, dimension(:), optional, intent(out) :: xgmax !< The x-position of the global maximum in the
!! units of G%geoLonT, often [degrees_E] or [km] or [m]
real, dimension(:), optional, intent(out) :: ygmax !< The y-position of the global maximum in the
!! units of G%geoLatT, often [degrees_N] or [km] or [m]
real, dimension(:), optional, intent(out) :: zgmax !< The z-position of the global maximum [layer]
integer :: MOM_generic_tracer_min_max !< Return value, the
!! number of tracers done here.
! Local variables
type(g_tracer_type), pointer :: g_tracer, g_tracer_next
real, dimension(:,:,:,:), pointer :: tr_field ! The tracer array whose extrema are being sought [conc]
real, dimension(:,:,:), pointer :: tr_ptr ! The tracer array whose extrema are being sought [conc]
real :: x_min ! The x-position of the global minimum in the units of G%geoLonT, often [degrees_E] or [km] or [m]
real :: y_min ! The y-position of the global minimum in the units of G%geoLatT, often [degrees_N] or [km] or [m]
real :: z_min ! The z-position of the global minimum [layer]
real :: x_max ! The x-position of the global maximum in the units of G%geoLonT, often [degrees_E] or [km] or [m]
real :: y_max ! The y-position of the global maximum in the units of G%geoLatT, often [degrees_N] or [km] or [m]
real :: z_max ! The z-position of the global maximum [layer]
character(len=128), parameter :: sub_name = 'MOM_generic_tracer_min_max'
logical :: find_location
integer :: isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau
integer :: k, is, ie, js, je, m
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
MOM_generic_tracer_min_max = 0
if (.not.associated(CS)) return
if (.NOT. associated(CS%g_tracer_list)) return ! No stocks.
call g_tracer_get_common(isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau)
find_location = present(xgmin) .or. present(ygmin) .or. present(zgmin) .or. &
present(xgmax) .or. present(ygmax) .or. present(zgmax)
m=ind_start ; g_tracer=>CS%g_tracer_list
do
call g_tracer_get_alias(g_tracer,names(m))
call g_tracer_get_values(g_tracer,names(m),'units',units(m))
call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field)
gmin(m) = -1.0
gmax(m) = -1.0
tr_ptr => tr_field(:,:,:,1)
if (find_location) then
call array_global_min_max(tr_ptr, G, nk, gmin(m), gmax(m), &
x_min, y_min, z_min, x_max, y_max, z_max)