-
Notifications
You must be signed in to change notification settings - Fork 227
/
MOM_diag_mediator.F90
4661 lines (4189 loc) · 241 KB
/
MOM_diag_mediator.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 subroutines here provide convenient wrappers to the fms diag_manager
!! interfaces with additional diagnostic capabilies.
module MOM_diag_mediator
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_checksums, only : chksum0, zchksum
use MOM_checksums, only : hchksum, uchksum, vchksum, Bchksum
use MOM_coms, only : PE_here
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE
use MOM_diag_manager_infra, only : MOM_diag_manager_init, MOM_diag_manager_end
use MOM_diag_manager_infra, only : diag_axis_init=>MOM_diag_axis_init, get_MOM_diag_axis_name
use MOM_diag_manager_infra, only : send_data_infra, MOM_diag_field_add_attribute, EAST, NORTH
use MOM_diag_manager_infra, only : register_diag_field_infra, register_static_field_infra
use MOM_diag_manager_infra, only : get_MOM_diag_field_id, DIAG_FIELD_NOT_FOUND
use MOM_diag_remap, only : diag_remap_ctrl, diag_remap_update, diag_remap_calc_hmask
use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap
use MOM_diag_remap, only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field
use MOM_diag_remap, only : horizontally_average_diag_field, diag_remap_get_axes_info
use MOM_diag_remap, only : diag_remap_configure_axes, diag_remap_axes_configured
use MOM_diag_remap, only : diag_remap_diag_registration_closed, diag_remap_set_active
use MOM_EOS, only : EOS_type
use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe, assert, callTree_showQuery
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_interface_heights, only : thickness_to_dz
use MOM_io, only : slasher, vardesc, query_vardesc, MOM_read_data
use MOM_io, only : get_filename_appendix
use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc
use MOM_string_functions, only : lowercase
use MOM_time_manager, only : time_type
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
implicit none ; private
#undef __DO_SAFETY_CHECKS__
#define IMPLIES(A, B) ((.not. (A)) .or. (B))
#define MAX_DSAMP_LEV 2
public set_axes_info, post_data, register_diag_field, time_type
public post_product_u, post_product_sum_u, post_product_v, post_product_sum_v
public set_masks_for_axes
! post_data_1d_k is a deprecated interface that can be replaced by a call to post_data, but
! it is being retained for backward compatibility to older versions of the ocean_BGC code.
public post_data_1d_k
public safe_alloc_ptr, safe_alloc_alloc
public enable_averaging, enable_averages, disable_averaging, query_averaging_enabled
public diag_mediator_init, diag_mediator_end, set_diag_mediator_grid
public diag_mediator_infrastructure_init
public diag_mediator_close_registration, get_diag_time_end
public diag_axis_init, ocean_register_diag, register_static_field
public register_scalar_field
public define_axes_group, diag_masks_set
public diag_register_area_ids
public register_cell_measure, diag_associate_volume_cell_measure
public diag_get_volume_cell_measure_dm_id
public diag_set_state_ptrs, diag_update_remap_grids
public diag_grid_storage_init, diag_grid_storage_end
public diag_copy_diag_to_storage, diag_copy_storage_to_diag
public diag_save_grids, diag_restore_grids
public found_in_diagtable
!> Make a diagnostic available for averaging or output.
interface post_data
module procedure post_data_3d, post_data_2d, post_data_1d_k, post_data_0d
end interface post_data
!> Down sample a field
interface downsample_field
module procedure downsample_field_2d, downsample_field_3d
end interface downsample_field
!> Down sample the mask of a field
interface downsample_mask
module procedure downsample_mask_2d, downsample_mask_3d
end interface downsample_mask
!> Down sample a diagnostic field
interface downsample_diag_field
module procedure downsample_diag_field_2d, downsample_diag_field_3d
end interface downsample_diag_field
!> Contained for down sampled masks
type, private :: diag_dsamp
real, pointer, dimension(:,:) :: mask2d => null() !< Mask for 2d (x-y) axes [nondim]
real, pointer, dimension(:,:,:) :: mask3d => null() !< Mask for 3d axes [nondim]
end type diag_dsamp
!> A group of 1D axes that comprise a 1D/2D/3D mesh
type, public :: axes_grp
character(len=15) :: id !< The id string for this particular combination of handles.
integer :: rank !< Number of dimensions in the list of axes.
integer, dimension(:), allocatable :: handles !< Handles to 1D axes.
type(diag_ctrl), pointer :: diag_cs => null() !< Circular link back to the main diagnostics control structure
!! (Used to avoid passing said structure into every possible call).
! ID's for cell_methods
character(len=9) :: x_cell_method = '' !< Default nature of data representation, if axes group
!! includes x-direction.
character(len=9) :: y_cell_method = '' !< Default nature of data representation, if axes group
!! includes y-direction.
character(len=9) :: v_cell_method = '' !< Default nature of data representation, if axes group
!! includes vertical direction.
! For remapping
integer :: nz = 0 !< Vertical dimension of diagnostic
integer :: vertical_coordinate_number = 0 !< Index of the corresponding diag_remap_ctrl for this axis group
! For detecting position on the grid
logical :: is_h_point = .false. !< If true, indicates that this axes group is for an h-point located field.
logical :: is_q_point = .false. !< If true, indicates that this axes group is for a q-point located field.
logical :: is_u_point = .false. !< If true, indicates that this axes group is for a u-point located field.
logical :: is_v_point = .false. !< If true, indicates that this axes group is for a v-point located field.
logical :: is_layer = .false. !< If true, indicates that this axes group is for a layer vertically-located field.
logical :: is_interface = .false. !< If true, indicates that this axes group is for an interface
!! vertically-located field.
logical :: is_native = .true. !< If true, indicates that this axes group is for a native model grid.
!! False for any other grid. Used for rank>2.
logical :: needs_remapping = .false. !< If true, indicates that this axes group is for a intensive layer-located
!! field that must be remapped to these axes. Used for rank>2.
logical :: needs_interpolating = .false. !< If true, indicates that this axes group is for a sampled
!! interface-located field that must be interpolated to
!! these axes. Used for rank>2.
integer :: downsample_level = 1 !< If greater than 1, the factor by which this diagnostic/axes/masks be downsampled
! For horizontally averaged diagnositcs (applies to 2d and 3d fields only)
type(axes_grp), pointer :: xyave_axes => null() !< The associated 1d axes for horizontall area-averaged diagnostics
! ID's for cell_measures
integer :: id_area = -1 !< The diag_manager id for area to be used for cell_measure of variables with this axes_grp.
integer :: id_volume = -1 !< The diag_manager id for volume to be used for cell_measure of variables
!! with this axes_grp.
! For masking
real, pointer, dimension(:,:) :: mask2d => null() !< Mask for 2d (x-y) axes [nondim]
real, pointer, dimension(:,:,:) :: mask3d => null() !< Mask for 3d axes [nondim]
type(diag_dsamp), dimension(2:MAX_DSAMP_LEV) :: dsamp !< Downsample container
end type axes_grp
!> Contains an array to store a diagnostic target grid
type, private :: diag_grids_type
real, dimension(:,:,:), allocatable :: h !< Target grid for remapped coordinate [H ~> m or kg m-2] or [Z ~> m]
end type diag_grids_type
!> Stores all the remapping grids and the model's native space thicknesses
type, public :: diag_grid_storage
integer :: num_diag_coords !< Number of target coordinates
real, dimension(:,:,:), allocatable :: h_state !< Layer thicknesses in native
!! space [H ~> m or kg m-2]
type(diag_grids_type), dimension(:), allocatable :: diag_grids !< Primarily empty, except h field
end type diag_grid_storage
! Integers to encode the total cell methods
!integer :: PPP=111 ! x:point,y:point,z:point, this kind of diagnostic is not currently present in diag_table.MOM6
!integer :: PPS=112 ! x:point,y:point,z:sum , this kind of diagnostic is not currently present in diag_table.MOM6
!integer :: PPM=113 ! x:point,y:point,z:mean , this kind of diagnostic is not currently present in diag_table.MOM6
integer :: PSP=121 !< x:point,y:sum,z:point
integer :: PSS=122 !< x:point,y:sum,z:point
integer :: PSM=123 !< x:point,y:sum,z:mean
integer :: PMP=131 !< x:point,y:mean,z:point
integer :: PMM=133 !< x:point,y:mean,z:mean
integer :: SPP=211 !< x:sum,y:point,z:point
integer :: SPS=212 !< x:sum,y:point,z:sum
integer :: SSP=221 !< x:sum;y:sum,z:point
integer :: MPP=311 !< x:mean,y:point,z:point
integer :: MPM=313 !< x:mean,y:point,z:mean
integer :: MMP=331 !< x:mean,y:mean,z:point
integer :: MMS=332 !< x:mean,y:mean,z:sum
integer :: SSS=222 !< x:sum,y:sum,z:sum
integer :: MMM=333 !< x:mean,y:mean,z:mean
integer :: MSK=-1 !< Use the downsample method of a mask
!> This type is used to represent a diagnostic at the diag_mediator level.
!!
!! There can be both 'primary' and 'seconday' diagnostics. The primaries
!! reside in the diag_cs%diags array. They have an id which is an index
!! into this array. The secondaries are 'variations' on the primary diagnostic.
!! For example the CMOR diagnostics are secondary. The secondary diagnostics
!! are kept in a list with the primary diagnostic as the head.
type, private :: diag_type
logical :: in_use !< True if this entry is being used.
integer :: fms_diag_id !< Underlying FMS diag_manager id.
integer :: fms_xyave_diag_id = -1 !< For a horizontally area-averaged diagnostic.
integer :: downsample_diag_id = -1 !< For a horizontally area-downsampled diagnostic.
character(64) :: debug_str = '' !< For FATAL errors and debugging.
type(axes_grp), pointer :: axes => null() !< The axis group for this diagnostic
type(diag_type), pointer :: next => null() !< Pointer to the next diagnostic
real :: conversion_factor = 0. !< If non-zero, a factor to multiply data by before posting to FMS,
!! often including factors to undo internal scaling in units of [a A-1 ~> 1]
logical :: v_extensive = .false. !< True for vertically extensive fields (vertically integrated).
!! False for intensive (concentrations).
integer :: xyz_method = 0 !< A 3 digit integer encoding the diagnostics cell method
!! It can be used to determine the downsample algorithm
end type diag_type
!> Container for down sampling information
type diagcs_dsamp
integer :: isc !< The start i-index of cell centers within the computational domain
integer :: iec !< The end i-index of cell centers within the computational domain
integer :: jsc !< The start j-index of cell centers within the computational domain
integer :: jec !< The end j-index of cell centers within the computational domain
integer :: isd !< The start i-index of cell centers within the data domain
integer :: ied !< The end i-index of cell centers within the data domain
integer :: jsd !< The start j-index of cell centers within the data domain
integer :: jed !< The end j-index of cell centers within the data domain
integer :: isg !< The start i-index of cell centers within the global domain
integer :: ieg !< The end i-index of cell centers within the global domain
integer :: jsg !< The start j-index of cell centers within the global domain
integer :: jeg !< The end j-index of cell centers within the global domain
integer :: isgB !< The start i-index of cell corners within the global domain
integer :: iegB !< The end i-index of cell corners within the global domain
integer :: jsgB !< The start j-index of cell corners within the global domain
integer :: jegB !< The end j-index of cell corners within the global domain
!>@{ Axes for each location on a diagnostic grid
type(axes_grp) :: axesBL, axesTL, axesCuL, axesCvL
type(axes_grp) :: axesBi, axesTi, axesCui, axesCvi
type(axes_grp) :: axesB1, axesT1, axesCu1, axesCv1
type(axes_grp), dimension(:), allocatable :: remap_axesTL, remap_axesBL, remap_axesCuL, remap_axesCvL
type(axes_grp), dimension(:), allocatable :: remap_axesTi, remap_axesBi, remap_axesCui, remap_axesCvi
!>@}
real, dimension(:,:), pointer :: mask2dT => null() !< 2D mask array for cell-center points [nondim]
real, dimension(:,:), pointer :: mask2dBu => null() !< 2D mask array for cell-corner points [nondim]
real, dimension(:,:), pointer :: mask2dCu => null() !< 2D mask array for east-face points [nondim]
real, dimension(:,:), pointer :: mask2dCv => null() !< 2D mask array for north-face points [nondim]
!>@{ 3D mask arrays for diagnostics at layers (mask...L) and interfaces (mask...i), all [nondim]
real, dimension(:,:,:), pointer :: mask3dTL => null()
real, dimension(:,:,:), pointer :: mask3dBL => null()
real, dimension(:,:,:), pointer :: mask3dCuL => null()
real, dimension(:,:,:), pointer :: mask3dCvL => null()
real, dimension(:,:,:), pointer :: mask3dTi => null()
real, dimension(:,:,:), pointer :: mask3dBi => null()
real, dimension(:,:,:), pointer :: mask3dCui => null()
real, dimension(:,:,:), pointer :: mask3dCvi => null()
!>@}
end type diagcs_dsamp
!> The following data type a list of diagnostic fields an their variants,
!! as well as variables that control the handling of model output.
type, public :: diag_ctrl
integer :: available_diag_doc_unit = -1 !< The unit number of a diagnostic documentation file.
!! This file is open if available_diag_doc_unit is > 0.
integer :: chksum_iounit = -1 !< The unit number of a diagnostic documentation file.
!! This file is open if available_diag_doc_unit is > 0.
logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics
logical :: show_call_tree !< Display the call tree while running. Set by VERBOSITY level.
logical :: index_space_axes !< If true, diagnostic horizontal coordinates axes are in index space.
! The following fields are used for the output of the data.
integer :: is !< The start i-index of cell centers within the computational domain
integer :: ie !< The end i-index of cell centers within the computational domain
integer :: js !< The start j-index of cell centers within the computational domain
integer :: je !< The end j-index of cell centers within the computational domain
integer :: isd !< The start i-index of cell centers within the data domain
integer :: ied !< The end i-index of cell centers within the data domain
integer :: jsd !< The start j-index of cell centers within the data domain
integer :: jed !< The end j-index of cell centers within the data domain
real :: time_int !< The time interval for any fields
!! that are offered for averaging [s].
type(time_type) :: time_end !< The end time of the valid
!! interval for any offered field.
logical :: ave_enabled = .false. !< True if averaging is enabled.
!>@{ The following are 3D and 2D axis groups defined for output. The names
!! indicate the horizontal (B, T, Cu, or Cv) and vertical (L, i, or 1) locations.
type(axes_grp) :: axesBL, axesTL, axesCuL, axesCvL
type(axes_grp) :: axesBi, axesTi, axesCui, axesCvi
type(axes_grp) :: axesB1, axesT1, axesCu1, axesCv1
!>@}
type(axes_grp) :: axesZi !< A 1-D z-space axis at interfaces
type(axes_grp) :: axesZL !< A 1-D z-space axis at layer centers
type(axes_grp) :: axesNull !< An axis group for scalars
real, dimension(:,:), pointer :: mask2dT => null() !< 2D mask array for cell-center points [nondim]
real, dimension(:,:), pointer :: mask2dBu => null() !< 2D mask array for cell-corner points [nondim]
real, dimension(:,:), pointer :: mask2dCu => null() !< 2D mask array for east-face points [nondim]
real, dimension(:,:), pointer :: mask2dCv => null() !< 2D mask array for north-face points [nondim]
!>@{ 3D mask arrays for diagnostics at layers (mask...L) and interfaces (mask...i) all [nondim]
real, dimension(:,:,:), pointer :: mask3dTL => null()
real, dimension(:,:,:), pointer :: mask3dBL => null()
real, dimension(:,:,:), pointer :: mask3dCuL => null()
real, dimension(:,:,:), pointer :: mask3dCvL => null()
real, dimension(:,:,:), pointer :: mask3dTi => null()
real, dimension(:,:,:), pointer :: mask3dBi => null()
real, dimension(:,:,:), pointer :: mask3dCui => null()
real, dimension(:,:,:), pointer :: mask3dCvi => null()
type(diagcs_dsamp), dimension(2:MAX_DSAMP_LEV) :: dsamp !< Downsample control container
!>@}
! Space for diagnostics is dynamically allocated as it is needed.
! The chunk size is how much the array should grow on each new allocation.
#define DIAG_ALLOC_CHUNK_SIZE 100
type(diag_type), dimension(:), allocatable :: diags !< The list of diagnostics
integer :: next_free_diag_id !< The next unused diagnostic ID
!> default missing value to be sent to ALL diagnostics registrations [various]
real :: missing_value = -1.0e+34
!> Number of diagnostic vertical coordinates (remapped)
integer :: num_diag_coords
!> Control structure for each possible coordinate
type(diag_remap_ctrl), dimension(:), allocatable :: diag_remap_cs
type(diag_grid_storage) :: diag_grid_temp !< Stores the remapped diagnostic grid
logical :: diag_grid_overridden = .false. !< True if the diagnostic grids have been overriden
type(axes_grp), dimension(:), allocatable :: &
remap_axesZL, & !< The 1-D z-space cell-centered axis for remapping
remap_axesZi !< The 1-D z-space interface axis for remapping
!>@{ Axes used for remapping
type(axes_grp), dimension(:), allocatable :: remap_axesTL, remap_axesBL, remap_axesCuL, remap_axesCvL
type(axes_grp), dimension(:), allocatable :: remap_axesTi, remap_axesBi, remap_axesCui, remap_axesCvi
!>@}
! Pointer to H, G and T&S needed for remapping
real, dimension(:,:,:), pointer :: h => null() !< The thicknesses needed for remapping [H ~> m or kg m-2]
real, dimension(:,:,:), pointer :: T => null() !< The temperatures needed for remapping [C ~> degC]
real, dimension(:,:,:), pointer :: S => null() !< The salinities needed for remapping [S ~> ppt]
type(EOS_type), pointer :: eqn_of_state => null() !< The equation of state type
type(thermo_var_ptrs), pointer :: tv => null() !< A sturcture with thermodynamic variables that are
!! are used to convert thicknesses to vertical extents
type(ocean_grid_type), pointer :: G => null() !< The ocean grid type
type(verticalGrid_type), pointer :: GV => null() !< The model's vertical ocean grid
type(unit_scale_type), pointer :: US => null() !< A dimensional unit scaling type
!> The volume cell measure (special diagnostic) manager id
integer :: volume_cell_measure_dm_id = -1
#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)
! Keep a copy of h so that we know whether it has changed [H ~> m or kg m-2]. If it has then
! need the target grid for vertical remapping needs to have been updated.
real, dimension(:,:,:), allocatable :: h_old
#endif
!> Number of checksum-only diagnostics
integer :: num_chksum_diags
real, dimension(:,:,:), allocatable :: h_begin !< Layer thicknesses at the beginning of the timestep used
!! for remapping of extensive variables [H ~> m or kg m-2]
end type diag_ctrl
!>@{ CPU clocks
integer :: id_clock_diag_mediator, id_clock_diag_remap, id_clock_diag_grid_updates
!>@}
contains
!> Sets up diagnostics axes
subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical)
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(param_file_type), intent(in) :: param_file !< Parameter file structure
type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
logical, optional, intent(in) :: set_vertical !< If true or missing, set up
!! vertical axes
! Local variables
integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, id_null
integer :: id_zl_native, id_zi_native
integer :: i, j, nz
real :: zlev(GV%ke) ! Numerical values for layer vertical coordinates, in unscaled units
! that might be [m], [kg m-3] or [nondim], depending on the coordinate.
real :: zinter(GV%ke+1) ! Numerical values for interface vertical coordinates, in unscaled units
! that might be [m], [kg m-3] or [nondim], depending on the coordinate.
logical :: set_vert
real, allocatable, dimension(:) :: IaxB, iax ! Index-based integer and half-integer i-axis labels [nondim]
real, allocatable, dimension(:) :: JaxB, jax ! Index-based integer and half-integer j-axis labels [nondim]
set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical
if (diag_cs%index_space_axes) then
allocate(IaxB(G%IsgB:G%IegB))
do i=G%IsgB, G%IegB
Iaxb(i)=real(i)
enddo
allocate(iax(G%isg:G%ieg))
do i=G%isg, G%ieg
iax(i)=real(i)-0.5
enddo
allocate(JaxB(G%JsgB:G%JegB))
do j=G%JsgB, G%JegB
JaxB(j)=real(j)
enddo
allocate(jax(G%jsg:G%jeg))
do j=G%jsg, G%jeg
jax(j)=real(j)-0.5
enddo
endif
! Horizontal axes for the native grids
if (G%symmetric) then
if (diag_cs%index_space_axes) then
id_xq = diag_axis_init('iq', IaxB(G%isgB:G%iegB), 'none', 'x', &
'q point grid-space longitude', G%Domain, position=EAST)
id_yq = diag_axis_init('jq', JaxB(G%jsgB:G%jegB), 'none', 'y', &
'q point grid space latitude', G%Domain, position=NORTH)
else
id_xq = diag_axis_init('xq', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', &
'q point nominal longitude', G%Domain, position=EAST)
id_yq = diag_axis_init('yq', G%gridLatB(G%jsgB:G%jegB), G%y_axis_units, 'y', &
'q point nominal latitude', G%Domain, position=NORTH)
endif
else
if (diag_cs%index_space_axes) then
id_xq = diag_axis_init('Iq', IaxB(G%isg:G%ieg), 'none', 'x', &
'q point grid-space longitude', G%Domain, position=EAST)
id_yq = diag_axis_init('Jq', JaxB(G%jsg:G%jeg), 'none', 'y', &
'q point grid space latitude', G%Domain, position=NORTH)
else
id_xq = diag_axis_init('xq', G%gridLonB(G%isg:G%ieg), G%x_axis_units, 'x', &
'q point nominal longitude', G%Domain, position=EAST)
id_yq = diag_axis_init('yq', G%gridLatB(G%jsg:G%jeg), G%y_axis_units, 'y', &
'q point nominal latitude', G%Domain, position=NORTH)
endif
endif
if (diag_cs%index_space_axes) then
id_xh = diag_axis_init('ih', iax(G%isg:G%ieg), 'none', 'x', &
'h point grid-space longitude', G%Domain)
id_yh = diag_axis_init('jh', jax(G%jsg:G%jeg), 'none', 'y', &
'h point grid space latitude', G%Domain)
else
id_xh = diag_axis_init('xh', G%gridLonT(G%isg:G%ieg), G%x_axis_units, 'x', &
'h point nominal longitude', G%Domain)
id_yh = diag_axis_init('yh', G%gridLatT(G%jsg:G%jeg), G%y_axis_units, 'y', &
'h point nominal latitude', G%Domain)
endif
if (set_vert) then
nz = GV%ke
zinter(1:nz+1) = GV%sInterface(1:nz+1)
zlev(1:nz) = GV%sLayer(1:nz)
id_zl = diag_axis_init('zl', zlev, trim(GV%zAxisUnits), 'z', &
'Layer '//trim(GV%zAxisLongName), direction=GV%direction)
id_zi = diag_axis_init('zi', zinter, trim(GV%zAxisUnits), 'z', &
'Interface '//trim(GV%zAxisLongName), direction=GV%direction)
else
id_zl = -1 ; id_zi = -1
endif
id_zl_native = id_zl ; id_zi_native = id_zi
! Vertical axes for the interfaces and layers
call define_axes_group(diag_cs, (/ id_zi /), diag_cs%axesZi, &
v_cell_method='point', is_interface=.true.)
call define_axes_group(diag_cs, (/ id_zL /), diag_cs%axesZL, &
v_cell_method='mean', is_layer=.true.)
! Axis groupings for the model layers
call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zL /), diag_cs%axesTL, &
x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &
is_h_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zL /), diag_cs%axesBL, &
x_cell_method='point', y_cell_method='point', v_cell_method='mean', &
is_q_point=.true., is_layer=.true.)
call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zL /), diag_cs%axesCuL, &
x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &
is_u_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zL /), diag_cs%axesCvL, &
x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &
is_v_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
! Axis groupings for the model interfaces
call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%axesTi, &
x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
is_h_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%axesBi, &
x_cell_method='point', y_cell_method='point', v_cell_method='point', &
is_q_point=.true., is_interface=.true.)
call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%axesCui, &
x_cell_method='point', y_cell_method='mean', v_cell_method='point', &
is_u_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%axesCvi, &
x_cell_method='mean', y_cell_method='point', v_cell_method='point', &
is_v_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
! Axis groupings for 2-D arrays
call define_axes_group(diag_cs, (/ id_xh, id_yh /), diag_cs%axesT1, &
x_cell_method='mean', y_cell_method='mean', is_h_point=.true.)
call define_axes_group(diag_cs, (/ id_xq, id_yq /), diag_cs%axesB1, &
x_cell_method='point', y_cell_method='point', is_q_point=.true.)
call define_axes_group(diag_cs, (/ id_xq, id_yh /), diag_cs%axesCu1, &
x_cell_method='point', y_cell_method='mean', is_u_point=.true.)
call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, &
x_cell_method='mean', y_cell_method='point', is_v_point=.true.)
! Axis group for special null axis from diag manager.
id_null = diag_axis_init('scalar_axis', (/0./), 'none', 'N', 'none', null_axis=.true.)
call define_axes_group(diag_cs, (/ id_null /), diag_cs%axesNull)
!Non-native Non-downsampled
if (diag_cs%num_diag_coords>0) then
allocate(diag_cs%remap_axesZL(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesTL(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesBL(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesCuL(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesCvL(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesZi(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesTi(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesBi(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesCui(diag_cs%num_diag_coords))
allocate(diag_cs%remap_axesCvi(diag_cs%num_diag_coords))
endif
do i=1, diag_cs%num_diag_coords
! For each possible diagnostic coordinate
call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), GV, US, param_file)
! Allocate these arrays since the size of the diagnostic array is now known
allocate(diag_cs%diag_remap_cs(i)%h(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz))
allocate(diag_cs%diag_remap_cs(i)%h_extensive(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz))
! This vertical coordinate has been configured so can be used.
if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then
! This fetches the 1D-axis id for layers and interfaces and overwrite
! id_zl and id_zi from above. It also returns the number of layers.
call diag_remap_get_axes_info(diag_cs%diag_remap_cs(i), nz, id_zL, id_zi)
! Axes for z layers
call define_axes_group(diag_cs, (/ id_zL /), diag_cs%remap_axesZL(i), &
nz=nz, vertical_coordinate_number=i, &
v_cell_method='mean', &
is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.)
call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zL /), diag_cs%remap_axesTL(i), &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &
is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
xyave_axes=diag_cs%remap_axesZL(i))
!! \note Remapping for B points is not yet implemented so needs_remapping is not
!! provided for remap_axesBL
call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zL /), diag_cs%remap_axesBL(i), &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='point', y_cell_method='point', v_cell_method='mean', &
is_q_point=.true., is_layer=.true., is_native=.false.)
call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zL /), diag_cs%remap_axesCuL(i), &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &
is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
xyave_axes=diag_cs%remap_axesZL(i))
call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zL /), diag_cs%remap_axesCvL(i), &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &
is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
xyave_axes=diag_cs%remap_axesZL(i))
! Axes for z interfaces
call define_axes_group(diag_cs, (/ id_zi /), diag_cs%remap_axesZi(i), &
nz=nz, vertical_coordinate_number=i, &
v_cell_method='point', &
is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true.)
call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%remap_axesTi(i), &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true., &
xyave_axes=diag_cs%remap_axesZi(i))
!! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBi
call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%remap_axesBi(i), &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='point', y_cell_method='point', v_cell_method='point', &
is_q_point=.true., is_interface=.true., is_native=.false.)
call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%remap_axesCui(i), &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='point', y_cell_method='mean', v_cell_method='point', &
is_u_point=.true., is_interface=.true., is_native=.false., &
needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))
call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%remap_axesCvi(i), &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='mean', y_cell_method='point', v_cell_method='point', &
is_v_point=.true., is_interface=.true., is_native=.false., &
needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))
endif
enddo
if (diag_cs%index_space_axes) then
deallocate(IaxB, iax, JaxB, jax)
endif
!Define the downsampled axes
call set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_native)
call diag_grid_storage_init(diag_CS%diag_grid_temp, G, GV, diag_CS)
end subroutine set_axes_info
subroutine set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_native)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
type(param_file_type), intent(in) :: param_file !< Parameter file structure
type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
integer, intent(in) :: id_zl_native !< ID of native layers
integer, intent(in) :: id_zi_native !< ID of native interfaces
! Local variables
integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh
integer :: i, j, nz, dl
real, dimension(:), pointer :: gridLonT_dsamp =>NULL() ! The longitude of downsampled T points for labeling
! the output axes, often in units of [degrees_N] or
! [km] or [m] or [gridpoints].
real, dimension(:), pointer :: gridLatT_dsamp =>NULL() ! The latitude of downsampled T points for labeling
! the output axes, often in units of [degrees_N] or
! [km] or [m] or [gridpoints].
real, dimension(:), pointer :: gridLonB_dsamp =>NULL() ! The longitude of downsampled B points for labeling
! the output axes, often in units of [degrees_N] or
! [km] or [m] or [gridpoints].
real, dimension(:), pointer :: gridLatB_dsamp =>NULL() ! The latitude of downsampled B points for labeling
! the output axes, often in units of [degrees_N] or
! [km] or [m] or [gridpoints].
id_zl = id_zl_native ; id_zi = id_zi_native
!Axes group for native downsampled diagnostics
do dl=2,MAX_DSAMP_LEV
if (dl /= 2) call MOM_error(FATAL, "set_axes_info_dsamp: Downsample level other than 2 is not supported yet!")
if (G%symmetric) then
allocate(gridLonB_dsamp(diag_cs%dsamp(dl)%isgB:diag_cs%dsamp(dl)%iegB))
allocate(gridLatB_dsamp(diag_cs%dsamp(dl)%jsgB:diag_cs%dsamp(dl)%jegB))
do i=diag_cs%dsamp(dl)%isgB,diag_cs%dsamp(dl)%iegB; gridLonB_dsamp(i) = G%gridLonB(G%isgB+dl*i); enddo
do j=diag_cs%dsamp(dl)%jsgB,diag_cs%dsamp(dl)%jegB; gridLatB_dsamp(j) = G%gridLatB(G%jsgB+dl*j); enddo
id_xq = diag_axis_init('xq', gridLonB_dsamp, G%x_axis_units, 'x', &
'q point nominal longitude', G%Domain, coarsen=2)
id_yq = diag_axis_init('yq', gridLatB_dsamp, G%y_axis_units, 'y', &
'q point nominal latitude', G%Domain, coarsen=2)
deallocate(gridLonB_dsamp,gridLatB_dsamp)
else
allocate(gridLonB_dsamp(diag_cs%dsamp(dl)%isg:diag_cs%dsamp(dl)%ieg))
allocate(gridLatB_dsamp(diag_cs%dsamp(dl)%jsg:diag_cs%dsamp(dl)%jeg))
do i=diag_cs%dsamp(dl)%isg,diag_cs%dsamp(dl)%ieg; gridLonB_dsamp(i) = G%gridLonB(G%isg+dl*i-2); enddo
do j=diag_cs%dsamp(dl)%jsg,diag_cs%dsamp(dl)%jeg; gridLatB_dsamp(j) = G%gridLatB(G%jsg+dl*j-2); enddo
id_xq = diag_axis_init('xq', gridLonB_dsamp, G%x_axis_units, 'x', &
'q point nominal longitude', G%Domain, coarsen=2)
id_yq = diag_axis_init('yq', gridLatB_dsamp, G%y_axis_units, 'y', &
'q point nominal latitude', G%Domain, coarsen=2)
deallocate(gridLonB_dsamp,gridLatB_dsamp)
endif
allocate(gridLonT_dsamp(diag_cs%dsamp(dl)%isg:diag_cs%dsamp(dl)%ieg))
allocate(gridLatT_dsamp(diag_cs%dsamp(dl)%jsg:diag_cs%dsamp(dl)%jeg))
do i=diag_cs%dsamp(dl)%isg,diag_cs%dsamp(dl)%ieg; gridLonT_dsamp(i) = G%gridLonT(G%isg+dl*i-2); enddo
do j=diag_cs%dsamp(dl)%jsg,diag_cs%dsamp(dl)%jeg; gridLatT_dsamp(j) = G%gridLatT(G%jsg+dl*j-2); enddo
id_xh = diag_axis_init('xh', gridLonT_dsamp, G%x_axis_units, 'x', &
'h point nominal longitude', G%Domain, coarsen=2)
id_yh = diag_axis_init('yh', gridLatT_dsamp, G%y_axis_units, 'y', &
'h point nominal latitude', G%Domain, coarsen=2)
deallocate(gridLonT_dsamp,gridLatT_dsamp)
! Axis groupings for the model layers
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zL /), diag_cs%dsamp(dl)%axesTL, dl, &
x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &
is_h_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zL /), diag_cs%dsamp(dl)%axesBL, dl, &
x_cell_method='point', y_cell_method='point', v_cell_method='mean', &
is_q_point=.true., is_layer=.true.)
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zL /), diag_cs%dsamp(dl)%axesCuL, dl, &
x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &
is_u_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zL /), diag_cs%dsamp(dl)%axesCvL, dl, &
x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &
is_v_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
! Axis groupings for the model interfaces
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%dsamp(dl)%axesTi, dl, &
x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
is_h_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%dsamp(dl)%axesBi, dl, &
x_cell_method='point', y_cell_method='point', v_cell_method='point', &
is_q_point=.true., is_interface=.true.)
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%dsamp(dl)%axesCui, dl, &
x_cell_method='point', y_cell_method='mean', v_cell_method='point', &
is_u_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%dsamp(dl)%axesCvi, dl, &
x_cell_method='mean', y_cell_method='point', v_cell_method='point', &
is_v_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
! Axis groupings for 2-D arrays
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh /), diag_cs%dsamp(dl)%axesT1, dl, &
x_cell_method='mean', y_cell_method='mean', is_h_point=.true.)
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq /), diag_cs%dsamp(dl)%axesB1, dl, &
x_cell_method='point', y_cell_method='point', is_q_point=.true.)
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh /), diag_cs%dsamp(dl)%axesCu1, dl, &
x_cell_method='point', y_cell_method='mean', is_u_point=.true.)
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq /), diag_cs%dsamp(dl)%axesCv1, dl, &
x_cell_method='mean', y_cell_method='point', is_v_point=.true.)
!Non-native axes
if (diag_cs%num_diag_coords>0) then
allocate(diag_cs%dsamp(dl)%remap_axesTL(diag_cs%num_diag_coords))
allocate(diag_cs%dsamp(dl)%remap_axesBL(diag_cs%num_diag_coords))
allocate(diag_cs%dsamp(dl)%remap_axesCuL(diag_cs%num_diag_coords))
allocate(diag_cs%dsamp(dl)%remap_axesCvL(diag_cs%num_diag_coords))
allocate(diag_cs%dsamp(dl)%remap_axesTi(diag_cs%num_diag_coords))
allocate(diag_cs%dsamp(dl)%remap_axesBi(diag_cs%num_diag_coords))
allocate(diag_cs%dsamp(dl)%remap_axesCui(diag_cs%num_diag_coords))
allocate(diag_cs%dsamp(dl)%remap_axesCvi(diag_cs%num_diag_coords))
endif
do i=1, diag_cs%num_diag_coords
! For each possible diagnostic coordinate
!call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), GV, param_file)
! This vertical coordinate has been configured so can be used.
if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then
! This fetches the 1D-axis id for layers and interfaces and overwrite
! id_zl and id_zi from above. It also returns the number of layers.
call diag_remap_get_axes_info(diag_cs%diag_remap_cs(i), nz, id_zL, id_zi)
! Axes for z layers
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zL /), diag_cs%dsamp(dl)%remap_axesTL(i), dl, &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &
is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
xyave_axes=diag_cs%remap_axesZL(i))
!! \note Remapping for B points is not yet implemented so needs_remapping is not
!! provided for remap_axesBL
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zL /), diag_cs%dsamp(dl)%remap_axesBL(i), dl, &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='point', y_cell_method='point', v_cell_method='mean', &
is_q_point=.true., is_layer=.true., is_native=.false.)
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zL /), diag_cs%dsamp(dl)%remap_axesCuL(i), dl, &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &
is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
xyave_axes=diag_cs%remap_axesZL(i))
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zL /), diag_cs%dsamp(dl)%remap_axesCvL(i), dl, &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &
is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
xyave_axes=diag_cs%remap_axesZL(i))
! Axes for z interfaces
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%dsamp(dl)%remap_axesTi(i), dl, &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true., &
xyave_axes=diag_cs%remap_axesZi(i))
!! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBi
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%dsamp(dl)%remap_axesBi(i), dl, &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='point', y_cell_method='point', v_cell_method='point', &
is_q_point=.true., is_interface=.true., is_native=.false.)
call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%dsamp(dl)%remap_axesCui(i), dl, &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='point', y_cell_method='mean', v_cell_method='point', &
is_u_point=.true., is_interface=.true., is_native=.false., &
needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))
call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%dsamp(dl)%remap_axesCvi(i), dl, &
nz=nz, vertical_coordinate_number=i, &
x_cell_method='mean', y_cell_method='point', v_cell_method='point', &
is_v_point=.true., is_interface=.true., is_native=.false., &
needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))
endif
enddo
enddo
end subroutine set_axes_info_dsamp
!> set_masks_for_axes sets up the 2d and 3d masks for diagnostics using the current grid
!! recorded after calling diag_update_remap_grids()
subroutine set_masks_for_axes(G, diag_cs)
type(ocean_grid_type), target, intent(in) :: G !< The ocean grid type.
type(diag_ctrl), pointer :: diag_cs !< A pointer to a type with many variables
!! used for diagnostics
! Local variables
integer :: c, nk, i, j, k
type(axes_grp), pointer :: axes => NULL(), h_axes => NULL() ! Current axes, for convenience
do c=1, diag_cs%num_diag_coords
! This vertical coordinate has been configured so can be used.
if (diag_remap_axes_configured(diag_cs%diag_remap_cs(c))) then
! Level/layer h-points in diagnostic coordinate
axes => diag_cs%remap_axesTL(c)
nk = axes%nz
allocate( axes%mask3d(G%isd:G%ied,G%jsd:G%jed,nk), source=0. )
call diag_remap_calc_hmask(diag_cs%diag_remap_cs(c), G, axes%mask3d)
h_axes => diag_cs%remap_axesTL(c) ! Use the h-point masks to generate the u-, v- and q- masks
! Level/layer u-points in diagnostic coordinate
axes => diag_cs%remap_axesCuL(c)
call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at u-layers')
call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
allocate( axes%mask3d(G%IsdB:G%IedB,G%jsd:G%jed,nk), source=0. )
do k = 1, nk ; do j=G%jsc,G%jec ; do I=G%isc-1,G%iec
if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j,k) > 0.) axes%mask3d(I,j,k) = 1.
enddo ; enddo ; enddo
! Level/layer v-points in diagnostic coordinate
axes => diag_cs%remap_axesCvL(c)
call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at v-layers')
call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
allocate( axes%mask3d(G%isd:G%ied,G%JsdB:G%JedB,nk), source=0. )
do k = 1, nk ; do J=G%jsc-1,G%jec ; do i=G%isc,G%iec
if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,J,k) = 1.
enddo ; enddo ; enddo
! Level/layer q-points in diagnostic coordinate
axes => diag_cs%remap_axesBL(c)
call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at q-layers')
call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
allocate( axes%mask3d(G%IsdB:G%IedB,G%JsdB:G%JedB,nk), source=0. )
do k = 1, nk ; do J=G%jsc-1,G%jec ; do I=G%isc-1,G%iec
if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j+1,k) + &
h_axes%mask3d(i+1,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(I,J,k) = 1.
enddo ; enddo ; enddo
! Interface h-points in diagnostic coordinate (w-point)
axes => diag_cs%remap_axesTi(c)
call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at h-interfaces')
call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
allocate( axes%mask3d(G%isd:G%ied,G%jsd:G%jed,nk+1), source=0. )
do J=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
if (h_axes%mask3d(i,j,1) > 0.) axes%mask3d(i,J,1) = 1.
do K = 2, nk
if (h_axes%mask3d(i,j,k-1) + h_axes%mask3d(i,j,k) > 0.) axes%mask3d(i,J,k) = 1.
enddo
if (h_axes%mask3d(i,j,nk) > 0.) axes%mask3d(i,J,nk+1) = 1.
enddo ; enddo
h_axes => diag_cs%remap_axesTi(c) ! Use the w-point masks to generate the u-, v- and q- masks
! Interface u-points in diagnostic coordinate
axes => diag_cs%remap_axesCui(c)
call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at u-interfaces')
call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
allocate( axes%mask3d(G%IsdB:G%IedB,G%jsd:G%jed,nk+1), source=0. )
do k = 1, nk+1 ; do j=G%jsc,G%jec ; do I=G%isc-1,G%iec
if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j,k) > 0.) axes%mask3d(I,j,k) = 1.
enddo ; enddo ; enddo
! Interface v-points in diagnostic coordinate
axes => diag_cs%remap_axesCvi(c)
call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at v-interfaces')
call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
allocate( axes%mask3d(G%isd:G%ied,G%JsdB:G%JedB,nk+1), source=0. )
do k = 1, nk+1 ; do J=G%jsc-1,G%jec ; do i=G%isc,G%iec
if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,J,k) = 1.
enddo ; enddo ; enddo
! Interface q-points in diagnostic coordinate
axes => diag_cs%remap_axesBi(c)
call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at q-interfaces')
call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
allocate( axes%mask3d(G%IsdB:G%IedB,G%JsdB:G%JedB,nk+1), source=0. )
do k = 1, nk ; do J=G%jsc-1,G%jec ; do I=G%isc-1,G%iec
if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j+1,k) + &
h_axes%mask3d(i+1,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(I,J,k) = 1.
enddo ; enddo ; enddo
endif
enddo
!Allocate and initialize the downsampled masks for the axes
call set_masks_for_axes_dsamp(G, diag_cs)
end subroutine set_masks_for_axes
subroutine set_masks_for_axes_dsamp(G, diag_cs)
type(ocean_grid_type), target, intent(in) :: G !< The ocean grid type.
type(diag_ctrl), pointer :: diag_cs !< A pointer to a type with many variables
!! used for diagnostics
! Local variables
integer :: c, dl
type(axes_grp), pointer :: axes => NULL() ! Current axes, for convenience
!Each downsampled axis needs both downsampled and non-downsampled mask
!The downsampled mask is needed for sending out the diagnostics output via diag_manager
!The non-downsampled mask is needed for downsampling the diagnostics field
do dl=2,MAX_DSAMP_LEV
if (dl /= 2) call MOM_error(FATAL, "set_masks_for_axes_dsamp: Downsample level other than 2 is not supported!")
do c=1, diag_cs%num_diag_coords
! Level/layer h-points in diagnostic coordinate
axes => diag_cs%remap_axesTL(c)
call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTL(c)%dsamp(dl)%mask3d, &
dl, G%isc, G%jsc, G%isd, G%jsd, &
G%HId2%isc, G%HId2%iec, G%HId2%jsc, G%HId2%jec, G%HId2%isd, G%HId2%ied, G%HId2%jsd, G%HId2%jed)
diag_cs%dsamp(dl)%remap_axesTL(c)%mask3d => axes%mask3d !set non-downsampled mask
! Level/layer u-points in diagnostic coordinate
axes => diag_cs%remap_axesCuL(c)
call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCuL(c)%dsamp(dl)%mask3d, &
dl, G%IscB, G%jsc, G%IsdB, G%jsd, &
G%HId2%IscB, G%HId2%IecB, G%HId2%jsc, G%HId2%jec, G%HId2%IsdB, G%HId2%IedB, G%HId2%jsd, G%HId2%jed)
diag_cs%dsamp(dl)%remap_axesCul(c)%mask3d => axes%mask3d !set non-downsampled mask
! Level/layer v-points in diagnostic coordinate
axes => diag_cs%remap_axesCvL(c)
call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvL(c)%dsamp(dl)%mask3d, &
dl, G%isc, G%JscB, G%isd, G%JsdB, &
G%HId2%isc, G%HId2%iec, G%HId2%JscB, G%HId2%JecB, G%HId2%isd, G%HId2%ied, G%HId2%JsdB, G%HId2%JedB)
diag_cs%dsamp(dl)%remap_axesCvL(c)%mask3d => axes%mask3d !set non-downsampled mask
! Level/layer q-points in diagnostic coordinate
axes => diag_cs%remap_axesBL(c)
call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBL(c)%dsamp(dl)%mask3d, &
dl, G%IscB, G%JscB, G%IsdB, G%JsdB, &
G%HId2%IscB, G%HId2%IecB, G%HId2%JscB, G%HId2%JecB, G%HId2%IsdB, G%HId2%IedB, G%HId2%JsdB, G%HId2%JedB)
diag_cs%dsamp(dl)%remap_axesBL(c)%mask3d => axes%mask3d !set non-downsampled mask
! Interface h-points in diagnostic coordinate (w-point)
axes => diag_cs%remap_axesTi(c)
call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTi(c)%dsamp(dl)%mask3d, &
dl, G%isc, G%jsc, G%isd, G%jsd, &
G%HId2%isc, G%HId2%iec, G%HId2%jsc, G%HId2%jec, G%HId2%isd, G%HId2%ied, G%HId2%jsd, G%HId2%jed)
diag_cs%dsamp(dl)%remap_axesTi(c)%mask3d => axes%mask3d !set non-downsampled mask
! Interface u-points in diagnostic coordinate
axes => diag_cs%remap_axesCui(c)
call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCui(c)%dsamp(dl)%mask3d, &
dl, G%IscB, G%jsc, G%IsdB, G%jsd, &
G%HId2%IscB, G%HId2%IecB, G%HId2%jsc, G%HId2%jec, G%HId2%IsdB, G%HId2%IedB, G%HId2%jsd, G%HId2%jed)
diag_cs%dsamp(dl)%remap_axesCui(c)%mask3d => axes%mask3d !set non-downsampled mask
! Interface v-points in diagnostic coordinate
axes => diag_cs%remap_axesCvi(c)
call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvi(c)%dsamp(dl)%mask3d, &
dl, G%isc, G%JscB, G%isd, G%JsdB, &
G%HId2%isc, G%HId2%iec, G%HId2%JscB, G%HId2%JecB, G%HId2%isd, G%HId2%ied, G%HId2%JsdB, G%HId2%JedB)
diag_cs%dsamp(dl)%remap_axesCvi(c)%mask3d => axes%mask3d !set non-downsampled mask
! Interface q-points in diagnostic coordinate
axes => diag_cs%remap_axesBi(c)
call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBi(c)%dsamp(dl)%mask3d, &
dl, G%IscB, G%JscB, G%IsdB, G%JsdB, &
G%HId2%IscB, G%HId2%IecB, G%HId2%JscB, G%HId2%JecB, G%HId2%IsdB, G%HId2%IedB, G%HId2%JsdB, G%HId2%JedB)
diag_cs%dsamp(dl)%remap_axesBi(c)%mask3d => axes%mask3d !set non-downsampled mask
enddo
enddo
end subroutine set_masks_for_axes_dsamp
!> Attaches the id of cell areas to axes groups for use with cell_measures
subroutine diag_register_area_ids(diag_cs, id_area_t, id_area_q)
type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
integer, optional, intent(in) :: id_area_t !< Diag_mediator id for area of h-cells
integer, optional, intent(in) :: id_area_q !< Diag_mediator id for area of q-cells
! Local variables
integer :: fms_id, i
if (present(id_area_t)) then
fms_id = diag_cs%diags(id_area_t)%fms_diag_id
diag_cs%axesT1%id_area = fms_id
diag_cs%axesTi%id_area = fms_id
diag_cs%axesTL%id_area = fms_id
do i=1, diag_cs%num_diag_coords
diag_cs%remap_axesTL(i)%id_area = fms_id
diag_cs%remap_axesTi(i)%id_area = fms_id
enddo
endif
if (present(id_area_q)) then
fms_id = diag_cs%diags(id_area_q)%fms_diag_id
diag_cs%axesB1%id_area = fms_id
diag_cs%axesBi%id_area = fms_id
diag_cs%axesBL%id_area = fms_id
do i=1, diag_cs%num_diag_coords
diag_cs%remap_axesBL(i)%id_area = fms_id
diag_cs%remap_axesBi(i)%id_area = fms_id
enddo
endif
end subroutine diag_register_area_ids
!> Sets a handle inside diagnostics mediator to associate 3d cell measures
subroutine register_cell_measure(G, diag, Time)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(diag_ctrl), target, intent(inout) :: diag !< Regulates diagnostic output
type(time_type), intent(in) :: Time !< Model time
! Local variables
integer :: id
id = register_diag_field('ocean_model', 'volcello', diag%axesTL, &
Time, 'Ocean grid-cell volume', 'm3', &
standard_name='ocean_volume', v_extensive=.true., &
x_cell_method='sum', y_cell_method='sum')
call diag_associate_volume_cell_measure(diag, id)
end subroutine register_cell_measure
!> Attaches the id of cell volumes to axes groups for use with cell_measures
subroutine diag_associate_volume_cell_measure(diag_cs, id_h_volume)
type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
integer, intent(in) :: id_h_volume !< Diag_manager id for volume of h-cells
! Local variables
type(diag_type), pointer :: tmp => NULL()
if (id_h_volume<=0) return ! Do nothing
diag_cs%volume_cell_measure_dm_id = id_h_volume ! Record for diag_get_volume_cell_measure_dm_id()
! Set the cell measure for this axes group to the FMS id in this coordinate system
diag_cs%diags(id_h_volume)%axes%id_volume = diag_cs%diags(id_h_volume)%fms_diag_id
tmp => diag_cs%diags(id_h_volume)%next ! First item in the list, if any
do while (associated(tmp))
! Set the cell measure for this axes group to the FMS id in this coordinate system
tmp%axes%id_volume = tmp%fms_diag_id
tmp => tmp%next ! Move to next axes group for this field
enddo