MOM6
MOM_diag_mediator.F90
Go to the documentation of this file.
1 !> The subroutines here provide convenient wrappers to the fms diag_manager
2 !! interfaces with additional diagnostic capabilies.
4 
5 ! This file is part of MOM6. See LICENSE.md for the license.
6 
7 use mom_checksums, only : chksum0, zchksum
9 use mom_coms, only : pe_here
10 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
11 use mom_cpu_clock, only : clock_module, clock_routine
12 use mom_error_handler, only : mom_error, fatal, warning, is_root_pe, assert
14 use mom_grid, only : ocean_grid_type
15 use mom_io, only : slasher, vardesc, query_vardesc, mom_read_data
16 use mom_io, only : get_filename_appendix
19 use mom_time_manager, only : time_type
22 use mom_eos, only : eos_type
32 
33 use diag_axis_mod, only : get_diag_axis_name
34 use diag_data_mod, only : null_axis_id
35 use diag_manager_mod, only : diag_manager_init, diag_manager_end
36 use diag_manager_mod, only : send_data, diag_axis_init, diag_field_add_attribute
37 ! The following module is needed for PGI since the following line does not compile with PGI 6.5.0
38 ! was: use diag_manager_mod, only : register_diag_field_fms=>register_diag_field
40 use diag_manager_mod, only : register_static_field_fms=>register_static_field
41 use diag_manager_mod, only : get_diag_field_id_fms=>get_diag_field_id
42 use diag_manager_mod, only : diag_field_not_found
43 
44 implicit none ; private
45 
46 #undef __DO_SAFETY_CHECKS__
47 #define IMPLIES(A, B) ((.not. (A)) .or. (B))
48 #define MAX_DSAMP_LEV 2
49 
51 public set_masks_for_axes
52 public post_data_1d_k
58 public diag_axis_init, ocean_register_diag, register_static_field
68 
69 !> Make a diagnostic available for averaging or output.
70 interface post_data
72 end interface post_data
73 
74 !> Down sample a field
77 end interface downsample_field
78 
79 !> Down sample the mask of a field
80 interface downsample_mask
81  module procedure downsample_mask_2d, downsample_mask_3d
82 end interface downsample_mask
83 
84 !> Down sample a diagnostic field
87 end interface downsample_diag_field
88 
89 !> Contained for down sampled masks
90 type, private :: diag_dsamp
91  real, pointer, dimension(:,:) :: mask2d => null() !< Mask for 2d (x-y) axes
92  real, pointer, dimension(:,:,:) :: mask3d => null() !< Mask for 3d axes
93 end type diag_dsamp
94 
95 !> A group of 1D axes that comprise a 1D/2D/3D mesh
96 type, public :: axes_grp
97  character(len=15) :: id !< The id string for this particular combination of handles.
98  integer :: rank !< Number of dimensions in the list of axes.
99  integer, dimension(:), allocatable :: handles !< Handles to 1D axes.
100  type(diag_ctrl), pointer :: diag_cs => null() !< Circular link back to the main diagnostics control structure
101  !! (Used to avoid passing said structure into every possible call).
102  ! ID's for cell_methods
103  character(len=9) :: x_cell_method = '' !< Default nature of data representation, if axes group
104  !! includes x-direction.
105  character(len=9) :: y_cell_method = '' !< Default nature of data representation, if axes group
106  !! includes y-direction.
107  character(len=9) :: v_cell_method = '' !< Default nature of data representation, if axes group
108  !! includes vertical direction.
109  ! For remapping
110  integer :: nz = 0 !< Vertical dimension of diagnostic
111  integer :: vertical_coordinate_number = 0 !< Index of the corresponding diag_remap_ctrl for this axis group
112  ! For detecting position on the grid
113  logical :: is_h_point = .false. !< If true, indicates that this axes group is for an h-point located field.
114  logical :: is_q_point = .false. !< If true, indicates that this axes group is for a q-point located field.
115  logical :: is_u_point = .false. !< If true, indicates that this axes group is for a u-point located field.
116  logical :: is_v_point = .false. !< If true, indicates that this axes group is for a v-point located field.
117  logical :: is_layer = .false. !< If true, indicates that this axes group is for a layer vertically-located field.
118  logical :: is_interface = .false. !< If true, indicates that this axes group is for an interface
119  !! vertically-located field.
120  logical :: is_native = .true. !< If true, indicates that this axes group is for a native model grid.
121  !! False for any other grid. Used for rank>2.
122  logical :: needs_remapping = .false. !< If true, indicates that this axes group is for a intensive layer-located
123  !! field that must be remapped to these axes. Used for rank>2.
124  logical :: needs_interpolating = .false. !< If true, indicates that this axes group is for a sampled
125  !! interface-located field that must be interpolated to
126  !! these axes. Used for rank>2.
127  integer :: downsample_level = 1 !< If greater than 1, the factor by which this diagnostic/axes/masks be downsampled
128  ! For horizontally averaged diagnositcs (applies to 2d and 3d fields only)
129  type(axes_grp), pointer :: xyave_axes => null() !< The associated 1d axes for horizontall area-averaged diagnostics
130  ! ID's for cell_measures
131  integer :: id_area = -1 !< The diag_manager id for area to be used for cell_measure of variables with this axes_grp.
132  integer :: id_volume = -1 !< The diag_manager id for volume to be used for cell_measure of variables
133  !! with this axes_grp.
134  ! For masking
135  real, pointer, dimension(:,:) :: mask2d => null() !< Mask for 2d (x-y) axes
136  real, pointer, dimension(:,:,:) :: mask3d => null() !< Mask for 3d axes
137  type(diag_dsamp), dimension(2:MAX_DSAMP_LEV) :: dsamp !< Downsample container
138 end type axes_grp
139 
140 !> Contains an array to store a diagnostic target grid
141 type, private :: diag_grids_type
142  real, dimension(:,:,:), allocatable :: h !< Target grid for remapped coordinate
143 end type diag_grids_type
144 
145 !> Stores all the remapping grids and the model's native space thicknesses
146 type, public :: diag_grid_storage
147  integer :: num_diag_coords !< Number of target coordinates
148  real, dimension(:,:,:), allocatable :: h_state !< Layer thicknesses in native space
149  type(diag_grids_type), dimension(:), allocatable :: diag_grids !< Primarily empty, except h field
150 end type diag_grid_storage
151 
152 ! Integers to encode the total cell methods
153 !integer :: PPP=111 ! x:point,y:point,z:point, this kind of diagnostic is not currently present in diag_table.MOM6
154 !integer :: PPS=112 ! x:point,y:point,z:sum , this kind of diagnostic is not currently present in diag_table.MOM6
155 !integer :: PPM=113 ! x:point,y:point,z:mean , this kind of diagnostic is not currently present in diag_table.MOM6
156 integer :: psp=121 !< x:point,y:sum,z:point
157 integer :: pss=122 !< x:point,y:sum,z:point
158 integer :: psm=123 !< x:point,y:sum,z:mean
159 integer :: pmp=131 !< x:point,y:mean,z:point
160 integer :: pmm=133 !< x:point,y:mean,z:mean
161 integer :: spp=211 !< x:sum,y:point,z:point
162 integer :: sps=212 !< x:sum,y:point,z:sum
163 integer :: ssp=221 !< x:sum;y:sum,z:point
164 integer :: mpp=311 !< x:mean,y:point,z:point
165 integer :: mpm=313 !< x:mean,y:point,z:mean
166 integer :: mmp=331 !< x:mean,y:mean,z:point
167 integer :: mms=332 !< x:mean,y:mean,z:sum
168 integer :: sss=222 !< x:sum,y:sum,z:sum
169 integer :: mmm=333 !< x:mean,y:mean,z:mean
170 integer :: msk=-1 !< Use the downsample method of a mask
171 
172 !> This type is used to represent a diagnostic at the diag_mediator level.
173 !!
174 !! There can be both 'primary' and 'seconday' diagnostics. The primaries
175 !! reside in the diag_cs%diags array. They have an id which is an index
176 !! into this array. The secondaries are 'variations' on the primary diagnostic.
177 !! For example the CMOR diagnostics are secondary. The secondary diagnostics
178 !! are kept in a list with the primary diagnostic as the head.
179 type, private :: diag_type
180  logical :: in_use !< True if this entry is being used.
181  integer :: fms_diag_id !< Underlying FMS diag_manager id.
182  integer :: fms_xyave_diag_id = -1 !< For a horizontally area-averaged diagnostic.
183  integer :: downsample_diag_id = -1 !< For a horizontally area-downsampled diagnostic.
184  character(64) :: debug_str = '' !< For FATAL errors and debugging.
185  type(axes_grp), pointer :: axes => null() !< The axis group for this diagnostic
186  type(diag_type), pointer :: next => null() !< Pointer to the next diagnostic
187  real :: conversion_factor = 0. !< A factor to multiply data by before posting to FMS, if non-zero.
188  logical :: v_extensive = .false. !< True for vertically extensive fields (vertically integrated).
189  !! False for intensive (concentrations).
190  integer :: xyz_method = 0 !< A 3 digit integer encoding the diagnostics cell method
191  !! It can be used to determine the downsample algorithm
192 end type diag_type
193 
194 !> Container for down sampling information
196  integer :: isc !< The start i-index of cell centers within the computational domain
197  integer :: iec !< The end i-index of cell centers within the computational domain
198  integer :: jsc !< The start j-index of cell centers within the computational domain
199  integer :: jec !< The end j-index of cell centers within the computational domain
200  integer :: isd !< The start i-index of cell centers within the data domain
201  integer :: ied !< The end i-index of cell centers within the data domain
202  integer :: jsd !< The start j-index of cell centers within the data domain
203  integer :: jed !< The end j-index of cell centers within the data domain
204  integer :: isg !< The start i-index of cell centers within the global domain
205  integer :: ieg !< The end i-index of cell centers within the global domain
206  integer :: jsg !< The start j-index of cell centers within the global domain
207  integer :: jeg !< The end j-index of cell centers within the global domain
208  integer :: isgb !< The start i-index of cell corners within the global domain
209  integer :: iegb !< The end i-index of cell corners within the global domain
210  integer :: jsgb !< The start j-index of cell corners within the global domain
211  integer :: jegb !< The end j-index of cell corners within the global domain
212 
213  !>@{ Axes for each location on a diagnostic grid
214  type(axes_grp) :: axesbl, axestl, axescul, axescvl
215  type(axes_grp) :: axesbi, axesti, axescui, axescvi
216  type(axes_grp) :: axesb1, axest1, axescu1, axescv1
217  type(axes_grp), dimension(:), allocatable :: remap_axestl, remap_axesbl, remap_axescul, remap_axescvl
218  type(axes_grp), dimension(:), allocatable :: remap_axesti, remap_axesbi, remap_axescui, remap_axescvi
219  !!@}
220 
221  real, dimension(:,:), pointer :: mask2dt => null() !< 2D mask array for cell-center points
222  real, dimension(:,:), pointer :: mask2dbu => null() !< 2D mask array for cell-corner points
223  real, dimension(:,:), pointer :: mask2dcu => null() !< 2D mask array for east-face points
224  real, dimension(:,:), pointer :: mask2dcv => null() !< 2D mask array for north-face points
225  !>@{ 3D mask arrays for diagnostics at layers (mask...L) and interfaces (mask...i)
226  real, dimension(:,:,:), pointer :: mask3dtl => null()
227  real, dimension(:,:,:), pointer :: mask3dbl => null()
228  real, dimension(:,:,:), pointer :: mask3dcul => null()
229  real, dimension(:,:,:), pointer :: mask3dcvl => null()
230  real, dimension(:,:,:), pointer :: mask3dti => null()
231  real, dimension(:,:,:), pointer :: mask3dbi => null()
232  real, dimension(:,:,:), pointer :: mask3dcui => null()
233  real, dimension(:,:,:), pointer :: mask3dcvi => null()
234  !!@}
235 end type diagcs_dsamp
236 
237 !> The following data type a list of diagnostic fields an their variants,
238 !! as well as variables that control the handling of model output.
239 type, public :: diag_ctrl
240  integer :: available_diag_doc_unit = -1 !< The unit number of a diagnostic documentation file.
241  !! This file is open if available_diag_doc_unit is > 0.
242  integer :: chksum_iounit = -1 !< The unit number of a diagnostic documentation file.
243  !! This file is open if available_diag_doc_unit is > 0.
244  logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics
245 
246 ! The following fields are used for the output of the data.
247  integer :: is !< The start i-index of cell centers within the computational domain
248  integer :: ie !< The end i-index of cell centers within the computational domain
249  integer :: js !< The start j-index of cell centers within the computational domain
250  integer :: je !< The end j-index of cell centers within the computational domain
251 
252  integer :: isd !< The start i-index of cell centers within the data domain
253  integer :: ied !< The end i-index of cell centers within the data domain
254  integer :: jsd !< The start j-index of cell centers within the data domain
255  integer :: jed !< The end j-index of cell centers within the data domain
256  real :: time_int !< The time interval for any fields
257  !! that are offered for averaging [s].
258  type(time_type) :: time_end !< The end time of the valid
259  !! interval for any offered field.
260  logical :: ave_enabled = .false. !< True if averaging is enabled.
261 
262  !>@{ The following are 3D and 2D axis groups defined for output. The names
263  !! indicate the horizontal (B, T, Cu, or Cv) and vertical (L, i, or 1) locations.
264  type(axes_grp) :: axesbl, axestl, axescul, axescvl
265  type(axes_grp) :: axesbi, axesti, axescui, axescvi
266  type(axes_grp) :: axesb1, axest1, axescu1, axescv1
267  !!@}
268  type(axes_grp) :: axeszi !< A 1-D z-space axis at interfaces
269  type(axes_grp) :: axeszl !< A 1-D z-space axis at layer centers
270  type(axes_grp) :: axesnull !< An axis group for scalars
271 
272  real, dimension(:,:), pointer :: mask2dt => null() !< 2D mask array for cell-center points
273  real, dimension(:,:), pointer :: mask2dbu => null() !< 2D mask array for cell-corner points
274  real, dimension(:,:), pointer :: mask2dcu => null() !< 2D mask array for east-face points
275  real, dimension(:,:), pointer :: mask2dcv => null() !< 2D mask array for north-face points
276  !>@{ 3D mask arrays for diagnostics at layers (mask...L) and interfaces (mask...i)
277  real, dimension(:,:,:), pointer :: mask3dtl => null()
278  real, dimension(:,:,:), pointer :: mask3dbl => null()
279  real, dimension(:,:,:), pointer :: mask3dcul => null()
280  real, dimension(:,:,:), pointer :: mask3dcvl => null()
281  real, dimension(:,:,:), pointer :: mask3dti => null()
282  real, dimension(:,:,:), pointer :: mask3dbi => null()
283  real, dimension(:,:,:), pointer :: mask3dcui => null()
284  real, dimension(:,:,:), pointer :: mask3dcvi => null()
285 
286  type(diagcs_dsamp), dimension(2:MAX_DSAMP_LEV) :: dsamp !< Downsample control container
287 
288  !!@}
289 
290 ! Space for diagnostics is dynamically allocated as it is needed.
291 ! The chunk size is how much the array should grow on each new allocation.
292 #define DIAG_ALLOC_CHUNK_SIZE 100
293  type(diag_type), dimension(:), allocatable :: diags !< The list of diagnostics
294  integer :: next_free_diag_id !< The next unused diagnostic ID
295 
296  !> default missing value to be sent to ALL diagnostics registrations
297  real :: missing_value = -1.0e+34
298 
299  !> Number of diagnostic vertical coordinates (remapped)
300  integer :: num_diag_coords
301  !> Control structure for each possible coordinate
302  type(diag_remap_ctrl), dimension(:), allocatable :: diag_remap_cs
303  type(diag_grid_storage) :: diag_grid_temp !< Stores the remapped diagnostic grid
304  logical :: diag_grid_overridden = .false. !< True if the diagnostic grids have been overriden
305 
306  type(axes_grp), dimension(:), allocatable :: &
307  remap_axeszl, & !< The 1-D z-space cell-centered axis for remapping
308  remap_axeszi !< The 1-D z-space interface axis for remapping
309  !!@{
310  type(axes_grp), dimension(:), allocatable :: remap_axestl, remap_axesbl, remap_axescul, remap_axescvl
311  type(axes_grp), dimension(:), allocatable :: remap_axesti, remap_axesbi, remap_axescui, remap_axescvi
312  !!@}
313 
314  ! Pointer to H, G and T&S needed for remapping
315  real, dimension(:,:,:), pointer :: h => null() !< The thicknesses needed for remapping
316  real, dimension(:,:,:), pointer :: t => null() !< The temperatures needed for remapping
317  real, dimension(:,:,:), pointer :: s => null() !< The salinities needed for remapping
318  type(eos_type), pointer :: eqn_of_state => null() !< The equation of state type
319  type(ocean_grid_type), pointer :: g => null() !< The ocean grid type
320  type(verticalgrid_type), pointer :: gv => null() !< The model's vertical ocean grid
321  type(unit_scale_type), pointer :: us => null() !< A dimensional unit scaling type
322 
323  !> The volume cell measure (special diagnostic) manager id
324  integer :: volume_cell_measure_dm_id = -1
325 
326 #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)
327  ! Keep a copy of h so that we know whether it has changed. If it has then
328  ! need the target grid for vertical remapping needs to have been updated.
329  real, dimension(:,:,:), allocatable :: h_old
330 #endif
331 
332  !> Number of checksum-only diagnostics
333  integer :: num_chksum_diags
334 
335 end type diag_ctrl
336 
337 ! CPU clocks
339 
340 contains
341 
342 !> Sets up diagnostics axes
343 subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical)
344  type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
345  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
346  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
347  type(param_file_type), intent(in) :: param_file !< Parameter file structure
348  type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
349  logical, optional, intent(in) :: set_vertical !< If true or missing, set up
350  !! vertical axes
351  ! Local variables
352  integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh
353  integer :: id_zl_native, id_zi_native
354  integer :: i, j, k, nz
355  real :: zlev(gv%ke), zinter(gv%ke+1)
356  logical :: set_vert
357 
358  set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical
359 
360  ! Horizontal axes for the native grids
361  if (g%symmetric) then
362  id_xq = diag_axis_init('xq', g%gridLonB(g%isgB:g%iegB), g%x_axis_units, 'x', &
363  'q point nominal longitude', domain2=g%Domain%mpp_domain)
364  id_yq = diag_axis_init('yq', g%gridLatB(g%jsgB:g%jegB), g%y_axis_units, 'y', &
365  'q point nominal latitude', domain2=g%Domain%mpp_domain)
366  else
367  id_xq = diag_axis_init('xq', g%gridLonB(g%isg:g%ieg), g%x_axis_units, 'x', &
368  'q point nominal longitude', domain2=g%Domain%mpp_domain)
369  id_yq = diag_axis_init('yq', g%gridLatB(g%jsg:g%jeg), g%y_axis_units, 'y', &
370  'q point nominal latitude', domain2=g%Domain%mpp_domain)
371  endif
372  id_xh = diag_axis_init('xh', g%gridLonT(g%isg:g%ieg), g%x_axis_units, 'x', &
373  'h point nominal longitude', domain2=g%Domain%mpp_domain)
374  id_yh = diag_axis_init('yh', g%gridLatT(g%jsg:g%jeg), g%y_axis_units, 'y', &
375  'h point nominal latitude', domain2=g%Domain%mpp_domain)
376 
377  if (set_vert) then
378  nz = gv%ke
379  zinter(1:nz+1) = gv%sInterface(1:nz+1)
380  zlev(1:nz) = gv%sLayer(1:nz)
381  id_zl = diag_axis_init('zl', zlev, trim(gv%zAxisUnits), 'z', &
382  'Layer '//trim(gv%zAxisLongName), &
383  direction=gv%direction)
384  id_zi = diag_axis_init('zi', zinter, trim(gv%zAxisUnits), 'z', &
385  'Interface '//trim(gv%zAxisLongName), &
386  direction=gv%direction)
387  else
388  id_zl = -1 ; id_zi = -1
389  endif
390  id_zl_native = id_zl ; id_zi_native = id_zi
391  ! Vertical axes for the interfaces and layers
392  call define_axes_group(diag_cs, (/ id_zi /), diag_cs%axesZi, &
393  v_cell_method='point', is_interface=.true.)
394  call define_axes_group(diag_cs, (/ id_zl /), diag_cs%axesZL, &
395  v_cell_method='mean', is_layer=.true.)
396 
397  ! Axis groupings for the model layers
398  call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zl /), diag_cs%axesTL, &
399  x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &
400  is_h_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
401  call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zl /), diag_cs%axesBL, &
402  x_cell_method='point', y_cell_method='point', v_cell_method='mean', &
403  is_q_point=.true., is_layer=.true.)
404  call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zl /), diag_cs%axesCuL, &
405  x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &
406  is_u_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
407  call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zl /), diag_cs%axesCvL, &
408  x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &
409  is_v_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
410 
411  ! Axis groupings for the model interfaces
412  call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%axesTi, &
413  x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
414  is_h_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
415  call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%axesBi, &
416  x_cell_method='point', y_cell_method='point', v_cell_method='point', &
417  is_q_point=.true., is_interface=.true.)
418  call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%axesCui, &
419  x_cell_method='point', y_cell_method='mean', v_cell_method='point', &
420  is_u_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
421  call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%axesCvi, &
422  x_cell_method='mean', y_cell_method='point', v_cell_method='point', &
423  is_v_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
424 
425  ! Axis groupings for 2-D arrays
426  call define_axes_group(diag_cs, (/ id_xh, id_yh /), diag_cs%axesT1, &
427  x_cell_method='mean', y_cell_method='mean', is_h_point=.true.)
428  call define_axes_group(diag_cs, (/ id_xq, id_yq /), diag_cs%axesB1, &
429  x_cell_method='point', y_cell_method='point', is_q_point=.true.)
430  call define_axes_group(diag_cs, (/ id_xq, id_yh /), diag_cs%axesCu1, &
431  x_cell_method='point', y_cell_method='mean', is_u_point=.true.)
432  call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, &
433  x_cell_method='mean', y_cell_method='point', is_v_point=.true.)
434 
435  ! Axis group for special null axis from diag manager
436  call define_axes_group(diag_cs, (/ null_axis_id /), diag_cs%axesNull)
437 
438 
439  !Non-native Non-downsampled
440  if (diag_cs%num_diag_coords>0) then
441  allocate(diag_cs%remap_axesZL(diag_cs%num_diag_coords))
442  allocate(diag_cs%remap_axesTL(diag_cs%num_diag_coords))
443  allocate(diag_cs%remap_axesBL(diag_cs%num_diag_coords))
444  allocate(diag_cs%remap_axesCuL(diag_cs%num_diag_coords))
445  allocate(diag_cs%remap_axesCvL(diag_cs%num_diag_coords))
446  allocate(diag_cs%remap_axesZi(diag_cs%num_diag_coords))
447  allocate(diag_cs%remap_axesTi(diag_cs%num_diag_coords))
448  allocate(diag_cs%remap_axesBi(diag_cs%num_diag_coords))
449  allocate(diag_cs%remap_axesCui(diag_cs%num_diag_coords))
450  allocate(diag_cs%remap_axesCvi(diag_cs%num_diag_coords))
451  endif
452 
453  do i=1, diag_cs%num_diag_coords
454  ! For each possible diagnostic coordinate
455  call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), gv, us, param_file)
456 
457  ! This vertical coordinate has been configured so can be used.
458  if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then
459 
460  ! This fetches the 1D-axis id for layers and interfaces and overwrite
461  ! id_zl and id_zi from above. It also returns the number of layers.
462  call diag_remap_get_axes_info(diag_cs%diag_remap_cs(i), nz, id_zl, id_zi)
463 
464  ! Axes for z layers
465  call define_axes_group(diag_cs, (/ id_zl /), diag_cs%remap_axesZL(i), &
466  nz=nz, vertical_coordinate_number=i, &
467  v_cell_method='mean', &
468  is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.)
469  call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zl /), diag_cs%remap_axesTL(i), &
470  nz=nz, vertical_coordinate_number=i, &
471  x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &
472  is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
473  xyave_axes=diag_cs%remap_axesZL(i))
474 
475  !! \note Remapping for B points is not yet implemented so needs_remapping is not
476  !! provided for remap_axesBL
477  call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zl /), diag_cs%remap_axesBL(i), &
478  nz=nz, vertical_coordinate_number=i, &
479  x_cell_method='point', y_cell_method='point', v_cell_method='mean', &
480  is_q_point=.true., is_layer=.true., is_native=.false.)
481 
482  call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zl /), diag_cs%remap_axesCuL(i), &
483  nz=nz, vertical_coordinate_number=i, &
484  x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &
485  is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
486  xyave_axes=diag_cs%remap_axesZL(i))
487 
488  call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zl /), diag_cs%remap_axesCvL(i), &
489  nz=nz, vertical_coordinate_number=i, &
490  x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &
491  is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
492  xyave_axes=diag_cs%remap_axesZL(i))
493 
494  ! Axes for z interfaces
495  call define_axes_group(diag_cs, (/ id_zi /), diag_cs%remap_axesZi(i), &
496  nz=nz, vertical_coordinate_number=i, &
497  v_cell_method='point', &
498  is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true.)
499  call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%remap_axesTi(i), &
500  nz=nz, vertical_coordinate_number=i, &
501  x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
502  is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true., &
503  xyave_axes=diag_cs%remap_axesZi(i))
504 
505  !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBi
506  call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%remap_axesBi(i), &
507  nz=nz, vertical_coordinate_number=i, &
508  x_cell_method='point', y_cell_method='point', v_cell_method='point', &
509  is_q_point=.true., is_interface=.true., is_native=.false.)
510 
511  call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%remap_axesCui(i), &
512  nz=nz, vertical_coordinate_number=i, &
513  x_cell_method='point', y_cell_method='mean', v_cell_method='point', &
514  is_u_point=.true., is_interface=.true., is_native=.false., &
515  needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))
516 
517  call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%remap_axesCvi(i), &
518  nz=nz, vertical_coordinate_number=i, &
519  x_cell_method='mean', y_cell_method='point', v_cell_method='point', &
520  is_v_point=.true., is_interface=.true., is_native=.false., &
521  needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))
522  endif
523  enddo
524 
525  !Define the downsampled axes
526  call set_axes_info_dsamp(g, gv, param_file, diag_cs, id_zl_native, id_zi_native)
527 
528  call diag_grid_storage_init(diag_cs%diag_grid_temp, g, diag_cs)
529 
530 end subroutine set_axes_info
531 
532 subroutine set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_native)
533  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
534  type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
535  type(param_file_type), intent(in) :: param_file !< Parameter file structure
536  type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
537  integer, intent(in) :: id_zl_native !< ID of native layers
538  integer, intent(in) :: id_zi_native !< ID of native interfaces
539 
540  ! Local variables
541  integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh
542  integer :: i, j, k, nz, dl
543  real, dimension(:), pointer :: gridLonT_dsamp =>null()
544  real, dimension(:), pointer :: gridLatT_dsamp =>null()
545  real, dimension(:), pointer :: gridLonB_dsamp =>null()
546  real, dimension(:), pointer :: gridLatB_dsamp =>null()
547 
548  id_zl = id_zl_native ; id_zi = id_zi_native
549  !Axes group for native downsampled diagnostics
550  do dl=2,max_dsamp_lev
551  if(dl .ne. 2) call mom_error(fatal, "set_axes_info_dsamp: Downsample level other than 2 is not supported yet!")
552  if (g%symmetric) then
553  allocate(gridlonb_dsamp(diag_cs%dsamp(dl)%isgB:diag_cs%dsamp(dl)%iegB))
554  allocate(gridlatb_dsamp(diag_cs%dsamp(dl)%jsgB:diag_cs%dsamp(dl)%jegB))
555  do i=diag_cs%dsamp(dl)%isgB,diag_cs%dsamp(dl)%iegB; gridlonb_dsamp(i) = g%gridLonB(g%isgB+dl*i); enddo
556  do j=diag_cs%dsamp(dl)%jsgB,diag_cs%dsamp(dl)%jegB; gridlatb_dsamp(j) = g%gridLatB(g%jsgB+dl*j); enddo
557  id_xq = diag_axis_init('xq', gridlonb_dsamp, g%x_axis_units, 'x', &
558  'q point nominal longitude', domain2=g%Domain%mpp_domain_d2)
559  id_yq = diag_axis_init('yq', gridlatb_dsamp, g%y_axis_units, 'y', &
560  'q point nominal latitude', domain2=g%Domain%mpp_domain_d2)
561  deallocate(gridlonb_dsamp,gridlatb_dsamp)
562  else
563  allocate(gridlonb_dsamp(diag_cs%dsamp(dl)%isg:diag_cs%dsamp(dl)%ieg))
564  allocate(gridlatb_dsamp(diag_cs%dsamp(dl)%jsg:diag_cs%dsamp(dl)%jeg))
565  do i=diag_cs%dsamp(dl)%isg,diag_cs%dsamp(dl)%ieg; gridlonb_dsamp(i) = g%gridLonB(g%isg+dl*i-2); enddo
566  do j=diag_cs%dsamp(dl)%jsg,diag_cs%dsamp(dl)%jeg; gridlatb_dsamp(j) = g%gridLatB(g%jsg+dl*j-2); enddo
567  id_xq = diag_axis_init('xq', gridlonb_dsamp, g%x_axis_units, 'x', &
568  'q point nominal longitude', domain2=g%Domain%mpp_domain_d2)
569  id_yq = diag_axis_init('yq', gridlatb_dsamp, g%y_axis_units, 'y', &
570  'q point nominal latitude', domain2=g%Domain%mpp_domain_d2)
571  deallocate(gridlonb_dsamp,gridlatb_dsamp)
572  endif
573 
574  allocate(gridlont_dsamp(diag_cs%dsamp(dl)%isg:diag_cs%dsamp(dl)%ieg))
575  allocate(gridlatt_dsamp(diag_cs%dsamp(dl)%jsg:diag_cs%dsamp(dl)%jeg))
576  do i=diag_cs%dsamp(dl)%isg,diag_cs%dsamp(dl)%ieg; gridlont_dsamp(i) = g%gridLonT(g%isg+dl*i-2); enddo
577  do j=diag_cs%dsamp(dl)%jsg,diag_cs%dsamp(dl)%jeg; gridlatt_dsamp(j) = g%gridLatT(g%jsg+dl*j-2); enddo
578  id_xh = diag_axis_init('xh', gridlont_dsamp, g%x_axis_units, 'x', &
579  'h point nominal longitude', domain2=g%Domain%mpp_domain_d2)
580  id_yh = diag_axis_init('yh', gridlatt_dsamp, g%y_axis_units, 'y', &
581  'h point nominal latitude', domain2=g%Domain%mpp_domain_d2)
582 
583  deallocate(gridlont_dsamp,gridlatt_dsamp)
584 
585  ! Axis groupings for the model layers
586  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zl /), diag_cs%dsamp(dl)%axesTL, dl, &
587  x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &
588  is_h_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
589  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zl /), diag_cs%dsamp(dl)%axesBL, dl, &
590  x_cell_method='point', y_cell_method='point', v_cell_method='mean', &
591  is_q_point=.true., is_layer=.true.)
592  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zl /), diag_cs%dsamp(dl)%axesCuL, dl, &
593  x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &
594  is_u_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
595  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zl /), diag_cs%dsamp(dl)%axesCvL, dl, &
596  x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &
597  is_v_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)
598 
599  ! Axis groupings for the model interfaces
600  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%dsamp(dl)%axesTi, dl, &
601  x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
602  is_h_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
603  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%dsamp(dl)%axesBi, dl, &
604  x_cell_method='point', y_cell_method='point', v_cell_method='point', &
605  is_q_point=.true., is_interface=.true.)
606  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%dsamp(dl)%axesCui, dl, &
607  x_cell_method='point', y_cell_method='mean', v_cell_method='point', &
608  is_u_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
609  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%dsamp(dl)%axesCvi, dl, &
610  x_cell_method='mean', y_cell_method='point', v_cell_method='point', &
611  is_v_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)
612 
613  ! Axis groupings for 2-D arrays
614  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh /), diag_cs%dsamp(dl)%axesT1, dl, &
615  x_cell_method='mean', y_cell_method='mean', is_h_point=.true.)
616  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq /), diag_cs%dsamp(dl)%axesB1, dl, &
617  x_cell_method='point', y_cell_method='point', is_q_point=.true.)
618  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh /), diag_cs%dsamp(dl)%axesCu1, dl, &
619  x_cell_method='point', y_cell_method='mean', is_u_point=.true.)
620  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq /), diag_cs%dsamp(dl)%axesCv1, dl, &
621  x_cell_method='mean', y_cell_method='point', is_v_point=.true.)
622 
623  !Non-native axes
624  if (diag_cs%num_diag_coords>0) then
625  allocate(diag_cs%dsamp(dl)%remap_axesTL(diag_cs%num_diag_coords))
626  allocate(diag_cs%dsamp(dl)%remap_axesBL(diag_cs%num_diag_coords))
627  allocate(diag_cs%dsamp(dl)%remap_axesCuL(diag_cs%num_diag_coords))
628  allocate(diag_cs%dsamp(dl)%remap_axesCvL(diag_cs%num_diag_coords))
629  allocate(diag_cs%dsamp(dl)%remap_axesTi(diag_cs%num_diag_coords))
630  allocate(diag_cs%dsamp(dl)%remap_axesBi(diag_cs%num_diag_coords))
631  allocate(diag_cs%dsamp(dl)%remap_axesCui(diag_cs%num_diag_coords))
632  allocate(diag_cs%dsamp(dl)%remap_axesCvi(diag_cs%num_diag_coords))
633  endif
634 
635  do i=1, diag_cs%num_diag_coords
636  ! For each possible diagnostic coordinate
637  !call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), GV, param_file)
638 
639  ! This vertical coordinate has been configured so can be used.
640  if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then
641 
642  ! This fetches the 1D-axis id for layers and interfaces and overwrite
643  ! id_zl and id_zi from above. It also returns the number of layers.
644  call diag_remap_get_axes_info(diag_cs%diag_remap_cs(i), nz, id_zl, id_zi)
645 
646  ! Axes for z layers
647  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zl /), diag_cs%dsamp(dl)%remap_axesTL(i), dl, &
648  nz=nz, vertical_coordinate_number=i, &
649  x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &
650  is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
651  xyave_axes=diag_cs%remap_axesZL(i))
652 
653  !! \note Remapping for B points is not yet implemented so needs_remapping is not
654  !! provided for remap_axesBL
655  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zl /), diag_cs%dsamp(dl)%remap_axesBL(i), dl, &
656  nz=nz, vertical_coordinate_number=i, &
657  x_cell_method='point', y_cell_method='point', v_cell_method='mean', &
658  is_q_point=.true., is_layer=.true., is_native=.false.)
659 
660  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zl /), diag_cs%dsamp(dl)%remap_axesCuL(i), dl, &
661  nz=nz, vertical_coordinate_number=i, &
662  x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &
663  is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
664  xyave_axes=diag_cs%remap_axesZL(i))
665 
666  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zl /), diag_cs%dsamp(dl)%remap_axesCvL(i), dl, &
667  nz=nz, vertical_coordinate_number=i, &
668  x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &
669  is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &
670  xyave_axes=diag_cs%remap_axesZL(i))
671 
672  ! Axes for z interfaces
673  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%dsamp(dl)%remap_axesTi(i), dl, &
674  nz=nz, vertical_coordinate_number=i, &
675  x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
676  is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true., &
677  xyave_axes=diag_cs%remap_axesZi(i))
678 
679  !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBi
680  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%dsamp(dl)%remap_axesBi(i), dl, &
681  nz=nz, vertical_coordinate_number=i, &
682  x_cell_method='point', y_cell_method='point', v_cell_method='point', &
683  is_q_point=.true., is_interface=.true., is_native=.false.)
684 
685  call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%dsamp(dl)%remap_axesCui(i), dl, &
686  nz=nz, vertical_coordinate_number=i, &
687  x_cell_method='point', y_cell_method='mean', v_cell_method='point', &
688  is_u_point=.true., is_interface=.true., is_native=.false., &
689  needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))
690 
691  call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%dsamp(dl)%remap_axesCvi(i), dl, &
692  nz=nz, vertical_coordinate_number=i, &
693  x_cell_method='mean', y_cell_method='point', v_cell_method='point', &
694  is_v_point=.true., is_interface=.true., is_native=.false., &
695  needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))
696  endif
697  enddo
698  enddo
699 
700 end subroutine set_axes_info_dsamp
701 
702 
703 !> set_masks_for_axes sets up the 2d and 3d masks for diagnostics using the current grid
704 !! recorded after calling diag_update_remap_grids()
705 subroutine set_masks_for_axes(G, diag_cs)
706  type(ocean_grid_type), target, intent(in) :: g !< The ocean grid type.
707  type(diag_ctrl), pointer :: diag_cs !< A pointer to a type with many variables
708  !! used for diagnostics
709  ! Local variables
710  integer :: c, nk, i, j, k, ii, jj
711  type(axes_grp), pointer :: axes => null(), h_axes => null() ! Current axes, for convenience
712 
713  do c=1, diag_cs%num_diag_coords
714  ! This vertical coordinate has been configured so can be used.
715  if (diag_remap_axes_configured(diag_cs%diag_remap_cs(c))) then
716 
717  ! Level/layer h-points in diagnostic coordinate
718  axes => diag_cs%remap_axesTL(c)
719  nk = axes%nz
720  allocate( axes%mask3d(g%isd:g%ied,g%jsd:g%jed,nk) ) ; axes%mask3d(:,:,:) = 0.
721  call diag_remap_calc_hmask(diag_cs%diag_remap_cs(c), g, axes%mask3d)
722 
723  h_axes => diag_cs%remap_axesTL(c) ! Use the h-point masks to generate the u-, v- and q- masks
724 
725  ! Level/layer u-points in diagnostic coordinate
726  axes => diag_cs%remap_axesCuL(c)
727  call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at u-layers')
728  call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
729  allocate( axes%mask3d(g%IsdB:g%IedB,g%jsd:g%jed,nk) ) ; axes%mask3d(:,:,:) = 0.
730  do k = 1, nk ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
731  if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j,k) > 0.) axes%mask3d(i,j,k) = 1.
732  enddo ; enddo ; enddo
733 
734  ! Level/layer v-points in diagnostic coordinate
735  axes => diag_cs%remap_axesCvL(c)
736  call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at v-layers')
737  call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
738  allocate( axes%mask3d(g%isd:g%ied,g%JsdB:g%JedB,nk) ) ; axes%mask3d(:,:,:) = 0.
739  do k = 1, nk ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
740  if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,j,k) = 1.
741  enddo ; enddo ; enddo
742 
743  ! Level/layer q-points in diagnostic coordinate
744  axes => diag_cs%remap_axesBL(c)
745  call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at q-layers')
746  call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
747  allocate( axes%mask3d(g%IsdB:g%IedB,g%JsdB:g%JedB,nk) ) ; axes%mask3d(:,:,:) = 0.
748  do k = 1, nk ; do j=g%jsc-1,g%jec ; do i=g%isc-1,g%iec
749  if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j+1,k) + &
750  h_axes%mask3d(i+1,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,j,k) = 1.
751  enddo ; enddo ; enddo
752 
753  ! Interface h-points in diagnostic coordinate (w-point)
754  axes => diag_cs%remap_axesTi(c)
755  call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at h-interfaces')
756  call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
757  allocate( axes%mask3d(g%isd:g%ied,g%jsd:g%jed,nk+1) ) ; axes%mask3d(:,:,:) = 0.
758  do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
759  if (h_axes%mask3d(i,j,1) > 0.) axes%mask3d(i,j,1) = 1.
760  do k = 2, nk
761  if (h_axes%mask3d(i,j,k-1) + h_axes%mask3d(i,j,k) > 0.) axes%mask3d(i,j,k) = 1.
762  enddo
763  if (h_axes%mask3d(i,j,nk) > 0.) axes%mask3d(i,j,nk+1) = 1.
764  enddo ; enddo
765 
766  h_axes => diag_cs%remap_axesTi(c) ! Use the w-point masks to generate the u-, v- and q- masks
767 
768  ! Interface u-points in diagnostic coordinate
769  axes => diag_cs%remap_axesCui(c)
770  call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at u-interfaces')
771  call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
772  allocate( axes%mask3d(g%IsdB:g%IedB,g%jsd:g%jed,nk+1) ) ; axes%mask3d(:,:,:) = 0.
773  do k = 1, nk+1 ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
774  if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j,k) > 0.) axes%mask3d(i,j,k) = 1.
775  enddo ; enddo ; enddo
776 
777  ! Interface v-points in diagnostic coordinate
778  axes => diag_cs%remap_axesCvi(c)
779  call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at v-interfaces')
780  call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
781  allocate( axes%mask3d(g%isd:g%ied,g%JsdB:g%JedB,nk+1) ) ; axes%mask3d(:,:,:) = 0.
782  do k = 1, nk+1 ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
783  if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,j,k) = 1.
784  enddo ; enddo ; enddo
785 
786  ! Interface q-points in diagnostic coordinate
787  axes => diag_cs%remap_axesBi(c)
788  call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at q-interfaces')
789  call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
790  allocate( axes%mask3d(g%IsdB:g%IedB,g%JsdB:g%JedB,nk+1) ) ; axes%mask3d(:,:,:) = 0.
791  do k = 1, nk ; do j=g%jsc-1,g%jec ; do i=g%isc-1,g%iec
792  if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j+1,k) + &
793  h_axes%mask3d(i+1,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,j,k) = 1.
794  enddo ; enddo ; enddo
795  endif
796  enddo
797 
798  !Allocate and initialize the downsampled masks for the axes
799  call set_masks_for_axes_dsamp(g, diag_cs)
800 
801 end subroutine set_masks_for_axes
802 
803 subroutine set_masks_for_axes_dsamp(G, diag_cs)
804  type(ocean_grid_type), target, intent(in) :: G !< The ocean grid type.
805  type(diag_ctrl), pointer :: diag_cs !< A pointer to a type with many variables
806  !! used for diagnostics
807  ! Local variables
808  integer :: c, nk, i, j, k, ii, jj
809  integer :: dl
810  type(axes_grp), pointer :: axes => null(), h_axes => null() ! Current axes, for convenience
811 
812  !Each downsampled axis needs both downsampled and non-downsampled mask
813  !The downsampled mask is needed for sending out the diagnostics output via diag_manager
814  !The non-downsampled mask is needed for downsampling the diagnostics field
815  do dl=2,max_dsamp_lev
816  if(dl .ne. 2) call mom_error(fatal, "set_masks_for_axes_dsamp: Downsample level other than 2 is not supported!")
817  do c=1, diag_cs%num_diag_coords
818  ! Level/layer h-points in diagnostic coordinate
819  axes => diag_cs%remap_axesTL(c)
820  call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTL(c)%dsamp(dl)%mask3d, dl,g%isc, g%jsc, &
821  g%HId2%isc, g%HId2%iec, g%HId2%jsc, g%HId2%jec, g%HId2%isd, g%HId2%ied, g%HId2%jsd, g%HId2%jed)
822  diag_cs%dsamp(dl)%remap_axesTL(c)%mask3d => axes%mask3d !set non-downsampled mask
823  ! Level/layer u-points in diagnostic coordinate
824  axes => diag_cs%remap_axesCuL(c)
825  call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCuL(c)%dsamp(dl)%mask3d, dl,g%IscB,g%JscB, &
826  g%HId2%IscB,g%HId2%IecB,g%HId2%jsc, g%HId2%jec,g%HId2%IsdB,g%HId2%IedB,g%HId2%jsd, g%HId2%jed)
827  diag_cs%dsamp(dl)%remap_axesCul(c)%mask3d => axes%mask3d !set non-downsampled mask
828  ! Level/layer v-points in diagnostic coordinate
829  axes => diag_cs%remap_axesCvL(c)
830  call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvL(c)%dsamp(dl)%mask3d, dl,g%isc ,g%JscB, &
831  g%HId2%isc ,g%HId2%iec, g%HId2%JscB,g%HId2%JecB,g%HId2%isd ,g%HId2%ied, g%HId2%JsdB,g%HId2%JedB)
832  diag_cs%dsamp(dl)%remap_axesCvL(c)%mask3d => axes%mask3d !set non-downsampled mask
833  ! Level/layer q-points in diagnostic coordinate
834  axes => diag_cs%remap_axesBL(c)
835  call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBL(c)%dsamp(dl)%mask3d, dl,g%IscB,g%JscB, &
836  g%HId2%IscB,g%HId2%IecB,g%HId2%JscB,g%HId2%JecB,g%HId2%IsdB,g%HId2%IedB,g%HId2%JsdB,g%HId2%JedB)
837  diag_cs%dsamp(dl)%remap_axesBL(c)%mask3d => axes%mask3d !set non-downsampled mask
838  ! Interface h-points in diagnostic coordinate (w-point)
839  axes => diag_cs%remap_axesTi(c)
840  call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTi(c)%dsamp(dl)%mask3d, dl,g%isc, g%jsc, &
841  g%HId2%isc, g%HId2%iec, g%HId2%jsc, g%HId2%jec, g%HId2%isd, g%HId2%ied, g%HId2%jsd, g%HId2%jed)
842  diag_cs%dsamp(dl)%remap_axesTi(c)%mask3d => axes%mask3d !set non-downsampled mask
843  ! Interface u-points in diagnostic coordinate
844  axes => diag_cs%remap_axesCui(c)
845  call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCui(c)%dsamp(dl)%mask3d, dl,g%IscB,g%JscB, &
846  g%HId2%IscB,g%HId2%IecB,g%HId2%jsc, g%HId2%jec,g%HId2%IsdB,g%HId2%IedB,g%HId2%jsd, g%HId2%jed)
847  diag_cs%dsamp(dl)%remap_axesCui(c)%mask3d => axes%mask3d !set non-downsampled mask
848  ! Interface v-points in diagnostic coordinate
849  axes => diag_cs%remap_axesCvi(c)
850  call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvi(c)%dsamp(dl)%mask3d, dl,g%isc ,g%JscB, &
851  g%HId2%isc ,g%HId2%iec, g%HId2%JscB,g%HId2%JecB,g%HId2%isd ,g%HId2%ied, g%HId2%JsdB,g%HId2%JedB)
852  diag_cs%dsamp(dl)%remap_axesCvi(c)%mask3d => axes%mask3d !set non-downsampled mask
853  ! Interface q-points in diagnostic coordinate
854  axes => diag_cs%remap_axesBi(c)
855  call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBi(c)%dsamp(dl)%mask3d, dl,g%IscB,g%JscB, &
856  g%HId2%IscB,g%HId2%IecB,g%HId2%JscB,g%HId2%JecB,g%HId2%IsdB,g%HId2%IedB,g%HId2%JsdB,g%HId2%JedB)
857  diag_cs%dsamp(dl)%remap_axesBi(c)%mask3d => axes%mask3d !set non-downsampled mask
858  enddo
859  enddo
860 end subroutine set_masks_for_axes_dsamp
861 
862 !> Attaches the id of cell areas to axes groups for use with cell_measures
863 subroutine diag_register_area_ids(diag_cs, id_area_t, id_area_q)
864  type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
865  integer, optional, intent(in) :: id_area_t !< Diag_mediator id for area of h-cells
866  integer, optional, intent(in) :: id_area_q !< Diag_mediator id for area of q-cells
867  ! Local variables
868  integer :: fms_id, i
869  if (present(id_area_t)) then
870  fms_id = diag_cs%diags(id_area_t)%fms_diag_id
871  diag_cs%axesT1%id_area = fms_id
872  diag_cs%axesTi%id_area = fms_id
873  diag_cs%axesTL%id_area = fms_id
874  do i=1, diag_cs%num_diag_coords
875  diag_cs%remap_axesTL(i)%id_area = fms_id
876  diag_cs%remap_axesTi(i)%id_area = fms_id
877  enddo
878  endif
879  if (present(id_area_q)) then
880  fms_id = diag_cs%diags(id_area_q)%fms_diag_id
881  diag_cs%axesB1%id_area = fms_id
882  diag_cs%axesBi%id_area = fms_id
883  diag_cs%axesBL%id_area = fms_id
884  do i=1, diag_cs%num_diag_coords
885  diag_cs%remap_axesBL(i)%id_area = fms_id
886  diag_cs%remap_axesBi(i)%id_area = fms_id
887  enddo
888  endif
889 end subroutine diag_register_area_ids
890 
891 !> Sets a handle inside diagnostics mediator to associate 3d cell measures
892 subroutine register_cell_measure(G, diag, Time)
893  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
894  type(diag_ctrl), target, intent(inout) :: diag !< Regulates diagnostic output
895  type(time_type), intent(in) :: time !< Model time
896  ! Local variables
897  integer :: id
898  id = register_diag_field('ocean_model', 'volcello', diag%axesTL, &
899  time, 'Ocean grid-cell volume', 'm3', &
900  standard_name='ocean_volume', v_extensive=.true., &
901  x_cell_method='sum', y_cell_method='sum')
903 
904 end subroutine register_cell_measure
905 
906 !> Attaches the id of cell volumes to axes groups for use with cell_measures
907 subroutine diag_associate_volume_cell_measure(diag_cs, id_h_volume)
908  type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
909  integer, intent(in) :: id_h_volume !< Diag_manager id for volume of h-cells
910  ! Local variables
911  type(diag_type), pointer :: tmp => null()
912 
913  if (id_h_volume<=0) return ! Do nothing
914  diag_cs%volume_cell_measure_dm_id = id_h_volume ! Record for diag_get_volume_cell_measure_dm_id()
915 
916  ! Set the cell measure for this axes group to the FMS id in this coordinate system
917  diag_cs%diags(id_h_volume)%axes%id_volume = diag_cs%diags(id_h_volume)%fms_diag_id
918 
919  tmp => diag_cs%diags(id_h_volume)%next ! First item in the list, if any
920  do while (associated(tmp))
921  ! Set the cell measure for this axes group to the FMS id in this coordinate system
922  tmp%axes%id_volume = tmp%fms_diag_id
923  tmp => tmp%next ! Move to next axes group for this field
924  enddo
925 
927 
928 !> Returns diag_manager id for cell measure of h-cells
929 integer function diag_get_volume_cell_measure_dm_id(diag_cs)
930  type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics control structure
931 
932  diag_get_volume_cell_measure_dm_id = diag_cs%volume_cell_measure_dm_id
933 
935 
936 !> Defines a group of "axes" from list of handles
937 subroutine define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_number, &
938  x_cell_method, y_cell_method, v_cell_method, &
939  is_h_point, is_q_point, is_u_point, is_v_point, &
940  is_layer, is_interface, &
941  is_native, needs_remapping, needs_interpolating, &
942  xyave_axes)
943  type(diag_ctrl), target, intent(in) :: diag_cs !< Diagnostics control structure
944  integer, dimension(:), intent(in) :: handles !< A list of 1D axis handles
945  type(axes_grp), intent(out) :: axes !< The group of 1D axes
946  integer, optional, intent(in) :: nz !< Number of layers in this diagnostic grid
947  integer, optional, intent(in) :: vertical_coordinate_number !< Index number for vertical coordinate
948  character(len=*), optional, intent(in) :: x_cell_method !< A x-direction cell method used to construct the
949  !! "cell_methods" attribute in CF convention
950  character(len=*), optional, intent(in) :: y_cell_method !< A y-direction cell method used to construct the
951  !! "cell_methods" attribute in CF convention
952  character(len=*), optional, intent(in) :: v_cell_method !< A vertical direction cell method used to construct
953  !! the "cell_methods" attribute in CF convention
954  logical, optional, intent(in) :: is_h_point !< If true, indicates this axes group for h-point
955  !! located fields
956  logical, optional, intent(in) :: is_q_point !< If true, indicates this axes group for q-point
957  !! located fields
958  logical, optional, intent(in) :: is_u_point !< If true, indicates this axes group for
959  !! u-point located fields
960  logical, optional, intent(in) :: is_v_point !< If true, indicates this axes group for
961  !! v-point located fields
962  logical, optional, intent(in) :: is_layer !< If true, indicates that this axes group is
963  !! for a layer vertically-located field.
964  logical, optional, intent(in) :: is_interface !< If true, indicates that this axes group
965  !! is for an interface vertically-located field.
966  logical, optional, intent(in) :: is_native !< If true, indicates that this axes group is
967  !! for a native model grid. False for any other grid.
968  logical, optional, intent(in) :: needs_remapping !< If true, indicates that this axes group is
969  !! for a intensive layer-located field that must
970  !! be remapped to these axes. Used for rank>2.
971  logical, optional, intent(in) :: needs_interpolating !< If true, indicates that this axes group
972  !! is for a sampled interface-located field that must
973  !! be interpolated to these axes. Used for rank>2.
974  type(axes_grp), optional, target :: xyave_axes !< The corresponding axes group for horizontally
975  !! area-average diagnostics
976  ! Local variables
977  integer :: n
978 
979  n = size(handles)
980  if (n<1 .or. n>3) call mom_error(fatal, "define_axes_group: wrong size for list of handles!")
981  allocate( axes%handles(n) )
982  axes%id = i2s(handles, n) ! Identifying string
983  axes%rank = n
984  axes%handles(:) = handles(:)
985  axes%diag_cs => diag_cs ! A [circular] link back to the diag_cs structure
986  if (present(x_cell_method)) then
987  if (axes%rank<2) call mom_error(fatal, 'define_axes_group: ' // &
988  'Can not set x_cell_method for rank<2.')
989  axes%x_cell_method = trim(x_cell_method)
990  else
991  axes%x_cell_method = ''
992  endif
993  if (present(y_cell_method)) then
994  if (axes%rank<2) call mom_error(fatal, 'define_axes_group: ' // &
995  'Can not set y_cell_method for rank<2.')
996  axes%y_cell_method = trim(y_cell_method)
997  else
998  axes%y_cell_method = ''
999  endif
1000  if (present(v_cell_method)) then
1001  if (axes%rank/=1 .and. axes%rank/=3) call mom_error(fatal, 'define_axes_group: ' // &
1002  'Can not set v_cell_method for rank<>1 or 3.')
1003  axes%v_cell_method = trim(v_cell_method)
1004  else
1005  axes%v_cell_method = ''
1006  endif
1007  if (present(nz)) axes%nz = nz
1008  if (present(vertical_coordinate_number)) axes%vertical_coordinate_number = vertical_coordinate_number
1009  if (present(is_h_point)) axes%is_h_point = is_h_point
1010  if (present(is_q_point)) axes%is_q_point = is_q_point
1011  if (present(is_u_point)) axes%is_u_point = is_u_point
1012  if (present(is_v_point)) axes%is_v_point = is_v_point
1013  if (present(is_layer)) axes%is_layer = is_layer
1014  if (present(is_interface)) axes%is_interface = is_interface
1015  if (present(is_native)) axes%is_native = is_native
1016  if (present(needs_remapping)) axes%needs_remapping = needs_remapping
1017  if (present(needs_interpolating)) axes%needs_interpolating = needs_interpolating
1018  if (present(xyave_axes)) axes%xyave_axes => xyave_axes
1019 
1020  ! Setup masks for this axes group
1021  axes%mask2d => null()
1022  if (axes%rank==2) then
1023  if (axes%is_h_point) axes%mask2d => diag_cs%mask2dT
1024  if (axes%is_u_point) axes%mask2d => diag_cs%mask2dCu
1025  if (axes%is_v_point) axes%mask2d => diag_cs%mask2dCv
1026  if (axes%is_q_point) axes%mask2d => diag_cs%mask2dBu
1027  endif
1028  ! A static 3d mask for non-native coordinates can only be setup when a grid is available
1029  axes%mask3d => null()
1030  if (axes%rank==3 .and. axes%is_native) then
1031  ! Native variables can/should use the native masks copied into diag_cs
1032  if (axes%is_layer) then
1033  if (axes%is_h_point) axes%mask3d => diag_cs%mask3dTL
1034  if (axes%is_u_point) axes%mask3d => diag_cs%mask3dCuL
1035  if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvL
1036  if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBL
1037  elseif (axes%is_interface) then
1038  if (axes%is_h_point) axes%mask3d => diag_cs%mask3dTi
1039  if (axes%is_u_point) axes%mask3d => diag_cs%mask3dCui
1040  if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvi
1041  if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBi
1042  endif
1043  endif
1044 
1045 end subroutine define_axes_group
1046 
1047 !> Defines a group of downsampled "axes" from list of handles
1048 subroutine define_axes_group_dsamp(diag_cs, handles, axes, dl, nz, vertical_coordinate_number, &
1049  x_cell_method, y_cell_method, v_cell_method, &
1050  is_h_point, is_q_point, is_u_point, is_v_point, &
1051  is_layer, is_interface, &
1052  is_native, needs_remapping, needs_interpolating, &
1053  xyave_axes)
1054  type(diag_ctrl), target, intent(in) :: diag_cs !< Diagnostics control structure
1055  integer, dimension(:), intent(in) :: handles !< A list of 1D axis handles
1056  type(axes_grp), intent(out) :: axes !< The group of 1D axes
1057  integer, intent(in) :: dl !< Downsample level
1058  integer, optional, intent(in) :: nz !< Number of layers in this diagnostic grid
1059  integer, optional, intent(in) :: vertical_coordinate_number !< Index number for vertical coordinate
1060  character(len=*), optional, intent(in) :: x_cell_method !< A x-direction cell method used to construct the
1061  !! "cell_methods" attribute in CF convention
1062  character(len=*), optional, intent(in) :: y_cell_method !< A y-direction cell method used to construct the
1063  !! "cell_methods" attribute in CF convention
1064  character(len=*), optional, intent(in) :: v_cell_method !< A vertical direction cell method used to construct
1065  !! the "cell_methods" attribute in CF convention
1066  logical, optional, intent(in) :: is_h_point !< If true, indicates this axes group for h-point
1067  !! located fields
1068  logical, optional, intent(in) :: is_q_point !< If true, indicates this axes group for q-point
1069  !! located fields
1070  logical, optional, intent(in) :: is_u_point !< If true, indicates this axes group for
1071  !! u-point located fields
1072  logical, optional, intent(in) :: is_v_point !< If true, indicates this axes group for
1073  !! v-point located fields
1074  logical, optional, intent(in) :: is_layer !< If true, indicates that this axes group is
1075  !! for a layer vertically-located field.
1076  logical, optional, intent(in) :: is_interface !< If true, indicates that this axes group
1077  !! is for an interface vertically-located field.
1078  logical, optional, intent(in) :: is_native !< If true, indicates that this axes group is
1079  !! for a native model grid. False for any other grid.
1080  logical, optional, intent(in) :: needs_remapping !< If true, indicates that this axes group is
1081  !! for a intensive layer-located field that must
1082  !! be remapped to these axes. Used for rank>2.
1083  logical, optional, intent(in) :: needs_interpolating !< If true, indicates that this axes group
1084  !! is for a sampled interface-located field that must
1085  !! be interpolated to these axes. Used for rank>2.
1086  type(axes_grp), optional, target :: xyave_axes !< The corresponding axes group for horizontally
1087  !! area-average diagnostics
1088  ! Local variables
1089  integer :: n
1090 
1091  n = size(handles)
1092  if (n<1 .or. n>3) call mom_error(fatal, "define_axes_group: wrong size for list of handles!")
1093  allocate( axes%handles(n) )
1094  axes%id = i2s(handles, n) ! Identifying string
1095  axes%rank = n
1096  axes%handles(:) = handles(:)
1097  axes%diag_cs => diag_cs ! A [circular] link back to the diag_cs structure
1098  if (present(x_cell_method)) then
1099  if (axes%rank<2) call mom_error(fatal, 'define_axes_group: ' // &
1100  'Can not set x_cell_method for rank<2.')
1101  axes%x_cell_method = trim(x_cell_method)
1102  else
1103  axes%x_cell_method = ''
1104  endif
1105  if (present(y_cell_method)) then
1106  if (axes%rank<2) call mom_error(fatal, 'define_axes_group: ' // &
1107  'Can not set y_cell_method for rank<2.')
1108  axes%y_cell_method = trim(y_cell_method)
1109  else
1110  axes%y_cell_method = ''
1111  endif
1112  if (present(v_cell_method)) then
1113  if (axes%rank/=1 .and. axes%rank/=3) call mom_error(fatal, 'define_axes_group: ' // &
1114  'Can not set v_cell_method for rank<>1 or 3.')
1115  axes%v_cell_method = trim(v_cell_method)
1116  else
1117  axes%v_cell_method = ''
1118  endif
1119  axes%downsample_level = dl
1120  if (present(nz)) axes%nz = nz
1121  if (present(vertical_coordinate_number)) axes%vertical_coordinate_number = vertical_coordinate_number
1122  if (present(is_h_point)) axes%is_h_point = is_h_point
1123  if (present(is_q_point)) axes%is_q_point = is_q_point
1124  if (present(is_u_point)) axes%is_u_point = is_u_point
1125  if (present(is_v_point)) axes%is_v_point = is_v_point
1126  if (present(is_layer)) axes%is_layer = is_layer
1127  if (present(is_interface)) axes%is_interface = is_interface
1128  if (present(is_native)) axes%is_native = is_native
1129  if (present(needs_remapping)) axes%needs_remapping = needs_remapping
1130  if (present(needs_interpolating)) axes%needs_interpolating = needs_interpolating
1131  if (present(xyave_axes)) axes%xyave_axes => xyave_axes
1132 
1133  ! Setup masks for this axes group
1134 
1135  axes%mask2d => null()
1136  if (axes%rank==2) then
1137  if (axes%is_h_point) axes%mask2d => diag_cs%mask2dT
1138  if (axes%is_u_point) axes%mask2d => diag_cs%mask2dCu
1139  if (axes%is_v_point) axes%mask2d => diag_cs%mask2dCv
1140  if (axes%is_q_point) axes%mask2d => diag_cs%mask2dBu
1141  endif
1142  ! A static 3d mask for non-native coordinates can only be setup when a grid is available
1143  axes%mask3d => null()
1144  if (axes%rank==3 .and. axes%is_native) then
1145  ! Native variables can/should use the native masks copied into diag_cs
1146  if (axes%is_layer) then
1147  if (axes%is_h_point) axes%mask3d => diag_cs%mask3dTL
1148  if (axes%is_u_point) axes%mask3d => diag_cs%mask3dCuL
1149  if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvL
1150  if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBL
1151  elseif (axes%is_interface) then
1152  if (axes%is_h_point) axes%mask3d => diag_cs%mask3dTi
1153  if (axes%is_u_point) axes%mask3d => diag_cs%mask3dCui
1154  if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvi
1155  if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBi
1156  endif
1157  endif
1158 
1159  axes%dsamp(dl)%mask2d => null()
1160  if (axes%rank==2) then
1161  if (axes%is_h_point) axes%dsamp(dl)%mask2d => diag_cs%dsamp(dl)%mask2dT
1162  if (axes%is_u_point) axes%dsamp(dl)%mask2d => diag_cs%dsamp(dl)%mask2dCu
1163  if (axes%is_v_point) axes%dsamp(dl)%mask2d => diag_cs%dsamp(dl)%mask2dCv
1164  if (axes%is_q_point) axes%dsamp(dl)%mask2d => diag_cs%dsamp(dl)%mask2dBu
1165  endif
1166  ! A static 3d mask for non-native coordinates can only be setup when a grid is available
1167  axes%dsamp(dl)%mask3d => null()
1168  if (axes%rank==3 .and. axes%is_native) then
1169  ! Native variables can/should use the native masks copied into diag_cs
1170  if (axes%is_layer) then
1171  if (axes%is_h_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dTL
1172  if (axes%is_u_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dCuL
1173  if (axes%is_v_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dCvL
1174  if (axes%is_q_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dBL
1175  elseif (axes%is_interface) then
1176  if (axes%is_h_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dTi
1177  if (axes%is_u_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dCui
1178  if (axes%is_v_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dCvi
1179  if (axes%is_q_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dBi
1180  endif
1181  endif
1182 
1183 end subroutine define_axes_group_dsamp
1184 
1185 !> Set up the array extents for doing diagnostics
1186 subroutine set_diag_mediator_grid(G, diag_cs)
1187  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
1188  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output
1189 
1190  diag_cs%is = g%isc - (g%isd-1) ; diag_cs%ie = g%iec - (g%isd-1)
1191  diag_cs%js = g%jsc - (g%jsd-1) ; diag_cs%je = g%jec - (g%jsd-1)
1192  diag_cs%isd = g%isd ; diag_cs%ied = g%ied
1193  diag_cs%jsd = g%jsd ; diag_cs%jed = g%jed
1194 
1195 end subroutine set_diag_mediator_grid
1196 
1197 !> Make a real scalar diagnostic available for averaging or output
1198 subroutine post_data_0d(diag_field_id, field, diag_cs, is_static)
1199  integer, intent(in) :: diag_field_id !< The id for an output variable returned by a
1200  !! previous call to register_diag_field.
1201  real, intent(in) :: field !< real value being offered for output or averaging
1202  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output
1203  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.
1204 
1205  ! Local variables
1206  logical :: used, is_stat
1207  type(diag_type), pointer :: diag => null()
1208 
1209  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)
1210  is_stat = .false. ; if (present(is_static)) is_stat = is_static
1211 
1212  ! Iterate over list of diag 'variants', e.g. CMOR aliases, call send_data
1213  ! for each one.
1214  call assert(diag_field_id < diag_cs%next_free_diag_id, &
1215  'post_data_0d: Unregistered diagnostic id')
1216  diag => diag_cs%diags(diag_field_id)
1217  do while (associated(diag))
1218  if (diag_cs%diag_as_chksum) then
1219  call chksum0(field, diag%debug_str, logunit=diag_cs%chksum_iounit)
1220  else if (is_stat) then
1221  used = send_data(diag%fms_diag_id, field)
1222  elseif (diag_cs%ave_enabled) then
1223  used = send_data(diag%fms_diag_id, field, diag_cs%time_end)
1224  endif
1225  diag => diag%next
1226  enddo
1227 
1228  if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator)
1229 end subroutine post_data_0d
1230 
1231 !> Make a real 1-d array diagnostic available for averaging or output
1232 subroutine post_data_1d_k(diag_field_id, field, diag_cs, is_static)
1233  integer, intent(in) :: diag_field_id !< The id for an output variable returned by a
1234  !! previous call to register_diag_field.
1235  real, target, intent(in) :: field(:) !< 1-d array being offered for output or averaging
1236  type(diag_ctrl), target, intent(in) :: diag_cs !< Structure used to regulate diagnostic output
1237  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.
1238 
1239  ! Local variables
1240  logical :: used ! The return value of send_data is not used for anything.
1241  real, dimension(:), pointer :: locfield => null()
1242  logical :: is_stat
1243  integer :: k, ks, ke
1244  type(diag_type), pointer :: diag => null()
1245 
1246  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)
1247  is_stat = .false. ; if (present(is_static)) is_stat = is_static
1248 
1249  ! Iterate over list of diag 'variants', e.g. CMOR aliases.
1250  call assert(diag_field_id < diag_cs%next_free_diag_id, &
1251  'post_data_1d_k: Unregistered diagnostic id')
1252  diag => diag_cs%diags(diag_field_id)
1253  do while (associated(diag))
1254 
1255  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) then
1256  ks = lbound(field,1) ; ke = ubound(field,1)
1257  allocate( locfield( ks:ke ) )
1258 
1259  do k=ks,ke
1260  if (field(k) == diag_cs%missing_value) then
1261  locfield(k) = diag_cs%missing_value
1262  else
1263  locfield(k) = field(k) * diag%conversion_factor
1264  endif
1265  enddo
1266  else
1267  locfield => field
1268  endif
1269 
1270  if (diag_cs%diag_as_chksum) then
1271  call zchksum(locfield, diag%debug_str, logunit=diag_cs%chksum_iounit)
1272  else if (is_stat) then
1273  used = send_data(diag%fms_diag_id, locfield)
1274  elseif (diag_cs%ave_enabled) then
1275  used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, weight=diag_cs%time_int)
1276  endif
1277  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield )
1278 
1279  diag => diag%next
1280  enddo
1281 
1282  if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator)
1283 end subroutine post_data_1d_k
1284 
1285 !> Make a real 2-d array diagnostic available for averaging or output
1286 subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask)
1287  integer, intent(in) :: diag_field_id !< The id for an output variable returned by a
1288  !! previous call to register_diag_field.
1289  real, intent(in) :: field(:,:) !< 2-d array being offered for output or averaging
1290  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output
1291  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.
1292  real, optional, intent(in) :: mask(:,:) !< If present, use this real array as the data mask.
1293 
1294  ! Local variables
1295  type(diag_type), pointer :: diag => null()
1296 
1297  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)
1298 
1299  ! Iterate over list of diag 'variants' (e.g. CMOR aliases) and post each.
1300  call assert(diag_field_id < diag_cs%next_free_diag_id, &
1301  'post_data_2d: Unregistered diagnostic id')
1302  diag => diag_cs%diags(diag_field_id)
1303  do while (associated(diag))
1304  call post_data_2d_low(diag, field, diag_cs, is_static, mask)
1305  diag => diag%next
1306  enddo
1307 
1308  if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator)
1309 end subroutine post_data_2d
1310 
1311 !> Make a real 2-d array diagnostic available for averaging or output
1312 !! using a diag_type instead of an integer id.
1313 subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask)
1314  type(diag_type), intent(in) :: diag !< A structure describing the diagnostic to post
1315  real, target, intent(in) :: field(:,:) !< 2-d array being offered for output or averaging
1316  type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output
1317  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.
1318  real, optional,target, intent(in) :: mask(:,:) !< If present, use this real array as the data mask.
1319 
1320  ! Local variables
1321  real, dimension(:,:), pointer :: locfield
1322  real, dimension(:,:), pointer :: locmask
1323  character(len=300) :: mesg
1324  logical :: used, is_stat
1325  integer :: cszi, cszj, dszi, dszj
1326  integer :: isv, iev, jsv, jev, i, j, chksum, isv_o,jsv_o
1327  real, dimension(:,:), allocatable, target :: locfield_dsamp
1328  real, dimension(:,:), allocatable, target :: locmask_dsamp
1329  integer :: dl
1330 
1331  locfield => null()
1332  locmask => null()
1333  is_stat = .false. ; if (present(is_static)) is_stat = is_static
1334 
1335  ! Determine the propery array indices, noting that because of the (:,:)
1336  ! declaration of field, symmetric arrays are using a SW-grid indexing,
1337  ! but non-symmetric arrays are using a NE-grid indexing. Send_data
1338  ! actually only uses the difference between ie and is to determine
1339  ! the output data size and assumes that halos are symmetric.
1340  isv = diag_cs%is ; iev = diag_cs%ie ; jsv = diag_cs%js ; jev = diag_cs%je
1341 
1342  cszi = diag_cs%ie-diag_cs%is +1 ; dszi = diag_cs%ied-diag_cs%isd +1
1343  cszj = diag_cs%je-diag_cs%js +1 ; dszj = diag_cs%jed-diag_cs%jsd +1
1344  if ( size(field,1) == dszi ) then
1345  isv = diag_cs%is ; iev = diag_cs%ie ! Data domain
1346  elseif ( size(field,1) == dszi + 1 ) then
1347  isv = diag_cs%is ; iev = diag_cs%ie+1 ! Symmetric data domain
1348  elseif ( size(field,1) == cszi) then
1349  isv = 1 ; iev = cszi ! Computational domain
1350  elseif ( size(field,1) == cszi + 1 ) then
1351  isv = 1 ; iev = cszi+1 ! Symmetric computational domain
1352  else
1353  write (mesg,*) " peculiar size ",size(field,1)," in i-direction\n"//&
1354  "does not match one of ", cszi, cszi+1, dszi, dszi+1
1355  call mom_error(fatal,"post_data_2d_low: "//trim(diag%debug_str)//trim(mesg))
1356  endif
1357 
1358  if ( size(field,2) == dszj ) then
1359  jsv = diag_cs%js ; jev = diag_cs%je ! Data domain
1360  elseif ( size(field,2) == dszj + 1 ) then
1361  jsv = diag_cs%js ; jev = diag_cs%je+1 ! Symmetric data domain
1362  elseif ( size(field,2) == cszj ) then
1363  jsv = 1 ; jev = cszj ! Computational domain
1364  elseif ( size(field,2) == cszj+1 ) then
1365  jsv = 1 ; jev = cszj+1 ! Symmetric computational domain
1366  else
1367  write (mesg,*) " peculiar size ",size(field,2)," in j-direction\n"//&
1368  "does not match one of ", cszj, cszj+1, dszj, dszj+1
1369  call mom_error(fatal,"post_data_2d_low: "//trim(diag%debug_str)//trim(mesg))
1370  endif
1371 
1372  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) then
1373  allocate( locfield( lbound(field,1):ubound(field,1), lbound(field,2):ubound(field,2) ) )
1374  do j=jsv,jev ; do i=isv,iev
1375  if (field(i,j) == diag_cs%missing_value) then
1376  locfield(i,j) = diag_cs%missing_value
1377  else
1378  locfield(i,j) = field(i,j) * diag%conversion_factor
1379  endif
1380  enddo ; enddo
1381  locfield(isv:iev,jsv:jev) = field(isv:iev,jsv:jev) * diag%conversion_factor
1382  else
1383  locfield => field
1384  endif
1385 
1386  if (present(mask)) then
1387  locmask => mask
1388  elseif(.NOT. is_stat) then
1389  if(associated(diag%axes%mask2d)) locmask => diag%axes%mask2d
1390  endif
1391 
1392  dl=1
1393  if(.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet
1394  !Downsample the diag field and mask (if present)
1395  if (dl > 1) then
1396  isv_o=isv ; jsv_o=jsv
1397  call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask)
1398  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield )
1399  locfield => locfield_dsamp
1400  if (present(mask)) then
1401  call downsample_field_2d(locmask, locmask_dsamp, dl, msk, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev)
1402  locmask => locmask_dsamp
1403  elseif(associated(diag%axes%dsamp(dl)%mask2d)) then
1404  locmask => diag%axes%dsamp(dl)%mask2d
1405  endif
1406  endif
1407 
1408  if (diag_cs%diag_as_chksum) then
1409  if (diag%axes%is_h_point) then
1410  call hchksum(locfield, diag%debug_str, diag_cs%G%HI, &
1411  logunit=diag_cs%chksum_iounit)
1412  else if (diag%axes%is_u_point) then
1413  call uchksum(locfield, diag%debug_str, diag_cs%G%HI, &
1414  logunit=diag_cs%chksum_iounit)
1415  else if (diag%axes%is_v_point) then
1416  call vchksum(locfield, diag%debug_str, diag_cs%G%HI, &
1417  logunit=diag_cs%chksum_iounit)
1418  else if (diag%axes%is_q_point) then
1419  call bchksum(locfield, diag%debug_str, diag_cs%G%HI, &
1420  logunit=diag_cs%chksum_iounit)
1421  else
1422  call mom_error(fatal, "post_data_2d_low: unknown axis type.")
1423  endif
1424  else
1425  if (is_stat) then
1426  if (present(mask)) then
1427  call assert(size(locfield) == size(locmask), &
1428  'post_data_2d_low is_stat: mask size mismatch: '//diag%debug_str)
1429  used = send_data(diag%fms_diag_id, locfield, &
1430  is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=locmask)
1431  !elseif (associated(diag%axes%mask2d)) then
1432  ! used = send_data(diag%fms_diag_id, locfield, &
1433  ! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%axes%mask2d)
1434  else
1435  used = send_data(diag%fms_diag_id, locfield, &
1436  is_in=isv, js_in=jsv, ie_in=iev, je_in=jev)
1437  endif
1438  elseif (diag_cs%ave_enabled) then
1439  if (associated(locmask)) then
1440  call assert(size(locfield) == size(locmask), &
1441  'post_data_2d_low: mask size mismatch: '//diag%debug_str)
1442  used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, &
1443  is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, &
1444  weight=diag_cs%time_int, rmask=locmask)
1445  else
1446  used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, &
1447  is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, &
1448  weight=diag_cs%time_int)
1449  endif
1450  endif
1451  endif
1452  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.) .and. dl<2) &
1453  deallocate( locfield )
1454 end subroutine post_data_2d_low
1455 
1456 !> Make a real 3-d array diagnostic available for averaging or output.
1457 subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h)
1459  integer, intent(in) :: diag_field_id !< The id for an output variable returned by a
1460  !! previous call to register_diag_field.
1461  real, intent(in) :: field(:,:,:) !< 3-d array being offered for output or averaging
1462  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output
1463  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.
1464  real, optional, intent(in) :: mask(:,:,:) !< If present, use this real array as the data mask.
1465  real, dimension(:,:,:), &
1466  target, optional, intent(in) :: alt_h !< An alternate thickness to use for vertically
1467  !! remapping this diagnostic [H ~> m or kg m-2].
1468 
1469  ! Local variables
1470  type(diag_type), pointer :: diag => null()
1471  integer :: nz, i, j, k
1472  real, dimension(:,:,:), allocatable :: remapped_field
1473  logical :: staggered_in_x, staggered_in_y
1474  real, dimension(:,:,:), pointer :: h_diag => null()
1475 
1476  if (present(alt_h)) then
1477  h_diag => alt_h
1478  else
1479  h_diag => diag_cs%h
1480  endif
1481 
1482  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)
1483 
1484  ! Iterate over list of diag 'variants', e.g. CMOR aliases, different vertical
1485  ! grids, and post each.
1486  call assert(diag_field_id < diag_cs%next_free_diag_id, &
1487  'post_data_3d: Unregistered diagnostic id')
1488  diag => diag_cs%diags(diag_field_id)
1489  do while (associated(diag))
1490  call assert(associated(diag%axes), 'post_data_3d: axes is not associated')
1491 
1492  staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point
1493  staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point
1494 
1495  if (diag%v_extensive .and. .not.diag%axes%is_native) then
1496  ! The field is vertically integrated and needs to be re-gridded
1497  if (present(mask)) then
1498  call mom_error(fatal,"post_data_3d: no mask for regridded field.")
1499  endif
1500 
1501  if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)
1502  allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz))
1503  call vertically_reintegrate_diag_field( &
1504  diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), &
1505  diag_cs%G, h_diag, staggered_in_x, staggered_in_y, &
1506  diag%axes%mask3d, diag_cs%missing_value, field, remapped_field)
1507  if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)
1508  if (associated(diag%axes%mask3d)) then
1509  ! Since 3d masks do not vary in the vertical, just use as much as is
1510  ! needed.
1511  call post_data_3d_low(diag, remapped_field, diag_cs, is_static, &
1512  mask=diag%axes%mask3d)
1513  else
1514  call post_data_3d_low(diag, remapped_field, diag_cs, is_static)
1515  endif
1516  if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)
1517  deallocate(remapped_field)
1518  if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)
1519  elseif (diag%axes%needs_remapping) then
1520  ! Remap this field to another vertical coordinate.
1521  if (present(mask)) then
1522  call mom_error(fatal,"post_data_3d: no mask for regridded field.")
1523  endif
1524 
1525  if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)
1526  allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz))
1527  call diag_remap_do_remap(diag_cs%diag_remap_cs( &
1528  diag%axes%vertical_coordinate_number), &
1529  diag_cs%G, diag_cs%GV, h_diag, staggered_in_x, staggered_in_y, &
1530  diag%axes%mask3d, diag_cs%missing_value, field, remapped_field)
1531  if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)
1532  if (associated(diag%axes%mask3d)) then
1533  ! Since 3d masks do not vary in the vertical, just use as much as is
1534  ! needed.
1535  call post_data_3d_low(diag, remapped_field, diag_cs, is_static, &
1536  mask=diag%axes%mask3d)
1537  else
1538  call post_data_3d_low(diag, remapped_field, diag_cs, is_static)
1539  endif
1540  if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)
1541  deallocate(remapped_field)
1542  if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)
1543  elseif (diag%axes%needs_interpolating) then
1544  ! Interpolate this field to another vertical coordinate.
1545  if (present(mask)) then
1546  call mom_error(fatal,"post_data_3d: no mask for regridded field.")
1547  endif
1548 
1549  if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)
1550  allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz+1))
1551  call vertically_interpolate_diag_field(diag_cs%diag_remap_cs( &
1552  diag%axes%vertical_coordinate_number), &
1553  diag_cs%G, h_diag, staggered_in_x, staggered_in_y, &
1554  diag%axes%mask3d, diag_cs%missing_value, field, remapped_field)
1555  if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)
1556  if (associated(diag%axes%mask3d)) then
1557  ! Since 3d masks do not vary in the vertical, just use as much as is
1558  ! needed.
1559  call post_data_3d_low(diag, remapped_field, diag_cs, is_static, &
1560  mask=diag%axes%mask3d)
1561  else
1562  call post_data_3d_low(diag, remapped_field, diag_cs, is_static)
1563  endif
1564  if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)
1565  deallocate(remapped_field)
1566  if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)
1567  else
1568  call post_data_3d_low(diag, field, diag_cs, is_static, mask)
1569  endif
1570  diag => diag%next
1571  enddo
1572  if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator)
1573 
1574 end subroutine post_data_3d
1575 
1576 !> Make a real 3-d array diagnostic available for averaging or output
1577 !! using a diag_type instead of an integer id.
1578 subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask)
1579  type(diag_type), intent(in) :: diag !< A structure describing the diagnostic to post
1580  real, target, intent(in) :: field(:,:,:) !< 3-d array being offered for output or averaging
1581  type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output
1582  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.
1583  real, optional,target, intent(in) :: mask(:,:,:) !< If present, use this real array as the data mask.
1584 
1585  ! Local variables
1586  real, dimension(:,:,:), pointer :: locfield
1587  real, dimension(:,:,:), pointer :: locmask
1588  character(len=300) :: mesg
1589  logical :: used ! The return value of send_data is not used for anything.
1590  logical :: staggered_in_x, staggered_in_y
1591  logical :: is_stat
1592  integer :: cszi, cszj, dszi, dszj
1593  integer :: isv, iev, jsv, jev, ks, ke, i, j, k, isv_c, jsv_c, isv_o,jsv_o
1594  integer :: chksum
1595  real, dimension(:,:,:), allocatable, target :: locfield_dsamp
1596  real, dimension(:,:,:), allocatable, target :: locmask_dsamp
1597  integer :: dl
1598 
1599  locfield => null()
1600  locmask => null()
1601  is_stat = .false. ; if (present(is_static)) is_stat = is_static
1602 
1603  ! Determine the proper array indices, noting that because of the (:,:)
1604  ! declaration of field, symmetric arrays are using a SW-grid indexing,
1605  ! but non-symmetric arrays are using a NE-grid indexing. Send_data
1606  ! actually only uses the difference between ie and is to determine
1607  ! the output data size and assumes that halos are symmetric.
1608  isv = diag_cs%is ; iev = diag_cs%ie ; jsv = diag_cs%js ; jev = diag_cs%je
1609 
1610  cszi = (diag_cs%ie-diag_cs%is) +1 ; dszi = (diag_cs%ied-diag_cs%isd) +1
1611  cszj = (diag_cs%je-diag_cs%js) +1 ; dszj = (diag_cs%jed-diag_cs%jsd) +1
1612  if ( size(field,1) == dszi ) then
1613  isv = diag_cs%is ; iev = diag_cs%ie ! Data domain
1614  elseif ( size(field,1) == dszi + 1 ) then
1615  isv = diag_cs%is ; iev = diag_cs%ie+1 ! Symmetric data domain
1616  elseif ( size(field,1) == cszi) then
1617  isv = 1 ; iev = cszi ! Computational domain
1618  elseif ( size(field,1) == cszi + 1 ) then
1619  isv = 1 ; iev = cszi+1 ! Symmetric computational domain
1620  else
1621  write (mesg,*) " peculiar size ",size(field,1)," in i-direction\n"//&
1622  "does not match one of ", cszi, cszi+1, dszi, dszi+1
1623  call mom_error(fatal,"post_data_3d_low: "//trim(diag%debug_str)//trim(mesg))
1624  endif
1625 
1626  if ( size(field,2) == dszj ) then
1627  jsv = diag_cs%js ; jev = diag_cs%je ! Data domain
1628  elseif ( size(field,2) == dszj + 1 ) then
1629  jsv = diag_cs%js ; jev = diag_cs%je+1 ! Symmetric data domain
1630  elseif ( size(field,2) == cszj ) then
1631  jsv = 1 ; jev = cszj ! Computational domain
1632  elseif ( size(field,2) == cszj+1 ) then
1633  jsv = 1 ; jev = cszj+1 ! Symmetric computational domain
1634  else
1635  write (mesg,*) " peculiar size ",size(field,2)," in j-direction\n"//&
1636  "does not match one of ", cszj, cszj+1, dszj, dszj+1
1637  call mom_error(fatal,"post_data_3d_low: "//trim(diag%debug_str)//trim(mesg))
1638  endif
1639 
1640  ks = lbound(field,3) ; ke = ubound(field,3)
1641  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) then
1642  allocate( locfield( lbound(field,1):ubound(field,1), lbound(field,2):ubound(field,2), ks:ke ) )
1643  ! locfield(:,:,:) = 0.0 ! Zeroing out this array would be a good idea, but it appears
1644  ! not to be necessary.
1645  isv_c = isv ; jsv_c = jsv
1646  if (diag%fms_xyave_diag_id>0) then
1647  staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point
1648  staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point
1649  ! When averaging a staggered field, edge points are always required.
1650  if (staggered_in_x) isv_c = iev - (diag_cs%ie - diag_cs%is) - 1
1651  if (staggered_in_y) jsv_c = jev - (diag_cs%je - diag_cs%js) - 1
1652  if (isv_c < lbound(locfield,1)) call mom_error(fatal, &
1653  "It is an error to average a staggered diagnostic field that does not "//&
1654  "have i-direction space to represent the symmetric computational domain.")
1655  if (jsv_c < lbound(locfield,2)) call mom_error(fatal, &
1656  "It is an error to average a staggered diagnostic field that does not "//&
1657  "have j-direction space to represent the symmetric computational domain.")
1658  endif
1659 
1660  do k=ks,ke ; do j=jsv,jev ; do i=isv,iev
1661  if (field(i,j,k) == diag_cs%missing_value) then
1662  locfield(i,j,k) = diag_cs%missing_value
1663  else
1664  locfield(i,j,k) = field(i,j,k) * diag%conversion_factor
1665  endif
1666  enddo ; enddo ; enddo
1667  else
1668  locfield => field
1669  endif
1670 
1671  if (present(mask)) then
1672  locmask => mask
1673  elseif(associated(diag%axes%mask3d)) then
1674  locmask => diag%axes%mask3d
1675  endif
1676 
1677  dl=1
1678  if(.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet
1679  !Downsample the diag field and mask (if present)
1680  if (dl > 1) then
1681  isv_o=isv ; jsv_o=jsv
1682  call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask)
1683  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield )
1684  locfield => locfield_dsamp
1685  if (present(mask)) then
1686  call downsample_field_3d(locmask, locmask_dsamp, dl, msk, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev)
1687  locmask => locmask_dsamp
1688  elseif(associated(diag%axes%dsamp(dl)%mask3d)) then
1689  locmask => diag%axes%dsamp(dl)%mask3d
1690  endif
1691  endif
1692 
1693  if (diag%fms_diag_id>0) then
1694  if (diag_cs%diag_as_chksum) then
1695  if (diag%axes%is_h_point) then
1696  call hchksum(locfield, diag%debug_str, diag_cs%G%HI, &
1697  logunit=diag_cs%chksum_iounit)
1698  else if (diag%axes%is_u_point) then
1699  call uchksum(locfield, diag%debug_str, diag_cs%G%HI, &
1700  logunit=diag_cs%chksum_iounit)
1701  else if (diag%axes%is_v_point) then
1702  call vchksum(locfield, diag%debug_str, diag_cs%G%HI, &
1703  logunit=diag_cs%chksum_iounit)
1704  else if (diag%axes%is_q_point) then
1705  call bchksum(locfield, diag%debug_str, diag_cs%G%HI, &
1706  logunit=diag_cs%chksum_iounit)
1707  else
1708  call mom_error(fatal, "post_data_3d_low: unknown axis type.")
1709  endif
1710  else
1711  if (is_stat) then
1712  if (present(mask)) then
1713  call assert(size(locfield) == size(locmask), &
1714  'post_data_3d_low is_stat: mask size mismatch: '//diag%debug_str)
1715  used = send_data(diag%fms_diag_id, locfield, &
1716  is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=locmask)
1717  !elseif (associated(diag%axes%mask2d)) then
1718  ! used = send_data(diag%fms_diag_id, locfield, &
1719  ! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%axes%mask2d)
1720  else
1721  used = send_data(diag%fms_diag_id, locfield, &
1722  is_in=isv, js_in=jsv, ie_in=iev, je_in=jev)
1723  endif
1724  elseif (diag_cs%ave_enabled) then
1725  if (associated(locmask)) then
1726  call assert(size(locfield) == size(locmask), &
1727  'post_data_3d_low: mask size mismatch: '//diag%debug_str)
1728  used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, &
1729  is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, &
1730  weight=diag_cs%time_int, rmask=locmask)
1731  else
1732  used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, &
1733  is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, &
1734  weight=diag_cs%time_int)
1735  endif
1736  endif
1737  endif
1738  endif
1739 
1740  if (diag%fms_xyave_diag_id>0) then
1741  call post_xy_average(diag_cs, diag, locfield)
1742  endif
1743 
1744  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.) .and. dl<2) &
1745  deallocate( locfield )
1746 
1747 end subroutine post_data_3d_low
1748 
1749 !> Post the horizontally area-averaged diagnostic
1750 subroutine post_xy_average(diag_cs, diag, field)
1751  type(diag_type), intent(in) :: diag !< This diagnostic
1752  real, target, intent(in) :: field(:,:,:) !< Diagnostic field
1753  type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics mediator control structure
1754  ! Local variable
1755  real, dimension(size(field,3)) :: averaged_field
1756  logical, dimension(size(field,3)) :: averaged_mask
1757  logical :: staggered_in_x, staggered_in_y, used
1758  integer :: nz, remap_nz, coord
1759 
1760  if (.not. diag_cs%ave_enabled) then
1761  return
1762  endif
1763 
1764  staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point
1765  staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point
1766 
1767  if (diag%axes%is_native) then
1768  call horizontally_average_diag_field(diag_cs%G, diag_cs%GV, diag_cs%h, &
1769  staggered_in_x, staggered_in_y, &
1770  diag%axes%is_layer, diag%v_extensive, &
1771  diag_cs%missing_value, field, &
1772  averaged_field, averaged_mask)
1773  else
1774  nz = size(field, 3)
1775  coord = diag%axes%vertical_coordinate_number
1776  remap_nz = diag_cs%diag_remap_cs(coord)%nz
1777 
1778  call assert(diag_cs%diag_remap_cs(coord)%initialized, &
1779  'post_xy_average: remap_cs not initialized.')
1780 
1781  call assert(implies(diag%axes%is_layer, nz == remap_nz), &
1782  'post_xy_average: layer field dimension mismatch.')
1783  call assert(implies(.not. diag%axes%is_layer, nz == remap_nz+1), &
1784  'post_xy_average: interface field dimension mismatch.')
1785 
1786  call horizontally_average_diag_field(diag_cs%G, diag_cs%GV, &
1787  diag_cs%diag_remap_cs(coord)%h, &
1788  staggered_in_x, staggered_in_y, &
1789  diag%axes%is_layer, diag%v_extensive, &
1790  diag_cs%missing_value, field, &
1791  averaged_field, averaged_mask)
1792  endif
1793 
1794  if (diag_cs%diag_as_chksum) then
1795  call zchksum(averaged_field, trim(diag%debug_str)//'_xyave', &
1796  logunit=diag_cs%chksum_iounit)
1797  else
1798  used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, &
1799  weight=diag_cs%time_int, mask=averaged_mask)
1800  endif
1801 end subroutine post_xy_average
1802 
1803 !> This subroutine enables the accumulation of time averages over the specified time interval.
1804 subroutine enable_averaging(time_int_in, time_end_in, diag_cs)
1805  real, intent(in) :: time_int_in !< The time interval [s] over which any
1806  !! values that are offered are valid.
1807  type(time_type), intent(in) :: time_end_in !< The end time of the valid interval
1808  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output
1809 
1810 ! This subroutine enables the accumulation of time averages over the specified time interval.
1811 
1812 ! if (num_file==0) return
1813  diag_cs%time_int = time_int_in
1814  diag_cs%time_end = time_end_in
1815  diag_cs%ave_enabled = .true.
1816 end subroutine enable_averaging
1817 
1818 !> Enable the accumulation of time averages over the specified time interval in time units.
1819 subroutine enable_averages(time_int, time_end, diag_CS, T_to_s)
1820  real, intent(in) :: time_int !< The time interval over which any values
1821  !! that are offered are valid [T ~> s].
1822  type(time_type), intent(in) :: time_end !< The end time of the valid interval.
1823  type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output
1824  real, optional, intent(in) :: t_to_s !< A conversion factor for time_int to [s].
1825 ! This subroutine enables the accumulation of time averages over the specified time interval.
1826 
1827  if (present(t_to_s)) then
1828  diag_cs%time_int = time_int*t_to_s
1829  elseif (associated(diag_cs%US)) then
1830  diag_cs%time_int = time_int*diag_cs%US%T_to_s
1831  else
1832  diag_cs%time_int = time_int
1833  endif
1834  diag_cs%time_end = time_end
1835  diag_cs%ave_enabled = .true.
1836 end subroutine enable_averages
1837 
1838 !> Call this subroutine to avoid averaging any offered fields.
1839 subroutine disable_averaging(diag_cs)
1840  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output
1841 
1842  diag_cs%time_int = 0.0
1843  diag_cs%ave_enabled = .false.
1844 
1845 end subroutine disable_averaging
1846 
1847 !> Call this subroutine to determine whether the averaging is
1848 !! currently enabled. .true. is returned if it is.
1849 function query_averaging_enabled(diag_cs, time_int, time_end)
1850  type(diag_ctrl), intent(in) :: diag_cs !< Structure used to regulate diagnostic output
1851  real, optional, intent(out) :: time_int !< Current setting of diag%time_int [s]
1852  type(time_type), optional, intent(out) :: time_end !< Current setting of diag%time_end
1853  logical :: query_averaging_enabled
1854 
1855  if (present(time_int)) time_int = diag_cs%time_int
1856  if (present(time_end)) time_end = diag_cs%time_end
1857  query_averaging_enabled = diag_cs%ave_enabled
1858 end function query_averaging_enabled
1859 
1860 !> This function returns the valid end time for use with diagnostics that are
1861 !! handled outside of the MOM6 diagnostics infrastructure.
1862 function get_diag_time_end(diag_cs)
1863  type(diag_ctrl), intent(in) :: diag_cs !< Structure used to regulate diagnostic output
1864  type(time_type) :: get_diag_time_end
1865  ! This function returns the valid end time for diagnostics that are handled
1866  ! outside of the MOM6 infrastructure, such as via the generic tracer code.
1867 
1868  get_diag_time_end = diag_cs%time_end
1869 end function get_diag_time_end
1870 
1871 !> Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics
1872 !! derived from one field.
1873 integer function register_diag_field(module_name, field_name, axes_in, init_time, &
1874  long_name, units, missing_value, range, mask_variant, standard_name, &
1875  verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, &
1876  cmor_long_name, cmor_units, cmor_standard_name, cell_methods, &
1877  x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
1878  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model"
1879  !! or "ice_shelf_model"
1880  character(len=*), intent(in) :: field_name !< Name of the diagnostic field
1881  type(axes_grp), target, intent(in) :: axes_in !< Container w/ up to 3 integer handles that
1882  !! indicates axes for this field
1883  type(time_type), intent(in) :: init_time !< Time at which a field is first available?
1884  character(len=*), optional, intent(in) :: long_name !< Long name of a field.
1885  character(len=*), optional, intent(in) :: units !< Units of a field.
1886  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field
1887  real, optional, intent(in) :: missing_value !< A value that indicates missing values.
1888  real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)
1889  logical, optional, intent(in) :: mask_variant !< If true a logical mask must be provided with
1890  !! post_data calls (not used in MOM?)
1891  logical, optional, intent(in) :: verbose !< If true, FMS is verbose (not used in MOM?)
1892  logical, optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?)
1893  character(len=*), optional, intent(out):: err_msg !< String into which an error message might be
1894  !! placed (not used in MOM?)
1895  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not
1896  !! be interpolated as a scalar
1897  integer, optional, intent(in) :: tile_count !< no clue (not used in MOM?)
1898  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field
1899  character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field
1900  character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field
1901  character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field
1902  character(len=*), optional, intent(in) :: cell_methods !< String to append as cell_methods attribute. Use '' to
1903  !! have no attribute. If present, this overrides the
1904  !! default constructed from the default for
1905  !! each individual axis direction.
1906  character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction.
1907  !! Use '' have no method.
1908  character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction.
1909  !! Use '' have no method.
1910  character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction.
1911  !! Use '' have no method.
1912  real, optional, intent(in) :: conversion !< A value to multiply data by before writing to file
1913  logical, optional, intent(in) :: v_extensive !< True for vertically extensive fields (vertically
1914  !! integrated). Default/absent for intensive.
1915  ! Local variables
1916  real :: mom_missing_value
1917  type(diag_ctrl), pointer :: diag_cs => null()
1918  type(axes_grp), pointer :: remap_axes => null()
1919  type(axes_grp), pointer :: axes => null()
1920  integer :: dm_id, i, dl
1921  character(len=256) :: new_module_name
1922  logical :: active
1923 
1924  axes => axes_in
1925  mom_missing_value = axes%diag_cs%missing_value
1926  if (present(missing_value)) mom_missing_value = missing_value
1927 
1928  diag_cs => axes%diag_cs
1929  dm_id = -1
1930 
1931  if (axes_in%id == diag_cs%axesTL%id) then
1932  axes => diag_cs%axesTL
1933  elseif (axes_in%id == diag_cs%axesBL%id) then
1934  axes => diag_cs%axesBL
1935  elseif (axes_in%id == diag_cs%axesCuL%id ) then
1936  axes => diag_cs%axesCuL
1937  elseif (axes_in%id == diag_cs%axesCvL%id) then
1938  axes => diag_cs%axesCvL
1939  elseif (axes_in%id == diag_cs%axesTi%id) then
1940  axes => diag_cs%axesTi
1941  elseif (axes_in%id == diag_cs%axesBi%id) then
1942  axes => diag_cs%axesBi
1943  elseif (axes_in%id == diag_cs%axesCui%id ) then
1944  axes => diag_cs%axesCui
1945  elseif (axes_in%id == diag_cs%axesCvi%id) then
1946  axes => diag_cs%axesCvi
1947  endif
1948 
1949  ! Register the native diagnostic
1950  active = register_diag_field_expand_cmor(dm_id, module_name, field_name, axes, &
1951  init_time, long_name=long_name, units=units, missing_value=mom_missing_value, &
1952  range=range, mask_variant=mask_variant, standard_name=standard_name, &
1953  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
1954  interp_method=interp_method, tile_count=tile_count, &
1955  cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, &
1956  cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, &
1957  cell_methods=cell_methods, x_cell_method=x_cell_method, &
1958  y_cell_method=y_cell_method, v_cell_method=v_cell_method, &
1959  conversion=conversion, v_extensive=v_extensive)
1960 
1961  ! For each diagnostic coordinate register the diagnostic again under a different module name
1962  do i=1,diag_cs%num_diag_coords
1963  new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)
1964 
1965  ! Register diagnostics remapped to z vertical coordinate
1966  if (axes_in%rank == 3) then
1967  remap_axes => null()
1968  if ((axes_in%id == diag_cs%axesTL%id)) then
1969  remap_axes => diag_cs%remap_axesTL(i)
1970  elseif (axes_in%id == diag_cs%axesBL%id) then
1971  remap_axes => diag_cs%remap_axesBL(i)
1972  elseif (axes_in%id == diag_cs%axesCuL%id ) then
1973  remap_axes => diag_cs%remap_axesCuL(i)
1974  elseif (axes_in%id == diag_cs%axesCvL%id) then
1975  remap_axes => diag_cs%remap_axesCvL(i)
1976  elseif (axes_in%id == diag_cs%axesTi%id) then
1977  remap_axes => diag_cs%remap_axesTi(i)
1978  elseif (axes_in%id == diag_cs%axesBi%id) then
1979  remap_axes => diag_cs%remap_axesBi(i)
1980  elseif (axes_in%id == diag_cs%axesCui%id ) then
1981  remap_axes => diag_cs%remap_axesCui(i)
1982  elseif (axes_in%id == diag_cs%axesCvi%id) then
1983  remap_axes => diag_cs%remap_axesCvi(i)
1984  endif
1985  ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will
1986  ! always exist but in the mean-time we have to do this check:
1987  ! call assert(associated(remap_axes), 'register_diag_field: remap_axes not set')
1988  if (associated(remap_axes)) then
1989  if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then
1990  active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, remap_axes, &
1991  init_time, long_name=long_name, units=units, missing_value=mom_missing_value, &
1992  range=range, mask_variant=mask_variant, standard_name=standard_name, &
1993  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
1994  interp_method=interp_method, tile_count=tile_count, &
1995  cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, &
1996  cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, &
1997  cell_methods=cell_methods, x_cell_method=x_cell_method, &
1998  y_cell_method=y_cell_method, v_cell_method=v_cell_method, &
1999  conversion=conversion, v_extensive=v_extensive)
2000  if (active) then
2001  call diag_remap_set_active(diag_cs%diag_remap_cs(i))
2002  endif
2003  endif ! remap_axes%needs_remapping
2004  endif ! associated(remap_axes)
2005  endif ! axes%rank == 3
2006  enddo ! i
2007 
2008  !Register downsampled diagnostics
2009  do dl=2,max_dsamp_lev
2010  ! Do not attempt to checksum the downsampled diagnostics
2011  if (diag_cs%diag_as_chksum) cycle
2012 
2013  new_module_name = trim(module_name)//'_d2'
2014 
2015  if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then
2016  axes => null()
2017  if (axes_in%id == diag_cs%axesTL%id) then
2018  axes => diag_cs%dsamp(dl)%axesTL
2019  elseif (axes_in%id == diag_cs%axesBL%id) then
2020  axes => diag_cs%dsamp(dl)%axesBL
2021  elseif (axes_in%id == diag_cs%axesCuL%id ) then
2022  axes => diag_cs%dsamp(dl)%axesCuL
2023  elseif (axes_in%id == diag_cs%axesCvL%id) then
2024  axes => diag_cs%dsamp(dl)%axesCvL
2025  elseif (axes_in%id == diag_cs%axesTi%id) then
2026  axes => diag_cs%dsamp(dl)%axesTi
2027  elseif (axes_in%id == diag_cs%axesBi%id) then
2028  axes => diag_cs%dsamp(dl)%axesBi
2029  elseif (axes_in%id == diag_cs%axesCui%id ) then
2030  axes => diag_cs%dsamp(dl)%axesCui
2031  elseif (axes_in%id == diag_cs%axesCvi%id) then
2032  axes => diag_cs%dsamp(dl)%axesCvi
2033  elseif (axes_in%id == diag_cs%axesT1%id) then
2034  axes => diag_cs%dsamp(dl)%axesT1
2035  elseif (axes_in%id == diag_cs%axesB1%id) then
2036  axes => diag_cs%dsamp(dl)%axesB1
2037  elseif (axes_in%id == diag_cs%axesCu1%id ) then
2038  axes => diag_cs%dsamp(dl)%axesCu1
2039  elseif (axes_in%id == diag_cs%axesCv1%id) then
2040  axes => diag_cs%dsamp(dl)%axesCv1
2041  else
2042  !Niki: Should we worry about these, e.g., diag_to_Z_CS?
2043  call mom_error(warning,"register_diag_field: Could not find a proper axes for " &
2044  //trim( new_module_name)//"-"//trim(field_name))
2045  endif
2046  endif
2047  ! Register the native diagnostic
2048  if (associated(axes)) then
2049  active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes, &
2050  init_time, long_name=long_name, units=units, missing_value=mom_missing_value, &
2051  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2052  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2053  interp_method=interp_method, tile_count=tile_count, &
2054  cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, &
2055  cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, &
2056  cell_methods=cell_methods, x_cell_method=x_cell_method, &
2057  y_cell_method=y_cell_method, v_cell_method=v_cell_method, &
2058  conversion=conversion, v_extensive=v_extensive)
2059  endif
2060 
2061  ! For each diagnostic coordinate register the diagnostic again under a different module name
2062  do i=1,diag_cs%num_diag_coords
2063  new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)//'_d2'
2064 
2065  ! Register diagnostics remapped to z vertical coordinate
2066  if (axes_in%rank == 3) then
2067  remap_axes => null()
2068  if ((axes_in%id == diag_cs%axesTL%id)) then
2069  remap_axes => diag_cs%dsamp(dl)%remap_axesTL(i)
2070  elseif (axes_in%id == diag_cs%axesBL%id) then
2071  remap_axes => diag_cs%dsamp(dl)%remap_axesBL(i)
2072  elseif (axes_in%id == diag_cs%axesCuL%id ) then
2073  remap_axes => diag_cs%dsamp(dl)%remap_axesCuL(i)
2074  elseif (axes_in%id == diag_cs%axesCvL%id) then
2075  remap_axes => diag_cs%dsamp(dl)%remap_axesCvL(i)
2076  elseif (axes_in%id == diag_cs%axesTi%id) then
2077  remap_axes => diag_cs%dsamp(dl)%remap_axesTi(i)
2078  elseif (axes_in%id == diag_cs%axesBi%id) then
2079  remap_axes => diag_cs%dsamp(dl)%remap_axesBi(i)
2080  elseif (axes_in%id == diag_cs%axesCui%id ) then
2081  remap_axes => diag_cs%dsamp(dl)%remap_axesCui(i)
2082  elseif (axes_in%id == diag_cs%axesCvi%id) then
2083  remap_axes => diag_cs%dsamp(dl)%remap_axesCvi(i)
2084  endif
2085 
2086  ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will
2087  ! always exist but in the mean-time we have to do this check:
2088  ! call assert(associated(remap_axes), 'register_diag_field: remap_axes not set')
2089  if (associated(remap_axes)) then
2090  if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then
2091  active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, remap_axes, &
2092  init_time, long_name=long_name, units=units, missing_value=mom_missing_value, &
2093  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2094  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2095  interp_method=interp_method, tile_count=tile_count, &
2096  cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, &
2097  cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, &
2098  cell_methods=cell_methods, x_cell_method=x_cell_method, &
2099  y_cell_method=y_cell_method, v_cell_method=v_cell_method, &
2100  conversion=conversion, v_extensive=v_extensive)
2101  if (active) then
2102  call diag_remap_set_active(diag_cs%diag_remap_cs(i))
2103  endif
2104  endif ! remap_axes%needs_remapping
2105  endif ! associated(remap_axes)
2106  endif ! axes%rank == 3
2107  enddo ! i
2108  enddo
2109 
2110  register_diag_field = dm_id
2111 
2112 end function register_diag_field
2113 
2114 !> Returns True if either the native or CMOr version of the diagnostic were registered. Updates 'dm_id'
2115 !! after calling register_diag_field_expand_axes() for both native and CMOR variants of the field.
2116 logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, axes, init_time, &
2117  long_name, units, missing_value, range, mask_variant, standard_name, &
2118  verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, &
2119  cmor_long_name, cmor_units, cmor_standard_name, cell_methods, &
2120  x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
2121  integer, intent(inout) :: dm_id !< The diag_mediator ID for this diagnostic group
2122  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" or "ice_shelf_model"
2123  character(len=*), intent(in) :: field_name !< Name of the diagnostic field
2124  type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes
2125  !! for this field
2126  type(time_type), intent(in) :: init_time !< Time at which a field is first available?
2127  character(len=*), optional, intent(in) :: long_name !< Long name of a field.
2128  character(len=*), optional, intent(in) :: units !< Units of a field.
2129  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field
2130  real, optional, intent(in) :: missing_value !< A value that indicates missing values.
2131  real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)
2132  logical, optional, intent(in) :: mask_variant !< If true a logical mask must be provided
2133  !! with post_data calls (not used in MOM?)
2134  logical, optional, intent(in) :: verbose !< If true, FMS is verbose (not used in MOM?)
2135  logical, optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?)
2136  character(len=*), optional, intent(out):: err_msg !< String into which an error message might be
2137  !! placed (not used in MOM?)
2138  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should
2139  !! not be interpolated as a scalar
2140  integer, optional, intent(in) :: tile_count !< no clue (not used in MOM?)
2141  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field
2142  character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field
2143  character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field
2144  character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field
2145  character(len=*), optional, intent(in) :: cell_methods !< String to append as cell_methods attribute.
2146  !! Use '' to have no attribute. If present, this
2147  !! overrides the default constructed from the default
2148  !! for each individual axis direction.
2149  character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction.
2150  !! Use '' have no method.
2151  character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction.
2152  !! Use '' have no method.
2153  character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction.
2154  !! Use '' have no method.
2155  real, optional, intent(in) :: conversion !< A value to multiply data by before writing to file
2156  logical, optional, intent(in) :: v_extensive !< True for vertically extensive fields (vertically
2157  !! integrated). Default/absent for intensive.
2158  ! Local variables
2159  real :: mom_missing_value
2160  type(diag_ctrl), pointer :: diag_cs => null()
2161  type(diag_type), pointer :: this_diag => null()
2162  integer :: fms_id, fms_xyave_id
2163  character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name, cm_string, msg
2164 
2165  mom_missing_value = axes%diag_cs%missing_value
2166  if (present(missing_value)) mom_missing_value = missing_value
2167 
2169  diag_cs => axes%diag_cs
2170 
2171  ! Set up the 'primary' diagnostic, first get an underlying FMS id
2172  fms_id = register_diag_field_expand_axes(module_name, field_name, axes, init_time, &
2173  long_name=long_name, units=units, missing_value=mom_missing_value, &
2174  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2175  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2176  interp_method=interp_method, tile_count=tile_count)
2177  if (.not. diag_cs%diag_as_chksum) &
2178  call attach_cell_methods(fms_id, axes, cm_string, cell_methods, &
2179  x_cell_method, y_cell_method, v_cell_method, &
2180  v_extensive=v_extensive)
2181  if (is_root_pe() .and. diag_cs%available_diag_doc_unit > 0) then
2182  msg = ''
2183  if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"'
2184  call log_available_diag(fms_id>0, module_name, field_name, cm_string, &
2185  msg, diag_cs, long_name, units, standard_name)
2186  endif
2187  ! Associated horizontally area-averaged diagnostic
2188  fms_xyave_id = diag_field_not_found
2189  if (associated(axes%xyave_axes)) then
2190  fms_xyave_id = register_diag_field_expand_axes(module_name, trim(field_name)//'_xyave', &
2191  axes%xyave_axes, init_time, &
2192  long_name=long_name, units=units, missing_value=mom_missing_value, &
2193  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2194  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2195  interp_method=interp_method, tile_count=tile_count)
2196  if (.not. diag_cs%diag_as_chksum) &
2197  call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, &
2198  cell_methods, v_cell_method, v_extensive=v_extensive)
2199  if (is_root_pe() .and. diag_cs%available_diag_doc_unit > 0) then
2200  msg = ''
2201  if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'_xyave"'
2202  call log_available_diag(fms_xyave_id>0, module_name, trim(field_name)//'_xyave', cm_string, &
2203  msg, diag_cs, long_name, units, standard_name)
2204  endif
2205  endif
2206  this_diag => null()
2207  if (fms_id /= diag_field_not_found .or. fms_xyave_id /= diag_field_not_found) then
2208  call add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name, msg)
2209  this_diag%fms_xyave_diag_id = fms_xyave_id
2210  !Encode and save the cell methods for this diag
2211  call add_xyz_method(this_diag, axes, x_cell_method, y_cell_method, v_cell_method, v_extensive)
2212  if (present(v_extensive)) this_diag%v_extensive = v_extensive
2213  if (present(conversion)) this_diag%conversion_factor = conversion
2215  endif
2216 
2217  ! For the CMOR variation of the above diagnostic
2218  if (present(cmor_field_name) .and. .not. diag_cs%diag_as_chksum) then
2219  ! Fallback values for strings set to "NULL"
2220  posted_cmor_units = "not provided" !
2221  posted_cmor_standard_name = "not provided" ! Values might be able to be replaced with a CS%missing field?
2222  posted_cmor_long_name = "not provided" !
2223 
2224  ! If attributes are present for MOM variable names, use them first for the register_diag_field
2225  ! call for CMOR verison of the variable
2226  if (present(units)) posted_cmor_units = units
2227  if (present(standard_name)) posted_cmor_standard_name = standard_name
2228  if (present(long_name)) posted_cmor_long_name = long_name
2229 
2230  ! If specified in the call to register_diag_field, override attributes with the CMOR versions
2231  if (present(cmor_units)) posted_cmor_units = cmor_units
2232  if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name
2233  if (present(cmor_long_name)) posted_cmor_long_name = cmor_long_name
2234 
2235  fms_id = register_diag_field_expand_axes(module_name, cmor_field_name, axes, init_time, &
2236  long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), &
2237  missing_value=mom_missing_value, range=range, mask_variant=mask_variant, &
2238  standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log, &
2239  err_msg=err_msg, interp_method=interp_method, tile_count=tile_count)
2240  call attach_cell_methods(fms_id, axes, cm_string, &
2241  cell_methods, x_cell_method, y_cell_method, v_cell_method, &
2242  v_extensive=v_extensive)
2243  if (is_root_pe() .and. diag_cs%available_diag_doc_unit > 0) then
2244  msg = 'native name is "'//trim(field_name)//'"'
2245  call log_available_diag(fms_id>0, module_name, cmor_field_name, cm_string, &
2246  msg, diag_cs, posted_cmor_long_name, posted_cmor_units, &
2247  posted_cmor_standard_name)
2248  endif
2249  ! Associated horizontally area-averaged diagnostic
2250  fms_xyave_id = diag_field_not_found
2251  if (associated(axes%xyave_axes)) then
2252  fms_xyave_id = register_diag_field_expand_axes(module_name, trim(cmor_field_name)//'_xyave', &
2253  axes%xyave_axes, init_time, &
2254  long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), &
2255  missing_value=mom_missing_value, range=range, mask_variant=mask_variant, &
2256  standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log, &
2257  err_msg=err_msg, interp_method=interp_method, tile_count=tile_count)
2258  call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, &
2259  cell_methods, v_cell_method, v_extensive=v_extensive)
2260  if (is_root_pe() .and. diag_cs%available_diag_doc_unit > 0) then
2261  msg = 'native name is "'//trim(field_name)//'_xyave"'
2262  call log_available_diag(fms_xyave_id>0, module_name, trim(cmor_field_name)//'_xyave', &
2263  cm_string, msg, diag_cs, posted_cmor_long_name, posted_cmor_units, &
2264  posted_cmor_standard_name)
2265  endif
2266  endif
2267  this_diag => null()
2268  if (fms_id /= diag_field_not_found .or. fms_xyave_id /= diag_field_not_found) then
2269  call add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name, msg)
2270  this_diag%fms_xyave_diag_id = fms_xyave_id
2271  !Encode and save the cell methods for this diag
2272  call add_xyz_method(this_diag, axes, x_cell_method, y_cell_method, v_cell_method, v_extensive)
2273  if (present(v_extensive)) this_diag%v_extensive = v_extensive
2274  if (present(conversion)) this_diag%conversion_factor = conversion
2276  endif
2277  endif
2278 
2280 
2281 !> Returns an FMS id from register_diag_field_fms (the diag_manager routine) after expanding axes
2282 !! (axes-group) into handles and conditionally adding an FMS area_id for cell_measures.
2283 integer function register_diag_field_expand_axes(module_name, field_name, axes, init_time, &
2284  long_name, units, missing_value, range, mask_variant, standard_name, &
2285  verbose, do_not_log, err_msg, interp_method, tile_count)
2286  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model"
2287  !! or "ice_shelf_model"
2288  character(len=*), intent(in) :: field_name !< Name of the diagnostic field
2289  type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that indicates
2290  !! axes for this field
2291  type(time_type), intent(in) :: init_time !< Time at which a field is first available?
2292  character(len=*), optional, intent(in) :: long_name !< Long name of a field.
2293  character(len=*), optional, intent(in) :: units !< Units of a field.
2294  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field
2295  real, optional, intent(in) :: missing_value !< A value that indicates missing values.
2296  real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)
2297  logical, optional, intent(in) :: mask_variant !< If true a logical mask must be provided
2298  !! with post_data calls (not used in MOM?)
2299  logical, optional, intent(in) :: verbose !< If true, FMS is verbose (not used in MOM?)
2300  logical, optional, intent(in) :: do_not_log !< If true, do not log something
2301  !! (not used in MOM?)
2302  character(len=*), optional, intent(out):: err_msg !< String into which an error message might be
2303  !! placed (not used in MOM?)
2304  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should
2305  !! not be interpolated as a scalar
2306  integer, optional, intent(in) :: tile_count !< no clue (not used in MOM?)
2307  ! Local variables
2308  integer :: fms_id, area_id, volume_id
2309 
2310  ! This gets the cell area associated with the grid location of this variable
2311  area_id = axes%id_area
2312  volume_id = axes%id_volume
2313 
2314  ! Get the FMS diagnostic id
2315  if (axes%diag_cs%diag_as_chksum) then
2316  fms_id = axes%diag_cs%num_chksum_diags + 1
2317  axes%diag_cs%num_chksum_diags = fms_id
2318  else if (present(interp_method) .or. axes%is_h_point) then
2319  ! If interp_method is provided we must use it
2320  if (area_id>0) then
2321  if (volume_id>0) then
2322  fms_id = register_diag_field_fms(module_name, field_name, axes%handles, &
2323  init_time, long_name=long_name, units=units, missing_value=missing_value, &
2324  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2325  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2326  interp_method=interp_method, tile_count=tile_count, area=area_id, volume=volume_id)
2327  else
2328  fms_id = register_diag_field_fms(module_name, field_name, axes%handles, &
2329  init_time, long_name=long_name, units=units, missing_value=missing_value, &
2330  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2331  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2332  interp_method=interp_method, tile_count=tile_count, area=area_id)
2333  endif
2334  else
2335  if (volume_id>0) then
2336  fms_id = register_diag_field_fms(module_name, field_name, axes%handles, &
2337  init_time, long_name=long_name, units=units, missing_value=missing_value, &
2338  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2339  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2340  interp_method=interp_method, tile_count=tile_count, volume=volume_id)
2341  else
2342  fms_id = register_diag_field_fms(module_name, field_name, axes%handles, &
2343  init_time, long_name=long_name, units=units, missing_value=missing_value, &
2344  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2345  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2346  interp_method=interp_method, tile_count=tile_count)
2347  endif
2348  endif
2349  else
2350  ! If interp_method is not provided and the field is not at an h-point then interp_method='none'
2351  if (area_id>0) then
2352  if (volume_id>0) then
2353  fms_id = register_diag_field_fms(module_name, field_name, axes%handles, &
2354  init_time, long_name=long_name, units=units, missing_value=missing_value, &
2355  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2356  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2357  interp_method='none', tile_count=tile_count, area=area_id, volume=volume_id)
2358  else
2359  fms_id = register_diag_field_fms(module_name, field_name, axes%handles, &
2360  init_time, long_name=long_name, units=units, missing_value=missing_value, &
2361  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2362  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2363  interp_method='none', tile_count=tile_count, area=area_id)
2364  endif
2365  else
2366  if (volume_id>0) then
2367  fms_id = register_diag_field_fms(module_name, field_name, axes%handles, &
2368  init_time, long_name=long_name, units=units, missing_value=missing_value, &
2369  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2370  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2371  interp_method='none', tile_count=tile_count, volume=volume_id)
2372  else
2373  fms_id = register_diag_field_fms(module_name, field_name, axes%handles, &
2374  init_time, long_name=long_name, units=units, missing_value=missing_value, &
2375  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2376  verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
2377  interp_method='none', tile_count=tile_count)
2378  endif
2379  endif
2380  endif
2381 
2383 
2385 
2386 !> Create a diagnostic type and attached to list
2387 subroutine add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name, msg)
2388  type(diag_ctrl), pointer :: diag_cs !< Diagnostics mediator control structure
2389  integer, intent(inout) :: dm_id !< The diag_mediator ID for this diagnostic group
2390  integer, intent(in) :: fms_id !< The FMS diag_manager ID for this diagnostic
2391  type(diag_type), pointer :: this_diag !< This diagnostic
2392  type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that
2393  !! indicates axes for this field
2394  character(len=*), intent(in) :: module_name !< Name of this module, usually
2395  !! "ocean_model" or "ice_shelf_model"
2396  character(len=*), intent(in) :: field_name !< Name of diagnostic
2397  character(len=*), intent(in) :: msg !< Message for errors
2398 
2399  ! If the diagnostic is needed obtain a diag_mediator ID (if needed)
2400  if (dm_id == -1) dm_id = get_new_diag_id(diag_cs)
2401  ! Create a new diag_type to store links in
2402  call alloc_diag_with_id(dm_id, diag_cs, this_diag)
2403  call assert(associated(this_diag), trim(msg)//': diag_type allocation failed')
2404  ! Record FMS id, masks and conversion factor, in diag_type
2405  this_diag%fms_diag_id = fms_id
2406  this_diag%debug_str = trim(module_name)//"-"//trim(field_name)
2407  this_diag%axes => axes
2408 
2409 end subroutine add_diag_to_list
2410 
2411 !> Adds the encoded "cell_methods" for a diagnostics as a diag% property
2412 !! This allows access to the cell_method for a given diagnostics at the time of sending
2413 subroutine add_xyz_method(diag, axes, x_cell_method, y_cell_method, v_cell_method, v_extensive)
2414  type(diag_type), pointer :: diag !< This diagnostic
2415  type(axes_grp), intent(in) :: axes !< Container w/ up to 3 integer handles that indicates
2416  !! axes for this field
2417  character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction.
2418  !! Use '' have no method.
2419  character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction.
2420  !! Use '' have no method.
2421  character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction.
2422  !! Use '' have no method.
2423  logical, optional, intent(in) :: v_extensive !< True for vertically extensive fields
2424  !! (vertically integrated). Default/absent for intensive.
2425  integer :: xyz_method
2426  character(len=9) :: mstr
2427 
2428  !This is a simple way to encode the cell method information made from 3 strings
2429  !(x_cell_method,y_cell_method,v_cell_method) in a 3 digit integer xyz
2430  !x_cell_method,y_cell_method,v_cell_method can each be 'point' or 'sum' or 'mean'
2431  !We can encode these with setting 1 for 'point', 2 for 'sum, 3 for 'mean' in
2432  !the 100s position for x, 10s position for y, 1s position for z
2433  !E.g., x:sum,y:point,z:mean is 213
2434 
2435  xyz_method = 111
2436 
2437  mstr = diag%axes%v_cell_method
2438  if (present(v_extensive)) then
2439  if (present(v_cell_method)) call mom_error(fatal, "attach_cell_methods: " // &
2440  'Vertical cell method was specified along with the vertically extensive flag.')
2441  if(v_extensive) then
2442  mstr='sum'
2443  else
2444  mstr='mean'
2445  endif
2446  elseif (present(v_cell_method)) then
2447  mstr = v_cell_method
2448  endif
2449  if (trim(mstr)=='sum') then
2450  xyz_method = xyz_method + 1
2451  elseif (trim(mstr)=='mean') then
2452  xyz_method = xyz_method + 2
2453  endif
2454 
2455  mstr = diag%axes%y_cell_method
2456  if (present(y_cell_method)) mstr = y_cell_method
2457  if (trim(mstr)=='sum') then
2458  xyz_method = xyz_method + 10
2459  elseif (trim(mstr)=='mean') then
2460  xyz_method = xyz_method + 20
2461  endif
2462 
2463  mstr = diag%axes%x_cell_method
2464  if (present(x_cell_method)) mstr = x_cell_method
2465  if (trim(mstr)=='sum') then
2466  xyz_method = xyz_method + 100
2467  elseif (trim(mstr)=='mean') then
2468  xyz_method = xyz_method + 200
2469  endif
2470 
2471  diag%xyz_method = xyz_method
2472 end subroutine add_xyz_method
2473 
2474 !> Attaches "cell_methods" attribute to a variable based on defaults for axes_grp or optional arguments.
2475 subroutine attach_cell_methods(id, axes, ostring, cell_methods, &
2476  x_cell_method, y_cell_method, v_cell_method, v_extensive)
2477  integer, intent(in) :: id !< Handle to diagnostic
2478  type(axes_grp), intent(in) :: axes !< Container w/ up to 3 integer handles that indicates
2479  !! axes for this field
2480  character(len=*), intent(out) :: ostring !< The cell_methods strings that would appear in the file
2481  character(len=*), optional, intent(in) :: cell_methods !< String to append as cell_methods attribute.
2482  !! Use '' to have no attribute. If present, this
2483  !! overrides the default constructed from the default
2484  !! for each individual axis direction.
2485  character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction.
2486  !! Use '' have no method.
2487  character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction.
2488  !! Use '' have no method.
2489  character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction.
2490  !! Use '' have no method.
2491  logical, optional, intent(in) :: v_extensive !< True for vertically extensive fields
2492  !! (vertically integrated). Default/absent for intensive.
2493  ! Local variables
2494  character(len=9) :: axis_name
2495  logical :: x_mean, y_mean, x_sum, y_sum
2496 
2497  x_mean = .false.
2498  y_mean = .false.
2499  x_sum = .false.
2500  y_sum = .false.
2501 
2502  ostring = ''
2503  if (present(cell_methods)) then
2504  if (present(x_cell_method) .or. present(y_cell_method) .or. present(v_cell_method) &
2505  .or. present(v_extensive)) then
2506  call mom_error(fatal, "attach_cell_methods: " // &
2507  'Individual direction cell method was specified along with a "cell_methods" string.')
2508  endif
2509  if (len(trim(cell_methods))>0) then
2510  call diag_field_add_attribute(id, 'cell_methods', trim(cell_methods))
2511  ostring = trim(cell_methods)
2512  endif
2513  else
2514  if (present(x_cell_method)) then
2515  if (len(trim(x_cell_method))>0) then
2516  call get_diag_axis_name(axes%handles(1), axis_name)
2517  call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(x_cell_method))
2518  ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(x_cell_method)
2519  if (trim(x_cell_method)=='mean') x_mean=.true.
2520  if (trim(x_cell_method)=='sum') x_sum=.true.
2521  endif
2522  else
2523  if (len(trim(axes%x_cell_method))>0) then
2524  call get_diag_axis_name(axes%handles(1), axis_name)
2525  call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%x_cell_method))
2526  ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%x_cell_method)
2527  if (trim(axes%x_cell_method)=='mean') x_mean=.true.
2528  if (trim(axes%x_cell_method)=='sum') x_sum=.true.
2529  endif
2530  endif
2531  if (present(y_cell_method)) then
2532  if (len(trim(y_cell_method))>0) then
2533  call get_diag_axis_name(axes%handles(2), axis_name)
2534  call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(y_cell_method))
2535  ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(y_cell_method)
2536  if (trim(y_cell_method)=='mean') y_mean=.true.
2537  if (trim(y_cell_method)=='sum') y_sum=.true.
2538  endif
2539  else
2540  if (len(trim(axes%y_cell_method))>0) then
2541  call get_diag_axis_name(axes%handles(2), axis_name)
2542  call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%y_cell_method))
2543  ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%y_cell_method)
2544  if (trim(axes%y_cell_method)=='mean') y_mean=.true.
2545  if (trim(axes%y_cell_method)=='sum') y_sum=.true.
2546  endif
2547  endif
2548  if (present(v_cell_method)) then
2549  if (present(v_extensive)) call mom_error(fatal, "attach_cell_methods: " // &
2550  'Vertical cell method was specified along with the vertically extensive flag.')
2551  if (len(trim(v_cell_method))>0) then
2552  if (axes%rank==1) then
2553  call get_diag_axis_name(axes%handles(1), axis_name)
2554  elseif (axes%rank==3) then
2555  call get_diag_axis_name(axes%handles(3), axis_name)
2556  endif
2557  call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(v_cell_method))
2558  ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(v_cell_method)
2559  endif
2560  elseif (present(v_extensive)) then
2561  if(v_extensive) then
2562  if (axes%rank==1) then
2563  call get_diag_axis_name(axes%handles(1), axis_name)
2564  elseif (axes%rank==3) then
2565  call get_diag_axis_name(axes%handles(3), axis_name)
2566  endif
2567  call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':sum')
2568  ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':sum'
2569  endif
2570  else
2571  if (len(trim(axes%v_cell_method))>0) then
2572  if (axes%rank==1) then
2573  call get_diag_axis_name(axes%handles(1), axis_name)
2574  elseif (axes%rank==3) then
2575  call get_diag_axis_name(axes%handles(3), axis_name)
2576  endif
2577  call diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%v_cell_method))
2578  ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%v_cell_method)
2579  endif
2580  endif
2581  if (x_mean .and. y_mean) then
2582  call diag_field_add_attribute(id, 'cell_methods', 'area:mean')
2583  ostring = trim(adjustl(ostring))//' area:mean'
2584  elseif (x_sum .and. y_sum) then
2585  call diag_field_add_attribute(id, 'cell_methods', 'area:sum')
2586  ostring = trim(adjustl(ostring))//' area:sum'
2587  endif
2588  endif
2589  ostring = adjustl(ostring)
2590 end subroutine attach_cell_methods
2591 
2592 function register_scalar_field(module_name, field_name, init_time, diag_cs, &
2593  long_name, units, missing_value, range, standard_name, &
2594  do_not_log, err_msg, interp_method, cmor_field_name, &
2595  cmor_long_name, cmor_units, cmor_standard_name)
2596  integer :: register_scalar_field !< An integer handle for a diagnostic array.
2597  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model"
2598  !! or "ice_shelf_model"
2599  character(len=*), intent(in) :: field_name !< Name of the diagnostic field
2600  type(time_type), intent(in) :: init_time !< Time at which a field is first available?
2601  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output
2602  character(len=*), optional, intent(in) :: long_name !< Long name of a field.
2603  character(len=*), optional, intent(in) :: units !< Units of a field.
2604  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field
2605  real, optional, intent(in) :: missing_value !< A value that indicates missing values.
2606  real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)
2607  logical, optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?)
2608  character(len=*), optional, intent(out):: err_msg !< String into which an error message might be
2609  !! placed (not used in MOM?)
2610  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not
2611  !! be interpolated as a scalar
2612  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field
2613  character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field
2614  character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field
2615  character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field
2616 
2617  ! Local variables
2618  real :: mom_missing_value
2619  integer :: dm_id, fms_id
2620  type(diag_type), pointer :: diag => null(), cmor_diag => null()
2621  character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name
2622 
2623  mom_missing_value = diag_cs%missing_value
2624  if (present(missing_value)) mom_missing_value = missing_value
2625 
2626  dm_id = -1
2627  diag => null()
2628  cmor_diag => null()
2629 
2630  if (diag_cs%diag_as_chksum) then
2631  fms_id = diag_cs%num_chksum_diags + 1
2632  diag_cs%num_chksum_diags = fms_id
2633  else
2634  fms_id = register_diag_field_fms(module_name, field_name, init_time, &
2635  long_name=long_name, units=units, missing_value=mom_missing_value, &
2636  range=range, standard_name=standard_name, do_not_log=do_not_log, &
2637  err_msg=err_msg)
2638  endif
2639 
2640  if (fms_id /= diag_field_not_found) then
2641  dm_id = get_new_diag_id(diag_cs)
2642  call alloc_diag_with_id(dm_id, diag_cs, diag)
2643  call assert(associated(diag), 'register_scalar_field: diag allocation failed')
2644  diag%fms_diag_id = fms_id
2645  diag%debug_str = trim(module_name)//"-"//trim(field_name)
2646  endif
2647 
2648  if (present(cmor_field_name)) then
2649  ! Fallback values for strings set to "not provided"
2650  posted_cmor_units = "not provided"
2651  posted_cmor_standard_name = "not provided"
2652  posted_cmor_long_name = "not provided"
2653 
2654  ! If attributes are present for MOM variable names, use them first for the register_static_field
2655  ! call for CMOR verison of the variable
2656  if (present(units)) posted_cmor_units = units
2657  if (present(standard_name)) posted_cmor_standard_name = standard_name
2658  if (present(long_name)) posted_cmor_long_name = long_name
2659 
2660  ! If specified in the call to register_static_field, override attributes with the CMOR versions
2661  if (present(cmor_units)) posted_cmor_units = cmor_units
2662  if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name
2663  if (present(cmor_long_name)) posted_cmor_long_name = cmor_long_name
2664 
2665  fms_id = register_diag_field_fms(module_name, cmor_field_name, init_time, &
2666  long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), &
2667  missing_value=mom_missing_value, range=range, &
2668  standard_name=trim(posted_cmor_standard_name), do_not_log=do_not_log, err_msg=err_msg)
2669  if (fms_id /= diag_field_not_found) then
2670  if (dm_id == -1) then
2671  dm_id = get_new_diag_id(diag_cs)
2672  endif
2673  call alloc_diag_with_id(dm_id, diag_cs, cmor_diag)
2674  cmor_diag%fms_diag_id = fms_id
2675  cmor_diag%debug_str = trim(module_name)//"-"//trim(cmor_field_name)
2676  endif
2677  endif
2678 
2679  ! Document diagnostics in list of available diagnostics
2680  if (is_root_pe() .and. diag_cs%available_diag_doc_unit > 0) then
2681  call log_available_diag(associated(diag), module_name, field_name, '', '', diag_cs, &
2682  long_name, units, standard_name)
2683  if (present(cmor_field_name)) then
2684  call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, &
2685  '', '', diag_cs, posted_cmor_long_name, posted_cmor_units, &
2686  posted_cmor_standard_name)
2687  endif
2688  endif
2689 
2690  register_scalar_field = dm_id
2691 
2692 end function register_scalar_field
2693 
2694 !> Registers a static diagnostic, returning an integer handle
2695 function register_static_field(module_name, field_name, axes, &
2696  long_name, units, missing_value, range, mask_variant, standard_name, &
2697  do_not_log, interp_method, tile_count, &
2698  cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, area, &
2699  x_cell_method, y_cell_method, area_cell_method, conversion)
2700  integer :: register_static_field !< An integer handle for a diagnostic array.
2701  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model"
2702  !! or "ice_shelf_model"
2703  character(len=*), intent(in) :: field_name !< Name of the diagnostic field
2704  type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that
2705  !! indicates axes for this field
2706  character(len=*), optional, intent(in) :: long_name !< Long name of a field.
2707  character(len=*), optional, intent(in) :: units !< Units of a field.
2708  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field
2709  real, optional, intent(in) :: missing_value !< A value that indicates missing values.
2710  real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)
2711  logical, optional, intent(in) :: mask_variant !< If true a logical mask must be provided with
2712  !! post_data calls (not used in MOM?)
2713  logical, optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?)
2714  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not
2715  !! be interpolated as a scalar
2716  integer, optional, intent(in) :: tile_count !< no clue (not used in MOM?)
2717  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field
2718  character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field
2719  character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field
2720  character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field
2721  integer, optional, intent(in) :: area !< fms_id for area_t
2722  character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction.
2723  character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction.
2724  character(len=*), optional, intent(in) :: area_cell_method !< Specifies the cell method for area
2725  real, optional, intent(in) :: conversion !< A value to multiply data by before writing to file
2726 
2727  ! Local variables
2728  real :: mom_missing_value
2729  type(diag_ctrl), pointer :: diag_cs => null()
2730  type(diag_type), pointer :: diag => null(), cmor_diag => null()
2731  integer :: dm_id, fms_id, cmor_id
2732  character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name
2733  character(len=9) :: axis_name
2734 
2735  mom_missing_value = axes%diag_cs%missing_value
2736  if (present(missing_value)) mom_missing_value = missing_value
2737 
2738  diag_cs => axes%diag_cs
2739  dm_id = -1
2740  diag => null()
2741  cmor_diag => null()
2742 
2743  if (diag_cs%diag_as_chksum) then
2744  fms_id = diag_cs%num_chksum_diags + 1
2745  diag_cs%num_chksum_diags = fms_id
2746  else
2747  fms_id = register_static_field_fms(module_name, field_name, axes%handles, &
2748  long_name=long_name, units=units, missing_value=mom_missing_value, &
2749  range=range, mask_variant=mask_variant, standard_name=standard_name, &
2750  do_not_log=do_not_log, &
2751  interp_method=interp_method, tile_count=tile_count, area=area)
2752  endif
2753 
2754  if (fms_id /= diag_field_not_found) then
2755  dm_id = get_new_diag_id(diag_cs)
2756  call alloc_diag_with_id(dm_id, diag_cs, diag)
2757  call assert(associated(diag), 'register_static_field: diag allocation failed')
2758  diag%fms_diag_id = fms_id
2759  diag%debug_str = trim(module_name)//"-"//trim(field_name)
2760  if (present(conversion)) diag%conversion_factor = conversion
2761 
2762  if (diag_cs%diag_as_chksum) then
2763  diag%axes => axes
2764  else
2765  if (present(x_cell_method)) then
2766  call get_diag_axis_name(axes%handles(1), axis_name)
2767  call diag_field_add_attribute(fms_id, 'cell_methods', &
2768  trim(axis_name)//':'//trim(x_cell_method))
2769  endif
2770  if (present(y_cell_method)) then
2771  call get_diag_axis_name(axes%handles(2), axis_name)
2772  call diag_field_add_attribute(fms_id, 'cell_methods', &
2773  trim(axis_name)//':'//trim(y_cell_method))
2774  endif
2775  if (present(area_cell_method)) then
2776  call diag_field_add_attribute(fms_id, 'cell_methods', &
2777  'area:'//trim(area_cell_method))
2778  endif
2779  endif
2780  endif
2781 
2782  if (present(cmor_field_name) .and. .not. diag_cs%diag_as_chksum) then
2783  ! Fallback values for strings set to "not provided"
2784  posted_cmor_units = "not provided"
2785  posted_cmor_standard_name = "not provided"
2786  posted_cmor_long_name = "not provided"
2787 
2788  ! If attributes are present for MOM variable names, use them first for the register_static_field
2789  ! call for CMOR verison of the variable
2790  if (present(units)) posted_cmor_units = units
2791  if (present(standard_name)) posted_cmor_standard_name = standard_name
2792  if (present(long_name)) posted_cmor_long_name = long_name
2793 
2794  ! If specified in the call to register_static_field, override attributes with the CMOR versions
2795  if (present(cmor_units)) posted_cmor_units = cmor_units
2796  if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name
2797  if (present(cmor_long_name)) posted_cmor_long_name = cmor_long_name
2798 
2799  fms_id = register_static_field_fms(module_name, cmor_field_name, &
2800  axes%handles, long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), &
2801  missing_value=mom_missing_value, range=range, mask_variant=mask_variant, &
2802  standard_name=trim(posted_cmor_standard_name), do_not_log=do_not_log, &
2803  interp_method=interp_method, tile_count=tile_count, area=area)
2804  if (fms_id /= diag_field_not_found) then
2805  if (dm_id == -1) then
2806  dm_id = get_new_diag_id(diag_cs)
2807  endif
2808  call alloc_diag_with_id(dm_id, diag_cs, cmor_diag)
2809  cmor_diag%fms_diag_id = fms_id
2810  cmor_diag%debug_str = trim(module_name)//"-"//trim(cmor_field_name)
2811  if (present(conversion)) cmor_diag%conversion_factor = conversion
2812  if (present(x_cell_method)) then
2813  call get_diag_axis_name(axes%handles(1), axis_name)
2814  call diag_field_add_attribute(fms_id, 'cell_methods', trim(axis_name)//':'//trim(x_cell_method))
2815  endif
2816  if (present(y_cell_method)) then
2817  call get_diag_axis_name(axes%handles(2), axis_name)
2818  call diag_field_add_attribute(fms_id, 'cell_methods', trim(axis_name)//':'//trim(y_cell_method))
2819  endif
2820  if (present(area_cell_method)) then
2821  call diag_field_add_attribute(fms_id, 'cell_methods', 'area:'//trim(area_cell_method))
2822  endif
2823  endif
2824  endif
2825 
2826  ! Document diagnostics in list of available diagnostics
2827  if (is_root_pe() .and. diag_cs%available_diag_doc_unit > 0) then
2828  call log_available_diag(associated(diag), module_name, field_name, '', '', diag_cs, &
2829  long_name, units, standard_name)
2830  if (present(cmor_field_name)) then
2831  call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, &
2832  '', '', diag_cs, posted_cmor_long_name, posted_cmor_units, &
2833  posted_cmor_standard_name)
2834  endif
2835  endif
2836 
2837  register_static_field = dm_id
2838 
2839 end function register_static_field
2840 
2841 !> Describe an option setting in the diagnostic files.
2842 subroutine describe_option(opt_name, value, diag_CS)
2843  character(len=*), intent(in) :: opt_name !< The name of the option
2844  character(len=*), intent(in) :: value !< A character string with the setting of the option.
2845  type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output
2846 
2847  character(len=240) :: mesg
2848  integer :: len_ind
2849 
2850  len_ind = len_trim(value) ! Add error handling for long values?
2851 
2852  mesg = " ! "//trim(opt_name)//": "//trim(value)
2853  write(diag_cs%available_diag_doc_unit, '(a)') trim(mesg)
2854 end subroutine describe_option
2855 
2856 !> Registers a diagnostic using the information encapsulated in the vardesc
2857 !! type argument and returns an integer handle to this diagostic. That
2858 !! integer handle is negative if the diagnostic is unused.
2859 function ocean_register_diag(var_desc, G, diag_CS, day)
2860  integer :: ocean_register_diag !< An integer handle to this diagnostic.
2861  type(vardesc), intent(in) :: var_desc !< The vardesc type describing the diagnostic
2862  type(ocean_grid_type), intent(in) :: g !< The ocean's grid type
2863  type(diag_ctrl), intent(in), target :: diag_cs !< The diagnotic control structure
2864  type(time_type), intent(in) :: day !< The current model time
2865 
2866  character(len=64) :: var_name ! A variable's name.
2867  character(len=48) :: units ! A variable's units.
2868  character(len=240) :: longname ! A variable's longname.
2869  character(len=8) :: hor_grid, z_grid ! Variable grid info.
2870  type(axes_grp), pointer :: axes => null()
2871 
2872  call query_vardesc(var_desc, units=units, longname=longname, hor_grid=hor_grid, &
2873  z_grid=z_grid, caller="ocean_register_diag")
2874 
2875  ! Use the hor_grid and z_grid components of vardesc to determine the
2876  ! desired axes to register the diagnostic field for.
2877  select case (z_grid)
2878 
2879  case ("L")
2880  select case (hor_grid)
2881  case ("q")
2882  axes => diag_cs%axesBL
2883  case ("h")
2884  axes => diag_cs%axesTL
2885  case ("u")
2886  axes => diag_cs%axesCuL
2887  case ("v")
2888  axes => diag_cs%axesCvL
2889  case ("Bu")
2890  axes => diag_cs%axesBL
2891  case ("T")
2892  axes => diag_cs%axesTL
2893  case ("Cu")
2894  axes => diag_cs%axesCuL
2895  case ("Cv")
2896  axes => diag_cs%axesCvL
2897  case ("z")
2898  axes => diag_cs%axeszL
2899  case default
2900  call mom_error(fatal, "ocean_register_diag: " // &
2901  "unknown hor_grid component "//trim(hor_grid))
2902  end select
2903 
2904  case ("i")
2905  select case (hor_grid)
2906  case ("q")
2907  axes => diag_cs%axesBi
2908  case ("h")
2909  axes => diag_cs%axesTi
2910  case ("u")
2911  axes => diag_cs%axesCui
2912  case ("v")
2913  axes => diag_cs%axesCvi
2914  case ("Bu")
2915  axes => diag_cs%axesBi
2916  case ("T")
2917  axes => diag_cs%axesTi
2918  case ("Cu")
2919  axes => diag_cs%axesCui
2920  case ("Cv")
2921  axes => diag_cs%axesCvi
2922  case ("z")
2923  axes => diag_cs%axeszi
2924  case default
2925  call mom_error(fatal, "ocean_register_diag: " // &
2926  "unknown hor_grid component "//trim(hor_grid))
2927  end select
2928 
2929  case ("1")
2930  select case (hor_grid)
2931  case ("q")
2932  axes => diag_cs%axesB1
2933  case ("h")
2934  axes => diag_cs%axesT1
2935  case ("u")
2936  axes => diag_cs%axesCu1
2937  case ("v")
2938  axes => diag_cs%axesCv1
2939  case ("Bu")
2940  axes => diag_cs%axesB1
2941  case ("T")
2942  axes => diag_cs%axesT1
2943  case ("Cu")
2944  axes => diag_cs%axesCu1
2945  case ("Cv")
2946  axes => diag_cs%axesCv1
2947  case default
2948  call mom_error(fatal, "ocean_register_diag: " // &
2949  "unknown hor_grid component "//trim(hor_grid))
2950  end select
2951 
2952  case default
2953  call mom_error(fatal,&
2954  "ocean_register_diag: unknown z_grid component "//trim(z_grid))
2955  end select
2956 
2957  ocean_register_diag = register_diag_field("ocean_model", trim(var_name), &
2958  axes, day, trim(longname), trim(units), missing_value=-1.0e+34)
2959 
2960 end function ocean_register_diag
2961 
2962 subroutine diag_mediator_infrastructure_init(err_msg)
2963  ! This subroutine initializes the FMS diag_manager.
2964  character(len=*), optional, intent(out) :: err_msg !< An error message
2965 
2966  call diag_manager_init(err_msg=err_msg)
2967 end subroutine diag_mediator_infrastructure_init
2968 
2969 !> diag_mediator_init initializes the MOM diag_mediator and opens the available
2970 !! diagnostics file, if appropriate.
2971 subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)
2972  type(ocean_grid_type), target, intent(inout) :: g !< The ocean grid type.
2973  type(verticalgrid_type), target, intent(in) :: gv !< The ocean vertical grid structure
2974  type(unit_scale_type), target, intent(in) :: us !< A dimensional unit scaling type
2975  integer, intent(in) :: nz !< The number of layers in the model's native grid.
2976  type(param_file_type), intent(in) :: param_file !< Parameter file structure
2977  type(diag_ctrl), intent(inout) :: diag_cs !< A pointer to a type with many variables
2978  !! used for diagnostics
2979  character(len=*), optional, intent(in) :: doc_file_dir !< A directory in which to create the
2980  !! file
2981 
2982  ! This subroutine initializes the diag_mediator and the diag_manager.
2983  ! The grid type should have its dimensions set by this point, but it
2984  ! is not necessary that the metrics and axis labels be set up yet.
2985 
2986  ! Local variables
2987  integer :: ios, i, new_unit
2988  logical :: opened, new_file
2989  character(len=8) :: this_pe
2990  character(len=240) :: doc_file, doc_file_dflt, doc_path
2991  character(len=240), allocatable :: diag_coords(:)
2992  ! This include declares and sets the variable "version".
2993 # include "version_variable.h"
2994  character(len=40) :: mdl = "MOM_diag_mediator" ! This module's name.
2995  character(len=32) :: filename_appendix = '' !fms appendix to filename for ensemble runs
2996 
2997  id_clock_diag_mediator = cpu_clock_id('(Ocean diagnostics framework)', grain=clock_module)
2998  id_clock_diag_remap = cpu_clock_id('(Ocean diagnostics remapping)', grain=clock_routine)
2999  id_clock_diag_grid_updates = cpu_clock_id('(Ocean diagnostics grid updates)', grain=clock_routine)
3000 
3001  ! Allocate and initialize list of all diagnostics (and variants)
3002  allocate(diag_cs%diags(diag_alloc_chunk_size))
3003  diag_cs%next_free_diag_id = 1
3004  do i=1, diag_alloc_chunk_size
3005  call initialize_diag_type(diag_cs%diags(i))
3006  enddo
3007 
3008  ! Read all relevant parameters and write them to the model log.
3009  call log_version(param_file, mdl, version, "")
3010 
3011  call get_param(param_file, mdl, 'NUM_DIAG_COORDS', diag_cs%num_diag_coords, &
3012  'The number of diagnostic vertical coordinates to use. '//&
3013  'For each coordinate, an entry in DIAG_COORDS must be provided.', &
3014  default=1)
3015  if (diag_cs%num_diag_coords>0) then
3016  allocate(diag_coords(diag_cs%num_diag_coords))
3017  if (diag_cs%num_diag_coords==1) then ! The default is to provide just one instance of Z*
3018  call get_param(param_file, mdl, 'DIAG_COORDS', diag_coords, &
3019  'A list of string tuples associating diag_table modules to '//&
3020  'a coordinate definition used for diagnostics. Each string '//&
3021  'is of the form "MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME".', &
3022  default='z Z ZSTAR')
3023  else ! If using more than 1 diagnostic coordinate, all must be explicitly defined
3024  call get_param(param_file, mdl, 'DIAG_COORDS', diag_coords, &
3025  'A list of string tuples associating diag_table modules to '//&
3026  'a coordinate definition used for diagnostics. Each string '//&
3027  'is of the form "MODULE_SUFFIX,PARAMETER_SUFFIX,COORDINATE_NAME".', &
3028  fail_if_missing=.true.)
3029  endif
3030  allocate(diag_cs%diag_remap_cs(diag_cs%num_diag_coords))
3031  ! Initialize each diagnostic vertical coordinate
3032  do i=1, diag_cs%num_diag_coords
3033  call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i))
3034  enddo
3035  deallocate(diag_coords)
3036  endif
3037 
3038  call get_param(param_file, mdl, 'DIAG_MISVAL', diag_cs%missing_value, &
3039  'Set the default missing value to use for diagnostics.', &
3040  default=1.e20)
3041  call get_param(param_file, mdl, 'DIAG_AS_CHKSUM', diag_cs%diag_as_chksum, &
3042  'Instead of writing diagnostics to the diag manager, write '//&
3043  'a text file containing the checksum (bitcount) of the array.', &
3044  default=.false.)
3045 
3046  if (diag_cs%diag_as_chksum) &
3047  diag_cs%num_chksum_diags = 0
3048 
3049  ! Keep pointers grid, h, T, S needed diagnostic remapping
3050  diag_cs%G => g
3051  diag_cs%GV => gv
3052  diag_cs%US => us
3053  diag_cs%h => null()
3054  diag_cs%T => null()
3055  diag_cs%S => null()
3056  diag_cs%eqn_of_state => null()
3057 
3058 #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)
3059  allocate(diag_cs%h_old(g%isd:g%ied,g%jsd:g%jed,nz))
3060  diag_cs%h_old(:,:,:) = 0.0
3061 #endif
3062 
3063  diag_cs%is = g%isc - (g%isd-1) ; diag_cs%ie = g%iec - (g%isd-1)
3064  diag_cs%js = g%jsc - (g%jsd-1) ; diag_cs%je = g%jec - (g%jsd-1)
3065  diag_cs%isd = g%isd ; diag_cs%ied = g%ied
3066  diag_cs%jsd = g%jsd ; diag_cs%jed = g%jed
3067 
3068  !Downsample indices for dl=2 (should be generalized to arbitrary dl, perhaps via a G array)
3069  diag_cs%dsamp(2)%isc = g%HId2%isc - (g%HId2%isd-1) ; diag_cs%dsamp(2)%iec = g%HId2%iec - (g%HId2%isd-1)
3070  diag_cs%dsamp(2)%jsc = g%HId2%jsc - (g%HId2%jsd-1) ; diag_cs%dsamp(2)%jec = g%HId2%jec - (g%HId2%jsd-1)
3071  diag_cs%dsamp(2)%isd = g%HId2%isd ; diag_cs%dsamp(2)%ied = g%HId2%ied
3072  diag_cs%dsamp(2)%jsd = g%HId2%jsd ; diag_cs%dsamp(2)%jed = g%HId2%jed
3073  diag_cs%dsamp(2)%isg = g%HId2%isg ; diag_cs%dsamp(2)%ieg = g%HId2%ieg
3074  diag_cs%dsamp(2)%jsg = g%HId2%jsg ; diag_cs%dsamp(2)%jeg = g%HId2%jeg
3075  diag_cs%dsamp(2)%isgB = g%HId2%isgB ; diag_cs%dsamp(2)%iegB = g%HId2%iegB
3076  diag_cs%dsamp(2)%jsgB = g%HId2%jsgB ; diag_cs%dsamp(2)%jegB = g%HId2%jegB
3077 
3078  ! Initialze available diagnostic log file
3079  if (is_root_pe() .and. (diag_cs%available_diag_doc_unit < 0)) then
3080  write(this_pe,'(i6.6)') pe_here()
3081  doc_file_dflt = "available_diags."//this_pe
3082  call get_param(param_file, mdl, "AVAILABLE_DIAGS_FILE", doc_file, &
3083  "A file into which to write a list of all available "//&
3084  "ocean diagnostics that can be included in a diag_table.", &
3085  default=doc_file_dflt, do_not_log=(diag_cs%available_diag_doc_unit/=-1))
3086  if (len_trim(doc_file) > 0) then
3087  new_file = .true. ; if (diag_cs%available_diag_doc_unit /= -1) new_file = .false.
3088  ! Find an unused unit number.
3089  do new_unit=512,42,-1
3090  inquire( new_unit, opened=opened)
3091  if (.not.opened) exit
3092  enddo
3093  if (opened) call mom_error(fatal, &
3094  "diag_mediator_init failed to find an unused unit number.")
3095 
3096  doc_path = doc_file
3097  if (present(doc_file_dir)) then ; if (len_trim(doc_file_dir) > 0) then
3098  doc_path = trim(slasher(doc_file_dir))//trim(doc_file)
3099  endif ; endif
3100 
3101  diag_cs%available_diag_doc_unit = new_unit
3102 
3103  if (new_file) then
3104  open(diag_cs%available_diag_doc_unit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', &
3105  action='WRITE', status='REPLACE', iostat=ios)
3106  else ! This file is being reopened, and should be appended.
3107  open(diag_cs%available_diag_doc_unit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', &
3108  action='WRITE', status='OLD', position='APPEND', iostat=ios)
3109  endif
3110  inquire(diag_cs%available_diag_doc_unit, opened=opened)
3111  if ((.not.opened) .or. (ios /= 0)) then
3112  call mom_error(fatal, "Failed to open available diags file "//trim(doc_path)//".")
3113  endif
3114  endif
3115  endif
3116 
3117  if (is_root_pe() .and. (diag_cs%chksum_iounit < 0) .and. diag_cs%diag_as_chksum) then
3118  !write(this_pe,'(i6.6)') PE_here()
3119  !doc_file_dflt = "chksum_diag."//this_pe
3120  doc_file_dflt = "chksum_diag"
3121  call get_param(param_file, mdl, "CHKSUM_DIAG_FILE", doc_file, &
3122  "A file into which to write all checksums of the "//&
3123  "diagnostics listed in the diag_table.", &
3124  default=doc_file_dflt, do_not_log=(diag_cs%chksum_iounit/=-1))
3125 
3126  call get_filename_appendix(filename_appendix)
3127  if (len_trim(filename_appendix) > 0) then
3128  doc_file = trim(doc_file) //'.'//trim(filename_appendix)
3129  endif
3130 #ifdef STATSLABEL
3131  doc_file = trim(doc_file)//"."//trim(adjustl(statslabel))
3132 #endif
3133 
3134  if (len_trim(doc_file) > 0) then
3135  new_file = .true. ; if (diag_cs%chksum_iounit /= -1) new_file = .false.
3136  ! Find an unused unit number.
3137  do new_unit=512,42,-1
3138  inquire( new_unit, opened=opened)
3139  if (.not.opened) exit
3140  enddo
3141  if (opened) call mom_error(fatal, &
3142  "diag_mediator_init failed to find an unused unit number.")
3143 
3144  doc_path = doc_file
3145  if (present(doc_file_dir)) then ; if (len_trim(doc_file_dir) > 0) then
3146  doc_path = trim(slasher(doc_file_dir))//trim(doc_file)
3147  endif ; endif
3148 
3149  diag_cs%chksum_iounit = new_unit
3150 
3151  if (new_file) then
3152  open(diag_cs%chksum_iounit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', &
3153  action='WRITE', status='REPLACE', iostat=ios)
3154  else ! This file is being reopened, and should be appended.
3155  open(diag_cs%chksum_iounit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', &
3156  action='WRITE', status='OLD', position='APPEND', iostat=ios)
3157  endif
3158  inquire(diag_cs%chksum_iounit, opened=opened)
3159  if ((.not.opened) .or. (ios /= 0)) then
3160  call mom_error(fatal, "Failed to open checksum diags file "//trim(doc_path)//".")
3161  endif
3162  endif
3163  endif
3164 
3165 end subroutine diag_mediator_init
3166 
3167 !> Set pointers to the default state fields used to remap diagnostics.
3168 subroutine diag_set_state_ptrs(h, T, S, eqn_of_state, diag_cs)
3169  real, dimension(:,:,:), target, intent(in ) :: h !< the model thickness array [H ~> m or kg m-2]
3170  real, dimension(:,:,:), target, intent(in ) :: t !< the model temperature array
3171  real, dimension(:,:,:), target, intent(in ) :: s !< the model salinity array
3172  type(eos_type), target, intent(in ) :: eqn_of_state !< Equation of state structure
3173  type(diag_ctrl), intent(inout) :: diag_cs !< diag mediator control structure
3174 
3175  ! Keep pointers to h, T, S needed for the diagnostic remapping
3176  diag_cs%h => h
3177  diag_cs%T => t
3178  diag_cs%S => s
3179  diag_cs%eqn_of_state => eqn_of_state
3180 
3181 end subroutine
3182 
3183 !> Build/update vertical grids for diagnostic remapping.
3184 !! \note The target grids need to be updated whenever sea surface
3185 !! height changes.
3186 subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S)
3187  type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
3188  real, target, optional, intent(in ) :: alt_h(:,:,:) !< Used if remapped grids should be something other than
3189  !! the current thicknesses [H ~> m or kg m-2]
3190  real, target, optional, intent(in ) :: alt_t(:,:,:) !< Used if remapped grids should be something other than
3191  !! the current temperatures
3192  real, target, optional, intent(in ) :: alt_s(:,:,:) !< Used if remapped grids should be something other than
3193  !! the current salinity
3194  ! Local variables
3195  integer :: i
3196  real, dimension(:,:,:), pointer :: h_diag => null()
3197  real, dimension(:,:,:), pointer :: t_diag => null(), s_diag => null()
3198 
3199  if (present(alt_h)) then
3200  h_diag => alt_h
3201  else
3202  h_diag => diag_cs%h
3203  endif
3204 
3205  if (present(alt_t)) then
3206  t_diag => alt_t
3207  else
3208  t_diag => diag_cs%T
3209  endif
3210 
3211  if (present(alt_s)) then
3212  s_diag => alt_s
3213  else
3214  s_diag => diag_cs%S
3215  endif
3216 
3217  if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates)
3218 
3219  if (diag_cs%diag_grid_overridden) then
3220  call mom_error(fatal, "diag_update_remap_grids was called, but current grids in "// &
3221  "diagnostic structure have been overridden")
3222  endif
3223 
3224  do i=1, diag_cs%num_diag_coords
3225  call diag_remap_update(diag_cs%diag_remap_cs(i), &
3226  diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, t_diag, s_diag, &
3227  diag_cs%eqn_of_state)
3228  enddo
3229 
3230 #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)
3231  ! Keep a copy of H - used to check whether grids are up-to-date
3232  ! when doing remapping.
3233  diag_cs%h_old(:,:,:) = diag_cs%h(:,:,:)
3234 #endif
3235 
3236  if (id_clock_diag_grid_updates>0) call cpu_clock_end(id_clock_diag_grid_updates)
3237 
3238 end subroutine diag_update_remap_grids
3239 
3240 !> Sets up the 2d and 3d masks for native diagnostics
3241 subroutine diag_masks_set(G, nz, diag_cs)
3242  type(ocean_grid_type), target, intent(in) :: g !< The ocean grid type.
3243  integer, intent(in) :: nz !< The number of layers in the model's native grid.
3244  type(diag_ctrl), pointer :: diag_cs !< A pointer to a type with many variables
3245  !! used for diagnostics
3246  ! Local variables
3247  integer :: k
3248 
3249  ! 2d masks point to the model masks since they are identical
3250  diag_cs%mask2dT => g%mask2dT
3251  diag_cs%mask2dBu => g%mask2dBu
3252  diag_cs%mask2dCu => g%mask2dCu
3253  diag_cs%mask2dCv => g%mask2dCv
3254 
3255  ! 3d native masks are needed by diag_manager but the native variables
3256  ! can only be masked 2d - for ocean points, all layers exists.
3257  allocate(diag_cs%mask3dTL(g%isd:g%ied,g%jsd:g%jed,1:nz))
3258  allocate(diag_cs%mask3dBL(g%IsdB:g%IedB,g%JsdB:g%JedB,1:nz))
3259  allocate(diag_cs%mask3dCuL(g%IsdB:g%IedB,g%jsd:g%jed,1:nz))
3260  allocate(diag_cs%mask3dCvL(g%isd:g%ied,g%JsdB:g%JedB,1:nz))
3261  do k=1,nz
3262  diag_cs%mask3dTL(:,:,k) = diag_cs%mask2dT(:,:)
3263  diag_cs%mask3dBL(:,:,k) = diag_cs%mask2dBu(:,:)
3264  diag_cs%mask3dCuL(:,:,k) = diag_cs%mask2dCu(:,:)
3265  diag_cs%mask3dCvL(:,:,k) = diag_cs%mask2dCv(:,:)
3266  enddo
3267  allocate(diag_cs%mask3dTi(g%isd:g%ied,g%jsd:g%jed,1:nz+1))
3268  allocate(diag_cs%mask3dBi(g%IsdB:g%IedB,g%JsdB:g%JedB,1:nz+1))
3269  allocate(diag_cs%mask3dCui(g%IsdB:g%IedB,g%jsd:g%jed,1:nz+1))
3270  allocate(diag_cs%mask3dCvi(g%isd:g%ied,g%JsdB:g%JedB,1:nz+1))
3271  do k=1,nz+1
3272  diag_cs%mask3dTi(:,:,k) = diag_cs%mask2dT(:,:)
3273  diag_cs%mask3dBi(:,:,k) = diag_cs%mask2dBu(:,:)
3274  diag_cs%mask3dCui(:,:,k) = diag_cs%mask2dCu(:,:)
3275  diag_cs%mask3dCvi(:,:,k) = diag_cs%mask2dCv(:,:)
3276  enddo
3277 
3278  !Allocate and initialize the downsampled masks
3279  call downsample_diag_masks_set(g, nz, diag_cs)
3280 
3281 end subroutine diag_masks_set
3282 
3283 subroutine diag_mediator_close_registration(diag_CS)
3284  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output
3285 
3286  integer :: i
3287 
3288  if (diag_cs%available_diag_doc_unit > -1) then
3289  close(diag_cs%available_diag_doc_unit) ; diag_cs%available_diag_doc_unit = -2
3290  endif
3291 
3292  do i=1, diag_cs%num_diag_coords
3293  call diag_remap_diag_registration_closed(diag_cs%diag_remap_cs(i))
3294  enddo
3295 
3296 end subroutine diag_mediator_close_registration
3297 
3298 subroutine diag_mediator_end(time, diag_CS, end_diag_manager)
3299  type(time_type), intent(in) :: time !< The current model time
3300  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output
3301  logical, optional, intent(in) :: end_diag_manager !< If true, call diag_manager_end()
3302 
3303  ! Local variables
3304  integer :: i
3305 
3306  if (diag_cs%available_diag_doc_unit > -1) then
3307  close(diag_cs%available_diag_doc_unit) ; diag_cs%available_diag_doc_unit = -3
3308  endif
3309  if (diag_cs%chksum_iounit > -1) then
3310  close(diag_cs%chksum_iounit) ; diag_cs%chksum_iounit = -3
3311  endif
3312 
3313  deallocate(diag_cs%diags)
3314 
3315  do i=1, diag_cs%num_diag_coords
3316  call diag_remap_end(diag_cs%diag_remap_cs(i))
3317  enddo
3318 
3319  call diag_grid_storage_end(diag_cs%diag_grid_temp)
3320  deallocate(diag_cs%mask3dTL)
3321  deallocate(diag_cs%mask3dBL)
3322  deallocate(diag_cs%mask3dCuL)
3323  deallocate(diag_cs%mask3dCvL)
3324  deallocate(diag_cs%mask3dTi)
3325  deallocate(diag_cs%mask3dBi)
3326  deallocate(diag_cs%mask3dCui)
3327  deallocate(diag_cs%mask3dCvi)
3328  do i=2,max_dsamp_lev
3329  deallocate(diag_cs%dsamp(i)%mask2dT)
3330  deallocate(diag_cs%dsamp(i)%mask2dBu)
3331  deallocate(diag_cs%dsamp(i)%mask2dCu)
3332  deallocate(diag_cs%dsamp(i)%mask2dCv)
3333  deallocate(diag_cs%dsamp(i)%mask3dTL)
3334  deallocate(diag_cs%dsamp(i)%mask3dBL)
3335  deallocate(diag_cs%dsamp(i)%mask3dCuL)
3336  deallocate(diag_cs%dsamp(i)%mask3dCvL)
3337  deallocate(diag_cs%dsamp(i)%mask3dTi)
3338  deallocate(diag_cs%dsamp(i)%mask3dBi)
3339  deallocate(diag_cs%dsamp(i)%mask3dCui)
3340  deallocate(diag_cs%dsamp(i)%mask3dCvi)
3341  enddo
3342 
3343 #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)
3344  deallocate(diag_cs%h_old)
3345 #endif
3346 
3347  if (present(end_diag_manager)) then
3348  if (end_diag_manager) call diag_manager_end(time)
3349  endif
3350 
3351 end subroutine diag_mediator_end
3352 
3353 !> Convert the first n elements (up to 3) of an integer array to an underscore delimited string.
3354 function i2s(a,n_in)
3355  ! "Convert the first n elements of an integer array to a string."
3356  ! Perhaps this belongs elsewhere in the MOM6 code?
3357  integer, dimension(:), intent(in) :: a !< The array of integers to translate
3358  integer, optional , intent(in) :: n_in !< The number of elements to translate, by default all
3359  character(len=15) :: i2s !< The returned string
3360 
3361  character(len=15) :: i2s_temp
3362  integer :: i,n
3363 
3364  n=size(a)
3365  if (present(n_in)) n = n_in
3366 
3367  i2s = ''
3368  do i=1,min(n,3)
3369  write (i2s_temp, '(I4.4)') a(i)
3370  i2s = trim(i2s) //'_'// trim(i2s_temp)
3371  enddo
3372  i2s = adjustl(i2s)
3373 end function i2s
3374 
3375 !> Returns a new diagnostic id, it may be necessary to expand the diagnostics array.
3376 integer function get_new_diag_id(diag_cs)
3377  type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
3378  ! Local variables
3379  type(diag_type), dimension(:), allocatable :: tmp
3380  integer :: i
3381 
3382  if (diag_cs%next_free_diag_id > size(diag_cs%diags)) then
3383  call assert(diag_cs%next_free_diag_id - size(diag_cs%diags) == 1, &
3384  'get_new_diag_id: inconsistent diag id')
3385 
3386  ! Increase the size of diag_cs%diags and copy data over.
3387  ! Do not use move_alloc() because it is not supported by Fortran 90
3388  allocate(tmp(size(diag_cs%diags)))
3389  tmp(:) = diag_cs%diags(:)
3390  deallocate(diag_cs%diags)
3391  allocate(diag_cs%diags(size(tmp) + diag_alloc_chunk_size))
3392  diag_cs%diags(1:size(tmp)) = tmp(:)
3393  deallocate(tmp)
3394 
3395  ! Initialize new part of the diag array.
3396  do i=diag_cs%next_free_diag_id, size(diag_cs%diags)
3397  call initialize_diag_type(diag_cs%diags(i))
3398  enddo
3399  endif
3400 
3401  get_new_diag_id = diag_cs%next_free_diag_id
3402  diag_cs%next_free_diag_id = diag_cs%next_free_diag_id + 1
3403 
3404 end function get_new_diag_id
3405 
3406 !> Initializes a diag_type (used after allocating new memory)
3407 subroutine initialize_diag_type(diag)
3408  type(diag_type), intent(inout) :: diag !< diag_type to be initialized
3409 
3410  diag%in_use = .false.
3411  diag%fms_diag_id = -1
3412  diag%axes => null()
3413  diag%next => null()
3414  diag%conversion_factor = 0.
3415 
3416 end subroutine initialize_diag_type
3417 
3418 !> Make a new diagnostic. Either use memory which is in the array of 'primary'
3419 !! diagnostics, or if that is in use, insert it to the list of secondary diags.
3420 subroutine alloc_diag_with_id(diag_id, diag_cs, diag)
3421  integer, intent(in ) :: diag_id !< id for the diagnostic
3422  type(diag_ctrl), target, intent(inout) :: diag_cs !< structure used to regulate diagnostic output
3423  type(diag_type), pointer :: diag !< structure representing a diagnostic (inout)
3424 
3425  type(diag_type), pointer :: tmp => null()
3426 
3427  if (.not. diag_cs%diags(diag_id)%in_use) then
3428  diag => diag_cs%diags(diag_id)
3429  else
3430  allocate(diag)
3431  tmp => diag_cs%diags(diag_id)%next
3432  diag_cs%diags(diag_id)%next => diag
3433  diag%next => tmp
3434  endif
3435  diag%in_use = .true.
3436 
3437 end subroutine alloc_diag_with_id
3438 
3439 !> Log a diagnostic to the available diagnostics file.
3440 subroutine log_available_diag(used, module_name, field_name, cell_methods_string, comment, &
3441  diag_CS, long_name, units, standard_name)
3442  logical, intent(in) :: used !< Whether this diagnostic was in the diag_table or not
3443  character(len=*), intent(in) :: module_name !< Name of the diagnostic module
3444  character(len=*), intent(in) :: field_name !< Name of this diagnostic field
3445  character(len=*), intent(in) :: cell_methods_string !< The spatial component of the CF cell_methods attribute
3446  character(len=*), intent(in) :: comment !< A comment to append after [Used|Unused]
3447  type(diag_ctrl), intent(in) :: diag_CS !< The diagnotics control structure
3448  character(len=*), optional, intent(in) :: long_name !< CF long name of diagnostic
3449  character(len=*), optional, intent(in) :: units !< Units for diagnostic
3450  character(len=*), optional, intent(in) :: standard_name !< CF standardized name of diagnostic
3451  ! Local variables
3452  character(len=240) :: mesg
3453 
3454  if (used) then
3455  mesg = '"'//trim(module_name)//'", "'//trim(field_name)//'" [Used]'
3456  else
3457  mesg = '"'//trim(module_name)//'", "'//trim(field_name)//'" [Unused]'
3458  endif
3459  if (len(trim((comment)))>0) then
3460  write(diag_cs%available_diag_doc_unit, '(a,x,"(",a,")")') trim(mesg),trim(comment)
3461  else
3462  write(diag_cs%available_diag_doc_unit, '(a)') trim(mesg)
3463  endif
3464  if (present(long_name)) call describe_option("long_name", long_name, diag_cs)
3465  if (present(units)) call describe_option("units", units, diag_cs)
3466  if (present(standard_name)) &
3467  call describe_option("standard_name", standard_name, diag_cs)
3468  if (len(trim((cell_methods_string)))>0) &
3469  call describe_option("cell_methods", trim(cell_methods_string), diag_cs)
3470 
3471 end subroutine log_available_diag
3472 
3473 !> Log the diagnostic chksum to the chksum diag file
3474 subroutine log_chksum_diag(docunit, description, chksum)
3475  integer, intent(in) :: docunit !< Handle of the log file
3476  character(len=*), intent(in) :: description !< Name of the diagnostic module
3477  integer, intent(in) :: chksum !< chksum of the diagnostic
3478 
3479  write(docunit, '(a,x,i9.8)') description, chksum
3480  flush(docunit)
3481 
3482 end subroutine log_chksum_diag
3483 
3484 !> Allocates fields necessary to store diagnostic remapping fields
3485 subroutine diag_grid_storage_init(grid_storage, G, diag)
3486  type(diag_grid_storage), intent(inout) :: grid_storage !< Structure containing a snapshot of the target grids
3487  type(ocean_grid_type), intent(in) :: g !< Horizontal grid
3488  type(diag_ctrl), intent(in) :: diag !< Diagnostic control structure used as the contructor
3489  !! template for this routine
3490 
3491  integer :: m, nz
3492  grid_storage%num_diag_coords = diag%num_diag_coords
3493 
3494  ! Don't do anything else if there are no remapped coordinates
3495  if (grid_storage%num_diag_coords < 1) return
3496 
3497  ! Allocate memory for the native space
3498  allocate(grid_storage%h_state(g%isd:g%ied,g%jsd:g%jed, g%ke))
3499  ! Allocate diagnostic remapping structures
3500  allocate(grid_storage%diag_grids(diag%num_diag_coords))
3501  ! Loop through and allocate memory for the grid on each target coordinate
3502  do m = 1, diag%num_diag_coords
3503  nz = diag%diag_remap_cs(m)%nz
3504  allocate(grid_storage%diag_grids(m)%h(g%isd:g%ied,g%jsd:g%jed, nz))
3505  enddo
3506 
3507 end subroutine diag_grid_storage_init
3508 
3509 !> Copy from the main diagnostic arrays to the grid storage as well as the native thicknesses
3510 subroutine diag_copy_diag_to_storage(grid_storage, h_state, diag)
3511  type(diag_grid_storage), intent(inout) :: grid_storage !< Structure containing a snapshot of the target grids
3512  real, dimension(:,:,:), intent(in) :: h_state !< Current model thicknesses
3513  type(diag_ctrl), intent(in) :: diag !< Diagnostic control structure used as the contructor
3514 
3515  integer :: m
3516 
3517  ! Don't do anything else if there are no remapped coordinates
3518  if (grid_storage%num_diag_coords < 1) return
3519 
3520  grid_storage%h_state(:,:,:) = h_state(:,:,:)
3521  do m = 1,grid_storage%num_diag_coords
3522  if (diag%diag_remap_cs(m)%nz > 0) &
3523  grid_storage%diag_grids(m)%h(:,:,:) = diag%diag_remap_cs(m)%h(:,:,:)
3524  enddo
3525 
3526 end subroutine diag_copy_diag_to_storage
3527 
3528 !> Copy from the stored diagnostic arrays to the main diagnostic grids
3529 subroutine diag_copy_storage_to_diag(diag, grid_storage)
3530  type(diag_ctrl), intent(inout) :: diag !< Diagnostic control structure used as the contructor
3531  type(diag_grid_storage), intent(in) :: grid_storage !< Structure containing a snapshot of the target grids
3532 
3533  integer :: m
3534 
3535  ! Don't do anything else if there are no remapped coordinates
3536  if (grid_storage%num_diag_coords < 1) return
3537 
3538  diag%diag_grid_overridden = .true.
3539  do m = 1,grid_storage%num_diag_coords
3540  if (diag%diag_remap_cs(m)%nz > 0) &
3541  diag%diag_remap_cs(m)%h(:,:,:) = grid_storage%diag_grids(m)%h(:,:,:)
3542  enddo
3543 
3544 end subroutine diag_copy_storage_to_diag
3545 
3546 !> Save the current diagnostic grids in the temporary structure within diag
3547 subroutine diag_save_grids(diag)
3548  type(diag_ctrl), intent(inout) :: diag !< Diagnostic control structure used as the contructor
3549 
3550  integer :: m
3551 
3552  ! Don't do anything else if there are no remapped coordinates
3553  if (diag%num_diag_coords < 1) return
3554 
3555  do m = 1,diag%num_diag_coords
3556  if (diag%diag_remap_cs(m)%nz > 0) &
3557  diag%diag_grid_temp%diag_grids(m)%h(:,:,:) = diag%diag_remap_cs(m)%h(:,:,:)
3558  enddo
3559 
3560 end subroutine diag_save_grids
3561 
3562 !> Restore the diagnostic grids from the temporary structure within diag
3563 subroutine diag_restore_grids(diag)
3564  type(diag_ctrl), intent(inout) :: diag !< Diagnostic control structure used as the contructor
3565 
3566  integer :: m
3567 
3568  ! Don't do anything else if there are no remapped coordinates
3569  if (diag%num_diag_coords < 1) return
3570 
3571  diag%diag_grid_overridden = .false.
3572  do m = 1,diag%num_diag_coords
3573  if (diag%diag_remap_cs(m)%nz > 0) &
3574  diag%diag_remap_cs(m)%h(:,:,:) = diag%diag_grid_temp%diag_grids(m)%h(:,:,:)
3575  enddo
3576 
3577 end subroutine diag_restore_grids
3578 
3579 !> Deallocates the fields in the remapping fields container
3580 subroutine diag_grid_storage_end(grid_storage)
3581  type(diag_grid_storage), intent(inout) :: grid_storage !< Structure containing a snapshot of the target grids
3582  ! Local variables
3583  integer :: m, nz
3584 
3585  ! Don't do anything else if there are no remapped coordinates
3586  if (grid_storage%num_diag_coords < 1) return
3587 
3588  ! Deallocate memory for the native space
3589  deallocate(grid_storage%h_state)
3590  ! Loop through and deallocate memory for the grid on each target coordinate
3591  do m = 1, grid_storage%num_diag_coords
3592  deallocate(grid_storage%diag_grids(m)%h)
3593  enddo
3594  ! Deallocate diagnostic remapping structures
3595  deallocate(grid_storage%diag_grids)
3596 end subroutine diag_grid_storage_end
3597 
3598 !< Allocate and initialize the masks for downsampled diagostics in diag_cs
3599 !! The downsampled masks in the axes would later "point" to these.
3600 subroutine downsample_diag_masks_set(G, nz, diag_cs)
3601  type(ocean_grid_type), target, intent(in) :: G !< The ocean grid type.
3602  integer, intent(in) :: nz !< The number of layers in the model's native grid.
3603  type(diag_ctrl), pointer :: diag_cs !< A pointer to a type with many variables
3604  !! used for diagnostics
3605  ! Local variables
3606  integer :: i,j,k,ii,jj,dl
3607 
3608 !print*,'original c extents ',G%isc,G%iec,G%jsc,G%jec
3609 !print*,'original c extents ',G%iscb,G%iecb,G%jscb,G%jecb
3610 !print*,'coarse c extents ',G%HId2%isc,G%HId2%iec,G%HId2%jsc,G%HId2%jec
3611 !print*,'original d extents ',G%isd,G%ied,G%jsd,G%jed
3612 !print*,'coarse d extents ',G%HId2%isd,G%HId2%ied,G%HId2%jsd,G%HId2%jed
3613 ! original c extents 5 52 5 52
3614 ! original cB-nonsym extents 5 52 5 52
3615 ! original cB-sym extents 4 52 4 52
3616 ! coarse c extents 3 26 3 26
3617 ! original d extents 1 56 1 56
3618 ! original dB-nonsym extents 1 56 1 56
3619 ! original dB-sym extents 0 56 0 56
3620 ! coarse d extents 1 28 1 28
3621 
3622  do dl=2,max_dsamp_lev
3623  ! 2d mask
3624  call downsample_mask(g%mask2dT, diag_cs%dsamp(dl)%mask2dT, dl,g%isc, g%jsc, &
3625  g%HId2%isc, g%HId2%iec, g%HId2%jsc, g%HId2%jec, g%HId2%isd, g%HId2%ied, g%HId2%jsd, g%HId2%jed)
3626  call downsample_mask(g%mask2dBu,diag_cs%dsamp(dl)%mask2dBu, dl,g%IscB,g%JscB, &
3627  g%HId2%IscB,g%HId2%IecB,g%HId2%JscB,g%HId2%JecB,g%HId2%IsdB,g%HId2%IedB,g%HId2%JsdB,g%HId2%JedB)
3628  call downsample_mask(g%mask2dCu,diag_cs%dsamp(dl)%mask2dCu, dl,g%IscB,g%JscB, &
3629  g%HId2%IscB,g%HId2%IecB,g%HId2%jsc, g%HId2%jec,g%HId2%IsdB,g%HId2%IedB,g%HId2%jsd, g%HId2%jed)
3630  call downsample_mask(g%mask2dCv,diag_cs%dsamp(dl)%mask2dCv, dl,g%isc ,g%JscB, &
3631  g%HId2%isc ,g%HId2%iec, g%HId2%JscB,g%HId2%JecB,g%HId2%isd ,g%HId2%ied, g%HId2%JsdB,g%HId2%JedB)
3632  ! 3d native masks are needed by diag_manager but the native variables
3633  ! can only be masked 2d - for ocean points, all layers exists.
3634  allocate(diag_cs%dsamp(dl)%mask3dTL(g%HId2%isd:g%HId2%ied,g%HId2%jsd:g%HId2%jed,1:nz))
3635  allocate(diag_cs%dsamp(dl)%mask3dBL(g%HId2%IsdB:g%HId2%IedB,g%HId2%JsdB:g%HId2%JedB,1:nz))
3636  allocate(diag_cs%dsamp(dl)%mask3dCuL(g%HId2%IsdB:g%HId2%IedB,g%HId2%jsd:g%HId2%jed,1:nz))
3637  allocate(diag_cs%dsamp(dl)%mask3dCvL(g%HId2%isd:g%HId2%ied,g%HId2%JsdB:g%HId2%JedB,1:nz))
3638  do k=1,nz
3639  diag_cs%dsamp(dl)%mask3dTL(:,:,k) = diag_cs%dsamp(dl)%mask2dT(:,:)
3640  diag_cs%dsamp(dl)%mask3dBL(:,:,k) = diag_cs%dsamp(dl)%mask2dBu(:,:)
3641  diag_cs%dsamp(dl)%mask3dCuL(:,:,k) = diag_cs%dsamp(dl)%mask2dCu(:,:)
3642  diag_cs%dsamp(dl)%mask3dCvL(:,:,k) = diag_cs%dsamp(dl)%mask2dCv(:,:)
3643  enddo
3644  allocate(diag_cs%dsamp(dl)%mask3dTi(g%HId2%isd:g%HId2%ied,g%HId2%jsd:g%HId2%jed,1:nz+1))
3645  allocate(diag_cs%dsamp(dl)%mask3dBi(g%HId2%IsdB:g%HId2%IedB,g%HId2%JsdB:g%HId2%JedB,1:nz+1))
3646  allocate(diag_cs%dsamp(dl)%mask3dCui(g%HId2%IsdB:g%HId2%IedB,g%HId2%jsd:g%HId2%jed,1:nz+1))
3647  allocate(diag_cs%dsamp(dl)%mask3dCvi(g%HId2%isd:g%HId2%ied,g%HId2%JsdB:g%HId2%JedB,1:nz+1))
3648  do k=1,nz+1
3649  diag_cs%dsamp(dl)%mask3dTi(:,:,k) = diag_cs%dsamp(dl)%mask2dT(:,:)
3650  diag_cs%dsamp(dl)%mask3dBi(:,:,k) = diag_cs%dsamp(dl)%mask2dBu(:,:)
3651  diag_cs%dsamp(dl)%mask3dCui(:,:,k) = diag_cs%dsamp(dl)%mask2dCu(:,:)
3652  diag_cs%dsamp(dl)%mask3dCvi(:,:,k) = diag_cs%dsamp(dl)%mask2dCv(:,:)
3653  enddo
3654  enddo
3655 end subroutine downsample_diag_masks_set
3656 
3657 !> Get the diagnostics-compute indices (to be passed to send_data) based on the shape of
3658 !! the diag field (the same way they are deduced for non-downsampled fields)
3659 subroutine downsample_diag_indices_get(fo1, fo2, dl, diag_cs, isv, iev, jsv, jev)
3660  integer, intent(in) :: fo1 !< The size of the diag field in x
3661  integer, intent(in) :: fo2 !< The size of the diag field in y
3662  integer, intent(in) :: dl !< Integer downsample level
3663  type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output
3664  integer, intent(out) :: isv !< i-start index for diagnostics
3665  integer, intent(out) :: iev !< i-end index for diagnostics
3666  integer, intent(out) :: jsv !< j-start index for diagnostics
3667  integer, intent(out) :: jev !< j-end index for diagnostics
3668  ! Local variables
3669  integer :: dszi,cszi,dszj,cszj,f1,f2
3670  character(len=500) :: mesg
3671  logical, save :: first_check = .true.
3672 
3673  !Check ONCE that the downsampled diag-compute domain is commensurate with the original
3674  !non-downsampled diag-compute domain.
3675  !This is a major limitation of the current implementation of the downsampled diagnostics.
3676  !We assume that the compute domain can be subdivided to dl*dl cells, hence avoiding the need of halo updates.
3677  !We want this check to error out only if there was a downsampled diagnostics requested and about to post that is
3678  !why the check is here and not in the init routines. This check need to be done only once, hence the outer if.
3679  if(first_check) then
3680  if(mod(diag_cs%ie-diag_cs%is+1, dl) .ne. 0 .OR. mod(diag_cs%je-diag_cs%js+1, dl) .ne. 0) then
3681  write (mesg,*) "Non-commensurate downsampled domain is not supported. "//&
3682  "Please choose a layout such that NIGLOBAL/Layout_X and NJGLOBAL/Layout_Y are both divisible by dl=",dl,&
3683  " Current domain extents: ", diag_cs%is,diag_cs%ie, diag_cs%js,diag_cs%je
3684  call mom_error(fatal,"downsample_diag_indices_get: "//trim(mesg))
3685  endif
3686  first_check = .false.
3687  endif
3688 
3689  cszi = diag_cs%dsamp(dl)%iec-diag_cs%dsamp(dl)%isc +1 ; dszi = diag_cs%dsamp(dl)%ied-diag_cs%dsamp(dl)%isd +1
3690  cszj = diag_cs%dsamp(dl)%jec-diag_cs%dsamp(dl)%jsc +1 ; dszj = diag_cs%dsamp(dl)%jed-diag_cs%dsamp(dl)%jsd +1
3691  isv = diag_cs%dsamp(dl)%isc ; iev = diag_cs%dsamp(dl)%iec
3692  jsv = diag_cs%dsamp(dl)%jsc ; jev = diag_cs%dsamp(dl)%jec
3693  f1 = fo1/dl
3694  f2 = fo2/dl
3695  !Correction for the symmetric case
3696  if (diag_cs%G%symmetric) then
3697  f1 = f1 + mod(fo1,dl)
3698  f2 = f2 + mod(fo2,dl)
3699  endif
3700  if ( f1 == dszi ) then
3701  isv = diag_cs%dsamp(dl)%isc ; iev = diag_cs%dsamp(dl)%iec ! field on Data domain, take compute domain indcies
3702  !The rest is not taken with the full MOM6 diag_table
3703  elseif ( f1 == dszi + 1 ) then
3704  isv = diag_cs%dsamp(dl)%isc ; iev = diag_cs%dsamp(dl)%iec+1 ! Symmetric data domain
3705  elseif ( f1 == cszi) then
3706  isv = 1 ; iev = (diag_cs%dsamp(dl)%iec-diag_cs%dsamp(dl)%isc) +1 ! Computational domain
3707  elseif ( f1 == cszi + 1 ) then
3708  isv = 1 ; iev = (diag_cs%dsamp(dl)%iec-diag_cs%dsamp(dl)%isc) +2 ! Symmetric computational domain
3709  else
3710  write (mesg,*) " peculiar size ",f1," in i-direction\n"//&
3711  "does not match one of ", cszi, cszi+1, dszi, dszi+1
3712  call mom_error(fatal,"downsample_diag_indices_get: "//trim(mesg))
3713  endif
3714  if ( f2 == dszj ) then
3715  jsv = diag_cs%dsamp(dl)%jsc ; jev = diag_cs%dsamp(dl)%jec ! Data domain
3716  elseif ( f2 == dszj + 1 ) then
3717  jsv = diag_cs%dsamp(dl)%jsc ; jev = diag_cs%dsamp(dl)%jec+1 ! Symmetric data domain
3718  elseif ( f2 == cszj) then
3719  jsv = 1 ; jev = (diag_cs%dsamp(dl)%jec-diag_cs%dsamp(dl)%jsc) +1 ! Computational domain
3720  elseif ( f2 == cszj + 1 ) then
3721  jsv = 1 ; jev = (diag_cs%dsamp(dl)%jec-diag_cs%dsamp(dl)%jsc) +2 ! Symmetric computational domain
3722  else
3723  write (mesg,*) " peculiar size ",f2," in j-direction\n"//&
3724  "does not match one of ", cszj, cszj+1, dszj, dszj+1
3725  call mom_error(fatal,"downsample_diag_indices_get: "//trim(mesg))
3726  endif
3727 end subroutine downsample_diag_indices_get
3728 
3729 !> This subroutine allocates and computes a downsampled array from an input array
3730 !! It also determines the diagnostics-compurte indices for the downsampled array
3731 !! 3d interface
3732 subroutine downsample_diag_field_3d(locfield, locfield_dsamp, dl, diag_cs, diag, isv, iev, jsv, jev, mask)
3733  real, dimension(:,:,:), pointer :: locfield !< Input array pointer
3734  real, dimension(:,:,:), allocatable, intent(inout) :: locfield_dsamp !< Output (downsampled) array
3735  type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output
3736  type(diag_type), intent(in) :: diag !< A structure describing the diagnostic to post
3737  integer, intent(in) :: dl !< Level of down sampling
3738  integer, intent(inout) :: isv !< i-start index for diagnostics
3739  integer, intent(inout) :: iev !< i-end index for diagnostics
3740  integer, intent(inout) :: jsv !< j-start index for diagnostics
3741  integer, intent(inout) :: jev !< j-end index for diagnostics
3742  real, optional,target, intent(in) :: mask(:,:,:) !< If present, use this real array as the data mask.
3743  ! Locals
3744  real, dimension(:,:,:), pointer :: locmask
3745  integer :: f1,f2,isv_o,jsv_o
3746 
3747  locmask => null()
3748  !Get the correct indices corresponding to input field
3749  !Shape of the input diag field
3750  f1=size(locfield,1)
3751  f2=size(locfield,2)
3752  !Save the extents of the original (fine) domain
3753  isv_o=isv;jsv_o=jsv
3754  !Get the shape of the downsampled field and overwrite isv,iev,jsv,jev with them
3755  call downsample_diag_indices_get(f1,f2, dl, diag_cs,isv,iev,jsv,jev)
3756  !Set the non-downsampled mask, it must be associated and initialized
3757  if (present(mask)) then
3758  locmask => mask
3759  elseif (associated(diag%axes%mask3d)) then
3760  locmask => diag%axes%mask3d
3761  else
3762  call mom_error(fatal, "downsample_diag_field_3d: Cannot downsample without a mask!!! ")
3763  endif
3764 
3765  call downsample_field(locfield, locfield_dsamp, dl, diag%xyz_method, locmask, diag_cs, diag, &
3766  isv_o,jsv_o,isv,iev,jsv,jev)
3767 
3768 end subroutine downsample_diag_field_3d
3769 
3770 !> This subroutine allocates and computes a downsampled array from an input array
3771 !! It also determines the diagnostics-compurte indices for the downsampled array
3772 !! 2d interface
3773 subroutine downsample_diag_field_2d(locfield, locfield_dsamp, dl, diag_cs, diag, isv, iev, jsv, jev, mask)
3774  real, dimension(:,:), pointer :: locfield !< Input array pointer
3775  real, dimension(:,:), allocatable, intent(inout) :: locfield_dsamp !< Output (downsampled) array
3776  type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output
3777  type(diag_type), intent(in) :: diag !< A structure describing the diagnostic to post
3778  integer, intent(in) :: dl !< Level of down sampling
3779  integer, intent(inout) :: isv !< i-start index for diagnostics
3780  integer, intent(inout) :: iev !< i-end index for diagnostics
3781  integer, intent(inout) :: jsv !< j-start index for diagnostics
3782  integer, intent(inout) :: jev !< j-end index for diagnostics
3783  real, optional,target, intent(in) :: mask(:,:) !< If present, use this real array as the data mask.
3784  ! Locals
3785  real, dimension(:,:), pointer :: locmask
3786  integer :: f1,f2,isv_o,jsv_o
3787 
3788  locmask => null()
3789  !Get the correct indices corresponding to input field
3790  !Shape of the input diag field
3791  f1=size(locfield,1)
3792  f2=size(locfield,2)
3793  !Save the extents of the original (fine) domain
3794  isv_o=isv;jsv_o=jsv
3795  !Get the shape of the downsampled field and overwrite isv,iev,jsv,jev with them
3796  call downsample_diag_indices_get(f1,f2, dl, diag_cs,isv,iev,jsv,jev)
3797  !Set the non-downsampled mask, it must be associated and initialized
3798  if (present(mask)) then
3799  locmask => mask
3800  elseif (associated(diag%axes%mask2d)) then
3801  locmask => diag%axes%mask2d
3802  else
3803  call mom_error(fatal, "downsample_diag_field_2d: Cannot downsample without a mask!!! ")
3804  endif
3805 
3806  call downsample_field(locfield, locfield_dsamp, dl, diag%xyz_method, locmask, diag_cs,diag, &
3807  isv_o,jsv_o,isv,iev,jsv,jev)
3808 
3809 end subroutine downsample_diag_field_2d
3810 
3811 !> \section downsampling The down sample algorithm
3812 !!
3813 !! The down sample method could be deduced (before send_data call)
3814 !! from the diag%x_cell_method, diag%y_cell_method and diag%v_cell_method
3815 !!
3816 !! This is the summary of the down sample algoritm for a diagnostic field f:
3817 !! \f[
3818 !! f(Id,Jd) = \sum_{i,j} f(Id+i,Jd+j) * weight(Id+i,Jd+j) / [ \sum_{i,j} weight(Id+i,Jd+j)]
3819 !! \f]
3820 !! Here, i and j run from 0 to dl-1 (dl being the down sample level).
3821 !! Id,Jd are the down sampled (coarse grid) indices run over the coarsened compute grid,
3822 !! if and jf are the original (fine grid) indices.
3823 !!
3824 !! \verbatim
3825 !! Example x_cell y_cell v_cell algorithm_id implemented weight(if,jf)
3826 !! ---------------------------------------------------------------------------------------
3827 !! theta mean mean mean MMM =222 G%areaT(if,jf)*h(if,jf)
3828 !! u point mean mean PMM =022 dyCu(if,jf)*h(if,jf)*delta(if,Id)
3829 !! v mean point mean MPM =202 dxCv(if,jf)*h(if,jf)*delta(jf,Jd)
3830 !! ? point sum mean PSM =012 h(if,jf)*delta(if,Id)
3831 !! volcello sum sum sum SSS =111 1
3832 !! T_dfxy_co sum sum point SSP =110 1
3833 !! umo point sum sum PSS =011 1*delta(if,Id)
3834 !! vmo sum point sum SPS =101 1*delta(jf,Jd)
3835 !! umo_2d point sum point PSP =010 1*delta(if,Id)
3836 !! vmo_2d sum point point SPP =100 1*delta(jf,Jd)
3837 !! ? point mean point PMP =020 dyCu(if,jf)*delta(if,Id)
3838 !! ? mean point point MPP =200 dxCv(if,jf)*delta(jf,Jd)
3839 !! w mean mean point MMP =220 G%areaT(if,jf)
3840 !! h*theta mean mean sum MMS =221 G%areaT(if,jf)
3841 !!
3842 !! delta is the Kronecker delta
3843 !! \endverbatim
3844 
3845 !> This subroutine allocates and computes a down sampled 3d array given an input array
3846 !! The down sample method is based on the "cell_methods" for the diagnostics as explained
3847 !! in the above table
3848 subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, diag,isv_o,jsv_o,isv_d,iev_d,jsv_d,jev_d)
3849  real, dimension(:,:,:), pointer :: field_in !< Original field to be down sampled
3850  real, dimension(:,:,:), allocatable :: field_out !< down sampled field
3851  integer, intent(in) :: dl !< Level of down sampling
3852  integer, intent(in) :: method !< Sampling method
3853  real, dimension(:,:,:), pointer :: mask !< Mask for field
3854  type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output
3855  type(diag_type), intent(in) :: diag !< A structure describing the diagnostic to post
3856  integer, intent(in) :: isv_o !< Original i-start index
3857  integer, intent(in) :: jsv_o !< Original j-start index
3858  integer, intent(in) :: isv_d !< i-start index of down sampled data
3859  integer, intent(in) :: iev_d !< i-end index of down sampled data
3860  integer, intent(in) :: jsv_d !< j-start index of down sampled data
3861  integer, intent(in) :: jev_d !< j-end index of down sampled data
3862  ! Locals
3863  character(len=240) :: mesg
3864  integer :: i,j,ii,jj,i0,j0,f1,f2,f_in1,f_in2
3865  integer :: k,ks,ke
3866  real :: ave,total_weight,weight
3867  real :: eps_vol ! A negligibly small volume or mass [H L2 ~> m3 or kg]
3868  real :: eps_area ! A negligibly small area [L2 ~> m2]
3869  real :: eps_face ! A negligibly small face area [H L ~> m2 or kg m-1]
3870 
3871  ks = 1 ; ke = size(field_in,3)
3872  eps_face = 1.0e-20 * diag_cs%G%US%m_to_L * diag_cs%GV%m_to_H
3873  eps_area = 1.0e-20 * diag_cs%G%US%m_to_L**2
3874  eps_vol = 1.0e-20 * diag_cs%G%US%m_to_L**2 * diag_cs%GV%m_to_H
3875 
3876  ! Allocate the down sampled field on the down sampled data domain
3877 ! allocate(field_out(diag_cs%dsamp(dl)%isd:diag_cs%dsamp(dl)%ied,diag_cs%dsamp(dl)%jsd:diag_cs%dsamp(dl)%jed,ks:ke))
3878 ! allocate(field_out(1:size(field_in,1)/dl,1:size(field_in,2)/dl,ks:ke))
3879  f_in1 = size(field_in,1)
3880  f_in2 = size(field_in,2)
3881  f1 = f_in1/dl
3882  f2 = f_in2/dl
3883  !Correction for the symmetric case
3884  if (diag_cs%G%symmetric) then
3885  f1 = f1 + mod(f_in1,dl)
3886  f2 = f2 + mod(f_in2,dl)
3887  endif
3888  allocate(field_out(1:f1,1:f2,ks:ke))
3889 
3890  ! Fill the down sampled field on the down sampled diagnostics (almost always compuate) domain
3891  if (method == mmm) then
3892  do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
3893  i0 = isv_o+dl*(i-isv_d)
3894  j0 = jsv_o+dl*(j-jsv_d)
3895  ave = 0.0
3896  total_weight = 0.0
3897  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
3898 ! do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1 !This seems to be faster!!!!
3899  weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj) * diag_cs%h(ii,jj,k)
3900  total_weight = total_weight + weight
3901  ave = ave+field_in(ii,jj,k) * weight
3902  enddo; enddo
3903  field_out(i,j,k) = ave/(total_weight + eps_vol) !Avoid zero mask at all aggregating cells where ave=0.0
3904  enddo; enddo; enddo
3905  elseif (method == sss) then !e.g., volcello
3906  do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
3907  i0 = isv_o+dl*(i-isv_d)
3908  j0 = jsv_o+dl*(j-jsv_d)
3909  ave = 0.0
3910  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
3911  weight = mask(ii,jj,k)
3912  ave = ave+field_in(ii,jj,k)*weight
3913  enddo; enddo
3914  field_out(i,j,k) = ave !Masked Sum (total_weight=1)
3915  enddo; enddo; enddo
3916  elseif(method == mmp .or. method == mms) then !e.g., T_advection_xy
3917  do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
3918  i0 = isv_o+dl*(i-isv_d)
3919  j0 = jsv_o+dl*(j-jsv_d)
3920  ave = 0.0
3921  total_weight = 0.0
3922  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
3923 ! do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1
3924  weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj)
3925  total_weight = total_weight + weight
3926  ave = ave+field_in(ii,jj,k)*weight
3927  enddo; enddo
3928  field_out(i,j,k) = ave / (total_weight+eps_area) !Avoid zero mask at all aggregating cells where ave=0.0
3929  enddo; enddo; enddo
3930  elseif(method == pmm) then
3931  do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
3932  i0 = isv_o+dl*(i-isv_d)
3933  j0 = jsv_o+dl*(j-jsv_d)
3934  ave = 0.0
3935  total_weight = 0.0
3936  ii=i0
3937  do jj=j0,j0+dl-1
3938  weight =mask(ii,jj,k) * diag_cs%G%dyCu(ii,jj) * diag_cs%h(ii,jj,k)
3939  total_weight = total_weight +weight
3940  ave=ave+field_in(ii,jj,k)*weight
3941  enddo
3942  field_out(i,j,k) = ave/(total_weight+eps_face) !Avoid zero mask at all aggregating cells where ave=0.0
3943  enddo; enddo; enddo
3944  elseif(method == pss) then !e.g. umo
3945  do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
3946  i0 = isv_o+dl*(i-isv_d)
3947  j0 = jsv_o+dl*(j-jsv_d)
3948  ave = 0.0
3949  ii=i0
3950  do jj=j0,j0+dl-1
3951  weight =mask(ii,jj,k)
3952  ave=ave+field_in(ii,jj,k)*weight
3953  enddo
3954  field_out(i,j,k) = ave !Masked Sum (total_weight=1)
3955  enddo; enddo; enddo
3956  elseif(method == sps) then !e.g. vmo
3957  do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
3958  i0 = isv_o+dl*(i-isv_d)
3959  j0 = jsv_o+dl*(j-jsv_d)
3960  ave = 0.0
3961  jj=j0
3962  do ii=i0,i0+dl-1
3963  weight =mask(ii,jj,k)
3964  ave=ave+field_in(ii,jj,k)*weight
3965  enddo
3966  field_out(i,j,k) = ave !Masked Sum (total_weight=1)
3967  enddo; enddo; enddo
3968  elseif(method == mpm) then
3969  do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
3970  i0 = isv_o+dl*(i-isv_d)
3971  j0 = jsv_o+dl*(j-jsv_d)
3972  ave = 0.0
3973  total_weight = 0.0
3974  jj=j0
3975  do ii=i0,i0+dl-1
3976  weight = mask(ii,jj,k) * diag_cs%G%dxCv(ii,jj) * diag_cs%h(ii,jj,k)
3977  total_weight = total_weight + weight
3978  ave=ave+field_in(ii,jj,k)*weight
3979  enddo
3980  field_out(i,j,k) = ave/(total_weight+eps_face) !Avoid zero mask at all aggregating cells where ave=0.0
3981  enddo; enddo; enddo
3982  elseif(method == msk) then !The input field is a mask, subsample
3983  field_out(:,:,:) = 0.0
3984  do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
3985  i0 = isv_o+dl*(i-isv_d)
3986  j0 = jsv_o+dl*(j-jsv_d)
3987  ave = 0.0
3988  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
3989  ave=ave+field_in(ii,jj,k)
3990  enddo; enddo
3991  if(ave > 0.0) field_out(i,j,k)=1.0
3992  enddo; enddo; enddo
3993  else
3994  write (mesg,*) " unknown sampling method: ",method
3995  call mom_error(fatal, "downsample_field_3d: "//trim(mesg)//" "//trim(diag%debug_str))
3996  endif
3997 
3998 end subroutine downsample_field_3d
3999 
4000 !> This subroutine allocates and computes a down sampled 2d array given an input array
4001 !! The down sample method is based on the "cell_methods" for the diagnostics as explained
4002 !! in the above table
4003 subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, diag, &
4004  isv_o, jsv_o, isv_d, iev_d, jsv_d, jev_d)
4005  real, dimension(:,:), pointer :: field_in !< Original field to be down sampled
4006  real, dimension(:,:), allocatable :: field_out !< Down sampled field
4007  integer, intent(in) :: dl !< Level of down sampling
4008  integer, intent(in) :: method !< Sampling method
4009  real, dimension(:,:), pointer :: mask !< Mask for field
4010  type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output
4011  type(diag_type), intent(in) :: diag !< A structure describing the diagnostic to post
4012  integer, intent(in) :: isv_o !< Original i-start index
4013  integer, intent(in) :: jsv_o !< Original j-start index
4014  integer, intent(in) :: isv_d !< i-start index of down sampled data
4015  integer, intent(in) :: iev_d !< i-end index of down sampled data
4016  integer, intent(in) :: jsv_d !< j-start index of down sampled data
4017  integer, intent(in) :: jev_d !< j-end index of down sampled data
4018  ! Locals
4019  character(len=240) :: mesg
4020  integer :: i,j,ii,jj,i0,j0,f1,f2,f_in1,f_in2
4021  real :: ave, total_weight, weight
4022  real :: epsilon = 1.0e-20 ! A negligibly small count of weights [nondim]
4023  real :: eps_area ! A negligibly small area [L2 ~> m2]
4024  real :: eps_len ! A negligibly small horizontal length [L ~> m]
4025 
4026  eps_len = 1.0e-20 * diag_cs%G%US%m_to_L
4027  eps_area = 1.0e-20 * diag_cs%G%US%m_to_L**2
4028 
4029  ! Allocate the down sampled field on the down sampled data domain
4030 ! allocate(field_out(diag_cs%dsamp(dl)%isd:diag_cs%dsamp(dl)%ied,diag_cs%dsamp(dl)%jsd:diag_cs%dsamp(dl)%jed))
4031 ! allocate(field_out(1:size(field_in,1)/dl,1:size(field_in,2)/dl))
4032  ! Fill the down sampled field on the down sampled diagnostics (almost always compuate) domain
4033  f_in1 = size(field_in,1)
4034  f_in2 = size(field_in,2)
4035  f1 = f_in1/dl
4036  f2 = f_in2/dl
4037  ! Correction for the symmetric case
4038  if (diag_cs%G%symmetric) then
4039  f1 = f1 + mod(f_in1,dl)
4040  f2 = f2 + mod(f_in2,dl)
4041  endif
4042  allocate(field_out(1:f1,1:f2))
4043 
4044  if (method == mmp) then
4045  do j=jsv_d,jev_d ; do i=isv_d,iev_d
4046  i0 = isv_o+dl*(i-isv_d)
4047  j0 = jsv_o+dl*(j-jsv_d)
4048  ave = 0.0
4049  total_weight = 0.0
4050  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
4051 ! do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1
4052  weight = mask(ii,jj)*diag_cs%G%areaT(ii,jj)
4053  total_weight = total_weight + weight
4054  ave = ave+field_in(ii,jj)*weight
4055  enddo; enddo
4056  field_out(i,j) = ave/(total_weight + eps_area) !Avoid zero mask at all aggregating cells where ave=0.0
4057  enddo; enddo
4058  elseif(method == ssp) then ! e.g., T_dfxy_cont_tendency_2d
4059  do j=jsv_d,jev_d ; do i=isv_d,iev_d
4060  i0 = isv_o+dl*(i-isv_d)
4061  j0 = jsv_o+dl*(j-jsv_d)
4062  ave = 0.0
4063  total_weight = 0.0
4064  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
4065 ! do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1
4066  weight = mask(ii,jj)
4067  total_weight = total_weight + weight
4068  ave=ave+field_in(ii,jj)*weight
4069  enddo; enddo
4070  field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0
4071  enddo; enddo
4072  elseif(method == psp) then ! e.g., umo_2d
4073  do j=jsv_d,jev_d ; do i=isv_d,iev_d
4074  i0 = isv_o+dl*(i-isv_d)
4075  j0 = jsv_o+dl*(j-jsv_d)
4076  ave = 0.0
4077  total_weight = 0.0
4078  ii=i0
4079  do jj=j0,j0+dl-1
4080  weight = mask(ii,jj)
4081  total_weight = total_weight +weight
4082  ave=ave+field_in(ii,jj)*weight
4083  enddo
4084  field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0
4085  enddo; enddo
4086  elseif(method == spp) then ! e.g., vmo_2d
4087  do j=jsv_d,jev_d ; do i=isv_d,iev_d
4088  i0 = isv_o+dl*(i-isv_d)
4089  j0 = jsv_o+dl*(j-jsv_d)
4090  ave = 0.0
4091  total_weight = 0.0
4092  jj=j0
4093  do ii=i0,i0+dl-1
4094  weight = mask(ii,jj)
4095  total_weight = total_weight +weight
4096  ave=ave+field_in(ii,jj)*weight
4097  enddo
4098  field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0
4099  enddo; enddo
4100  elseif(method == pmp) then
4101  do j=jsv_d,jev_d ; do i=isv_d,iev_d
4102  i0 = isv_o+dl*(i-isv_d)
4103  j0 = jsv_o+dl*(j-jsv_d)
4104  ave = 0.0
4105  total_weight = 0.0
4106  ii=i0
4107  do jj=j0,j0+dl-1
4108  weight = mask(ii,jj) * diag_cs%G%dyCu(ii,jj)!*diag_cs%h(ii,jj,1) !Niki?
4109  total_weight = total_weight +weight
4110  ave=ave+field_in(ii,jj)*weight
4111  enddo
4112  field_out(i,j) = ave/(total_weight+eps_len) !Avoid zero mask at all aggregating cells where ave=0.0
4113  enddo; enddo
4114  elseif(method == mpp) then
4115  do j=jsv_d,jev_d ; do i=isv_d,iev_d
4116  i0 = isv_o+dl*(i-isv_d)
4117  j0 = jsv_o+dl*(j-jsv_d)
4118  ave = 0.0
4119  total_weight = 0.0
4120  jj=j0
4121  do ii=i0,i0+dl-1
4122  weight = mask(ii,jj)* diag_cs%G%dxCv(ii,jj)!*diag_cs%h(ii,jj,1) !Niki?
4123  total_weight = total_weight +weight
4124  ave=ave+field_in(ii,jj)*weight
4125  enddo
4126  field_out(i,j) = ave/(total_weight+eps_len) !Avoid zero mask at all aggregating cells where ave=0.0
4127  enddo; enddo
4128  elseif(method == msk) then !The input field is a mask, subsample
4129  field_out(:,:) = 0.0
4130  do j=jsv_d,jev_d ; do i=isv_d,iev_d
4131  i0 = isv_o+dl*(i-isv_d)
4132  j0 = jsv_o+dl*(j-jsv_d)
4133  ave = 0.0
4134  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
4135  ave=ave+field_in(ii,jj)
4136  enddo; enddo
4137  if(ave > 0.0) field_out(i,j)=1.0
4138  enddo; enddo
4139  else
4140  write (mesg,*) " unknown sampling method: ",method
4141  call mom_error(fatal, "downsample_field_2d: "//trim(mesg)//" "//trim(diag%debug_str))
4142  endif
4143 
4144 end subroutine downsample_field_2d
4145 
4146 !> Allocate and compute the 2d down sampled mask
4147 !! The masks are down sampled based on a minority rule, i.e., a coarse cell is open (1)
4148 !! if at least one of the sub-cells are open, otherwise it's closed (0)
4149 subroutine downsample_mask_2d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_d, jsc_d, jec_d, &
4150  isd_d, ied_d, jsd_d, jed_d)
4151  real, dimension(:,:), intent(in) :: field_in !< Original field to be down sampled
4152  real, dimension(:,:), pointer :: field_out !< Down sampled field
4153  integer, intent(in) :: dl !< Level of down sampling
4154  integer, intent(in) :: isc_o !< Original i-start index
4155  integer, intent(in) :: jsc_o !< Original j-start index
4156  integer, intent(in) :: isc_d !< Computational i-start index of down sampled data
4157  integer, intent(in) :: iec_d !< Computational i-end index of down sampled data
4158  integer, intent(in) :: jsc_d !< Computational j-start index of down sampled data
4159  integer, intent(in) :: jec_d !< Computational j-end index of down sampled data
4160  integer, intent(in) :: isd_d !< Computational i-start index of down sampled data
4161  integer, intent(in) :: ied_d !< Computational i-end index of down sampled data
4162  integer, intent(in) :: jsd_d !< Computational j-start index of down sampled data
4163  integer, intent(in) :: jed_d !< Computational j-end index of down sampled data
4164  ! Locals
4165  integer :: i,j,ii,jj,i0,j0
4166  real :: tot_non_zero
4167  ! down sampled mask = 0 unless the mask value of one of the down sampling cells is 1
4168  allocate(field_out(isd_d:ied_d,jsd_d:jed_d))
4169  field_out(:,:) = 0.0
4170  do j=jsc_d,jec_d ; do i=isc_d,iec_d
4171  i0 = isc_o+dl*(i-isc_d)
4172  j0 = jsc_o+dl*(j-jsc_d)
4173  tot_non_zero = 0.0
4174  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
4175  tot_non_zero = tot_non_zero + field_in(ii,jj)
4176  enddo;enddo
4177  if(tot_non_zero > 0.0) field_out(i,j)=1.0
4178  enddo; enddo
4179 end subroutine downsample_mask_2d
4180 
4181 !> Allocate and compute the 3d down sampled mask
4182 !! The masks are down sampled based on a minority rule, i.e., a coarse cell is open (1)
4183 !! if at least one of the sub-cells are open, otherwise it's closed (0)
4184 subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_d, jsc_d, jec_d, &
4185  isd_d, ied_d, jsd_d, jed_d)
4186  real, dimension(:,:,:), intent(in) :: field_in !< Original field to be down sampled
4187  real, dimension(:,:,:), pointer :: field_out !< down sampled field
4188  integer, intent(in) :: dl !< Level of down sampling
4189  integer, intent(in) :: isc_o !< Original i-start index
4190  integer, intent(in) :: jsc_o !< Original j-start index
4191  integer, intent(in) :: isc_d !< Computational i-start index of down sampled data
4192  integer, intent(in) :: iec_d !< Computational i-end index of down sampled data
4193  integer, intent(in) :: jsc_d !< Computational j-start index of down sampled data
4194  integer, intent(in) :: jec_d !< Computational j-end index of down sampled data
4195  integer, intent(in) :: isd_d !< Computational i-start index of down sampled data
4196  integer, intent(in) :: ied_d !< Computational i-end index of down sampled data
4197  integer, intent(in) :: jsd_d !< Computational j-start index of down sampled data
4198  integer, intent(in) :: jed_d !< Computational j-end index of down sampled data
4199  ! Locals
4200  integer :: i,j,ii,jj,i0,j0,k,ks,ke
4201  real :: tot_non_zero
4202  ! down sampled mask = 0 unless the mask value of one of the down sampling cells is 1
4203  ks = lbound(field_in,3) ; ke = ubound(field_in,3)
4204  allocate(field_out(isd_d:ied_d,jsd_d:jed_d,ks:ke))
4205  field_out(:,:,:) = 0.0
4206  do k= ks,ke ; do j=jsc_d,jec_d ; do i=isc_d,iec_d
4207  i0 = isc_o+dl*(i-isc_d)
4208  j0 = jsc_o+dl*(j-jsc_d)
4209  tot_non_zero = 0.0
4210  do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
4211  tot_non_zero = tot_non_zero + field_in(ii,jj,k)
4212  enddo;enddo
4213  if(tot_non_zero > 0.0) field_out(i,j,k)=1.0
4214  enddo; enddo; enddo
4215 end subroutine downsample_mask_3d
4216 
4217 end module mom_diag_mediator
4218 
mom_diag_mediator::define_axes_group_dsamp
subroutine define_axes_group_dsamp(diag_cs, handles, axes, dl, nz, vertical_coordinate_number, x_cell_method, y_cell_method, v_cell_method, is_h_point, is_q_point, is_u_point, is_v_point, is_layer, is_interface, is_native, needs_remapping, needs_interpolating, xyave_axes)
Defines a group of downsampled "axes" from list of handles.
Definition: MOM_diag_mediator.F90:1054
mom_time_manager
Wraps the FMS time manager functions.
Definition: MOM_time_manager.F90:2
mom_diag_mediator::diag_restore_grids
subroutine, public diag_restore_grids(diag)
Restore the diagnostic grids from the temporary structure within diag.
Definition: MOM_diag_mediator.F90:3564
mom_diag_remap::diag_remap_diag_registration_closed
subroutine, public diag_remap_diag_registration_closed(remap_cs)
Inform that all diagnostics have been registered. If _set_active() has not been called on the remappi...
Definition: MOM_diag_remap.F90:159
mom_diag_mediator::diag_get_volume_cell_measure_dm_id
integer function, public diag_get_volume_cell_measure_dm_id(diag_cs)
Returns diag_manager id for cell measure of h-cells.
Definition: MOM_diag_mediator.F90:930
mom_diag_mediator::post_data_1d_k
subroutine, public post_data_1d_k(diag_field_id, field, diag_cs, is_static)
Make a real 1-d array diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:1233
mom_diag_remap::vertically_reintegrate_diag_field
subroutine, public vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, mask, missing_value, field, reintegrated_field)
Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid.
Definition: MOM_diag_remap.F90:479
mom_diag_mediator::mmm
integer mmm
x:mean,y:mean,z:mean
Definition: MOM_diag_mediator.F90:169
mom_diag_mediator::mpm
integer mpm
x:mean,y:point,z:mean
Definition: MOM_diag_mediator.F90:165
mom_diag_mediator::pmm
integer pmm
x:point,y:mean,z:mean
Definition: MOM_diag_mediator.F90:160
mom_safe_alloc
Convenience functions for safely allocating memory without accidentally reallocating pointer and caus...
Definition: MOM_safe_alloc.F90:3
mom_diag_mediator::diagcs_dsamp
Container for down sampling information.
Definition: MOM_diag_mediator.F90:195
mom_io::query_vardesc
subroutine, public query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, cmor_field_name, cmor_units, cmor_longname, conversion, caller)
This routine queries vardesc.
Definition: MOM_io.F90:699
mom_diag_mediator::register_diag_field_expand_cmor
logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, axes, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns True if either the native or CMOr version of the diagnostic were registered....
Definition: MOM_diag_mediator.F90:2121
mom_diag_mediator::downsample_field
Down sample a field.
Definition: MOM_diag_mediator.F90:75
mom_diag_mediator::query_averaging_enabled
logical function, public query_averaging_enabled(diag_cs, time_int, time_end)
Call this subroutine to determine whether the averaging is currently enabled. .true....
Definition: MOM_diag_mediator.F90:1850
mom_diag_mediator::ssp
integer ssp
x:sum;y:sum,z:point
Definition: MOM_diag_mediator.F90:163
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_checksums::bchksum
Checksums an array (2d or 3d) staggered at corner points.
Definition: MOM_checksums.F90:53
mom_checksums::vchksum
Checksums an array (2d or 3d) staggered at C-grid v points.
Definition: MOM_checksums.F90:38
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator::set_masks_for_axes_dsamp
subroutine set_masks_for_axes_dsamp(G, diag_cs)
Definition: MOM_diag_mediator.F90:804
mom_diag_mediator::post_data_3d_low
subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask)
Make a real 3-d array diagnostic available for averaging or output using a diag_type instead of an in...
Definition: MOM_diag_mediator.F90:1579
mom_diag_mediator::mpp
integer mpp
x:mean,y:point,z:point
Definition: MOM_diag_mediator.F90:164
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_diag_remap::diag_remap_end
subroutine, public diag_remap_end(remap_cs)
De-init a diagnostic remapping type. Free allocated memory.
Definition: MOM_diag_remap.F90:142
mom_diag_mediator::mmp
integer mmp
x:mean,y:mean,z:point
Definition: MOM_diag_mediator.F90:166
mom_diag_remap::diag_remap_init
subroutine, public diag_remap_init(remap_cs, coord_tuple)
Initialize a diagnostic remapping type with the given vertical coordinate.
Definition: MOM_diag_remap.F90:124
mom_diag_mediator::axes_grp
A group of 1D axes that comprise a 1D/2D/3D mesh.
Definition: MOM_diag_mediator.F90:96
mom_diag_mediator::downsample_diag_indices_get
subroutine downsample_diag_indices_get(fo1, fo2, dl, diag_cs, isv, iev, jsv, jev)
Get the diagnostics-compute indices (to be passed to send_data) based on the shape of the diag field ...
Definition: MOM_diag_mediator.F90:3660
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_diag_mediator::msk
integer msk
Use the downsample method of a mask.
Definition: MOM_diag_mediator.F90:170
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
mom_diag_mediator::diag_associate_volume_cell_measure
subroutine, public diag_associate_volume_cell_measure(diag_cs, id_h_volume)
Attaches the id of cell volumes to axes groups for use with cell_measures.
Definition: MOM_diag_mediator.F90:908
mom_diag_mediator::diag_save_grids
subroutine, public diag_save_grids(diag)
Save the current diagnostic grids in the temporary structure within diag.
Definition: MOM_diag_mediator.F90:3548
mom_diag_mediator::id_clock_diag_grid_updates
integer id_clock_diag_grid_updates
Sets up diagnostics axes.
Definition: MOM_diag_mediator.F90:338
mom_diag_mediator::alloc_diag_with_id
subroutine alloc_diag_with_id(diag_id, diag_cs, diag)
Make a new diagnostic. Either use memory which is in the array of 'primary' diagnostics,...
Definition: MOM_diag_mediator.F90:3421
mom_diag_mediator::set_axes_info
subroutine, public set_axes_info(G, GV, US, param_file, diag_cs, set_vertical)
Sets up diagnostics axes.
Definition: MOM_diag_mediator.F90:344
mom_diag_mediator::diag_mediator_init
subroutine, public diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)
diag_mediator_init initializes the MOM diag_mediator and opens the available diagnostics file,...
Definition: MOM_diag_mediator.F90:2972
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_diag_mediator::downsample_field_3d
subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, diag, isv_o, jsv_o, isv_d, iev_d, jsv_d, jev_d)
This subroutine allocates and computes a down sampled 3d array given an input array The down sample m...
Definition: MOM_diag_mediator.F90:3849
mom_diag_mediator::get_new_diag_id
integer function get_new_diag_id(diag_cs)
Returns a new diagnostic id, it may be necessary to expand the diagnostics array.
Definition: MOM_diag_mediator.F90:3377
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_checksums::chksum
This is an older interface for 1-, 2-, or 3-D checksums.
Definition: MOM_checksums.F90:63
mom_diag_remap::diag_remap_do_remap
subroutine, public diag_remap_do_remap(remap_cs, G, GV, h, staggered_in_x, staggered_in_y, mask, missing_value, field, remapped_field)
Remap diagnostic field to alternative vertical grid.
Definition: MOM_diag_remap.F90:340
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_diag_mediator::diag_mediator_close_registration
subroutine, public diag_mediator_close_registration(diag_CS)
Definition: MOM_diag_mediator.F90:3284
mom_diag_mediator::downsample_diag_masks_set
subroutine downsample_diag_masks_set(G, nz, diag_cs)
Definition: MOM_diag_mediator.F90:3601
mom_diag_mediator::id_clock_diag_mediator
integer id_clock_diag_mediator
Sets up diagnostics axes.
Definition: MOM_diag_mediator.F90:338
mom_diag_mediator::diag_grid_storage_end
subroutine, public diag_grid_storage_end(grid_storage)
Deallocates the fields in the remapping fields container.
Definition: MOM_diag_mediator.F90:3581
mom_diag_mediator::post_data_2d
subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask)
Make a real 2-d array diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:1287
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_diag_mediator::diag_register_area_ids
subroutine, public diag_register_area_ids(diag_cs, id_area_t, id_area_q)
Attaches the id of cell areas to axes groups for use with cell_measures.
Definition: MOM_diag_mediator.F90:864
mom_diag_mediator::attach_cell_methods
subroutine attach_cell_methods(id, axes, ostring, cell_methods, x_cell_method, y_cell_method, v_cell_method, v_extensive)
Attaches "cell_methods" attribute to a variable based on defaults for axes_grp or optional arguments.
Definition: MOM_diag_mediator.F90:2477
mom_diag_mediator::ocean_register_diag
integer function, public ocean_register_diag(var_desc, G, diag_CS, day)
Registers a diagnostic using the information encapsulated in the vardesc type argument and returns an...
Definition: MOM_diag_mediator.F90:2860
mom_diag_mediator::downsample_field_2d
subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, diag, isv_o, jsv_o, isv_d, iev_d, jsv_d, jev_d)
This subroutine allocates and computes a down sampled 2d array given an input array The down sample m...
Definition: MOM_diag_mediator.F90:4005
mom_diag_remap::diag_remap_get_axes_info
subroutine, public diag_remap_get_axes_info(remap_cs, nz, id_layer, id_interface)
Get layer and interface axes ids for this coordinate Needed when defining axes groups.
Definition: MOM_diag_remap.F90:239
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::downsample_mask_3d
subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_d, jsc_d, jec_d, isd_d, ied_d, jsd_d, jed_d)
Allocate and compute the 3d down sampled mask The masks are down sampled based on a minority rule,...
Definition: MOM_diag_mediator.F90:4186
mom_diag_mediator::diag_grids_type
Contains an array to store a diagnostic target grid.
Definition: MOM_diag_mediator.F90:141
mom_coms
Interfaces to non-domain-oriented communication subroutines, including the MOM6 reproducing sums faci...
Definition: MOM_coms.F90:3
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_safe_alloc::safe_alloc_alloc
Allocate a 2-d or 3-d allocatable array.
Definition: MOM_safe_alloc.F90:18
mom_diag_mediator::mms
integer mms
x:mean,y:mean,z:sum
Definition: MOM_diag_mediator.F90:167
mom_diag_mediator::downsample_diag_field_2d
subroutine downsample_diag_field_2d(locfield, locfield_dsamp, dl, diag_cs, diag, isv, iev, jsv, jev, mask)
This subroutine allocates and computes a downsampled array from an input array It also determines the...
Definition: MOM_diag_mediator.F90:3774
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
mom_diag_mediator::diag_masks_set
subroutine, public diag_masks_set(G, nz, diag_cs)
Sets up the 2d and 3d masks for native diagnostics.
Definition: MOM_diag_mediator.F90:3242
mom_diag_mediator::diag_mediator_end
subroutine, public diag_mediator_end(time, diag_CS, end_diag_manager)
Definition: MOM_diag_mediator.F90:3299
mom_diag_mediator::diag_grid_storage
Stores all the remapping grids and the model's native space thicknesses.
Definition: MOM_diag_mediator.F90:146
mom_diag_mediator::psm
integer psm
x:point,y:sum,z:mean
Definition: MOM_diag_mediator.F90:158
mom_diag_mediator::diag_type
This type is used to represent a diagnostic at the diag_mediator level.
Definition: MOM_diag_mediator.F90:179
mom_diag_mediator::diag_update_remap_grids
subroutine, public diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S)
Build/update vertical grids for diagnostic remapping.
Definition: MOM_diag_mediator.F90:3187
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_checksums
Routines to calculate checksums of various array and vector types.
Definition: MOM_checksums.F90:2
mom_diag_mediator::register_static_field
integer function, public register_static_field(module_name, field_name, axes, long_name, units, missing_value, range, mask_variant, standard_name, do_not_log, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, area, x_cell_method, y_cell_method, area_cell_method, conversion)
Registers a static diagnostic, returning an integer handle.
Definition: MOM_diag_mediator.F90:2700
mom_diag_mediator::log_available_diag
subroutine log_available_diag(used, module_name, field_name, cell_methods_string, comment, diag_CS, long_name, units, standard_name)
Log a diagnostic to the available diagnostics file.
Definition: MOM_diag_mediator.F90:3442
mom_io::mom_read_data
Read a data field from a file.
Definition: MOM_io.F90:74
mom_diag_mediator::sps
integer sps
x:sum,y:point,z:sum
Definition: MOM_diag_mediator.F90:162
mom_diag_mediator::downsample_mask
Down sample the mask of a field.
Definition: MOM_diag_mediator.F90:80
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_checksums::chksum0
subroutine, public chksum0(scalar, mesg, scale, logunit)
Checksum a scalar field (consistent with array checksums)
Definition: MOM_checksums.F90:88
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_diag_mediator::get_diag_time_end
type(time_type) function, public get_diag_time_end(diag_cs)
This function returns the valid end time for use with diagnostics that are handled outside of the MOM...
Definition: MOM_diag_mediator.F90:1863
mom_diag_mediator::initialize_diag_type
subroutine initialize_diag_type(diag)
Initializes a diag_type (used after allocating new memory)
Definition: MOM_diag_mediator.F90:3408
mom_diag_mediator::psp
integer psp
x:point,y:sum,z:point
Definition: MOM_diag_mediator.F90:156
mom_diag_manager_wrapper::register_diag_field_fms
A wrapper for register_diag_field_array()
Definition: MOM_diag_manager_wrapper.F90:14
mom_diag_remap
provides runtime remapping of diagnostics to z star, sigma and rho vertical coordinates.
Definition: MOM_diag_remap.F90:56
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_diag_mediator::disable_averaging
subroutine, public disable_averaging(diag_cs)
Call this subroutine to avoid averaging any offered fields.
Definition: MOM_diag_mediator.F90:1840
mom_diag_mediator::downsample_diag_field_3d
subroutine downsample_diag_field_3d(locfield, locfield_dsamp, dl, diag_cs, diag, isv, iev, jsv, jev, mask)
This subroutine allocates and computes a downsampled array from an input array It also determines the...
Definition: MOM_diag_mediator.F90:3733
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_diag_mediator::enable_averages
subroutine, public enable_averages(time_int, time_end, diag_CS, T_to_s)
Enable the accumulation of time averages over the specified time interval in time units.
Definition: MOM_diag_mediator.F90:1820
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_error_handler::assert
subroutine, public assert(logical_arg, msg)
Issues a FATAL error if the assertion fails, i.e. the first argument is false.
Definition: MOM_error_handler.F90:182
mom_diag_mediator::diag_set_state_ptrs
subroutine, public diag_set_state_ptrs(h, T, S, eqn_of_state, diag_cs)
Set pointers to the default state fields used to remap diagnostics.
Definition: MOM_diag_mediator.F90:3169
mom_diag_remap::horizontally_average_diag_field
subroutine, public horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_in_y, is_layer, is_extensive, missing_value, field, averaged_field, averaged_mask)
Horizontally average field.
Definition: MOM_diag_remap.F90:644
mom_diag_mediator::post_data_2d_low
subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask)
Make a real 2-d array diagnostic available for averaging or output using a diag_type instead of an in...
Definition: MOM_diag_mediator.F90:1314
mom_diag_mediator::diag_copy_diag_to_storage
subroutine, public diag_copy_diag_to_storage(grid_storage, h_state, diag)
Copy from the main diagnostic arrays to the grid storage as well as the native thicknesses.
Definition: MOM_diag_mediator.F90:3511
mom_io::vardesc
Type for describing a variable, typically a tracer.
Definition: MOM_io.F90:53
mom_diag_mediator::register_cell_measure
subroutine, public register_cell_measure(G, diag, Time)
Sets a handle inside diagnostics mediator to associate 3d cell measures.
Definition: MOM_diag_mediator.F90:893
mom_checksums::zchksum
subroutine, public zchksum(array, mesg, scale, logunit)
Checksum a 1d array (typically a column).
Definition: MOM_checksums.F90:121
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_diag_mediator::define_axes_group
subroutine, public define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_number, x_cell_method, y_cell_method, v_cell_method, is_h_point, is_q_point, is_u_point, is_v_point, is_layer, is_interface, is_native, needs_remapping, needs_interpolating, xyave_axes)
Defines a group of "axes" from list of handles.
Definition: MOM_diag_mediator.F90:943
mom_diag_mediator::spp
integer spp
x:sum,y:point,z:point
Definition: MOM_diag_mediator.F90:161
mom_checksums::uchksum
Checksums an array (2d or 3d) staggered at C-grid u points.
Definition: MOM_checksums.F90:33
mom_diag_mediator::downsample_diag_field
Down sample a diagnostic field.
Definition: MOM_diag_mediator.F90:85
mom_diag_mediator::diag_dsamp
Contained for down sampled masks.
Definition: MOM_diag_mediator.F90:90
mom_diag_remap::vertically_interpolate_diag_field
subroutine, public vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, mask, missing_value, field, interpolated_field)
Vertically interpolate diagnostic field to alternative vertical grid.
Definition: MOM_diag_remap.F90:560
mom_diag_mediator::register_scalar_field
integer function, public register_scalar_field(module_name, field_name, init_time, diag_cs, long_name, units, missing_value, range, standard_name, do_not_log, err_msg, interp_method, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name)
Definition: MOM_diag_mediator.F90:2596
mom_diag_remap::diag_remap_calc_hmask
subroutine, public diag_remap_calc_hmask(remap_cs, G, mask)
Calculate masks for target grid.
Definition: MOM_diag_remap.F90:435
mom_diag_mediator::describe_option
subroutine describe_option(opt_name, value, diag_CS)
Describe an option setting in the diagnostic files.
Definition: MOM_diag_mediator.F90:2843
mom_string_functions::lowercase
character(len=len(input_string)) function, public lowercase(input_string)
Return a string in which all uppercase letters have been replaced by their lowercase counterparts.
Definition: MOM_string_functions.F90:24
mom_diag_mediator::add_xyz_method
subroutine add_xyz_method(diag, axes, x_cell_method, y_cell_method, v_cell_method, v_extensive)
Adds the encoded "cell_methods" for a diagnostics as a diag% property This allows access to the cell_...
Definition: MOM_diag_mediator.F90:2414
mom_diag_remap::diag_remap_configure_axes
subroutine, public diag_remap_configure_axes(remap_cs, GV, US, param_file)
Configure the vertical axes for a diagnostic remapping control structure. Reads a configuration param...
Definition: MOM_diag_remap.F90:180
mom_checksums::hchksum
Checksums an array (2d or 3d) staggered at tracer points.
Definition: MOM_checksums.F90:48
mom_diag_mediator::pmp
integer pmp
x:point,y:mean,z:point
Definition: MOM_diag_mediator.F90:159
mom_diag_mediator::set_diag_mediator_grid
subroutine, public set_diag_mediator_grid(G, diag_cs)
Set up the array extents for doing diagnostics.
Definition: MOM_diag_mediator.F90:1187
mom_diag_mediator::post_data_3d
subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h)
Make a real 3-d array diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:1458
mom_diag_mediator::enable_averaging
subroutine, public enable_averaging(time_int_in, time_end_in, diag_cs)
This subroutine enables the accumulation of time averages over the specified time interval.
Definition: MOM_diag_mediator.F90:1805
mom_diag_remap::diag_remap_set_active
subroutine, public diag_remap_set_active(remap_cs)
Indicate that this remapping type is actually used by the diag manager. If this is never called then ...
Definition: MOM_diag_remap.F90:171
mom_diag_remap::diag_remap_update
subroutine, public diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state)
Build/update target vertical grids for diagnostic remapping.
Definition: MOM_diag_remap.F90:268
mom_diag_mediator::add_diag_to_list
subroutine add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name, msg)
Create a diagnostic type and attached to list.
Definition: MOM_diag_mediator.F90:2388
mom_diag_mediator::post_xy_average
subroutine post_xy_average(diag_cs, diag, field)
Post the horizontally area-averaged diagnostic.
Definition: MOM_diag_mediator.F90:1751
mom_diag_mediator::i2s
character(len=15) function i2s(a, n_in)
Convert the first n elements (up to 3) of an integer array to an underscore delimited string.
Definition: MOM_diag_mediator.F90:3355
mom_diag_mediator::log_chksum_diag
subroutine log_chksum_diag(docunit, description, chksum)
Log the diagnostic chksum to the chksum diag file.
Definition: MOM_diag_mediator.F90:3475
mom_diag_mediator::sss
integer sss
x:sum,y:sum,z:sum
Definition: MOM_diag_mediator.F90:168
mom_diag_remap::diag_remap_ctrl
Represents remapping of diagnostics to a particular vertical coordinate.
Definition: MOM_diag_remap.F90:103
mom_safe_alloc::safe_alloc_ptr
Allocate a pointer to a 1-d, 2-d or 3-d array.
Definition: MOM_safe_alloc.F90:12
mom_diag_remap::diag_remap_axes_configured
logical function, public diag_remap_axes_configured(remap_cs)
Whether or not the axes for this vertical coordinated has been configured. Configuration is complete ...
Definition: MOM_diag_remap.F90:255
mom_diag_manager_wrapper
A simple (very thin) wrapper for register_diag_field to avoid a compiler bug with PGI.
Definition: MOM_diag_manager_wrapper.F90:2
mom_diag_mediator::diag_mediator_infrastructure_init
subroutine, public diag_mediator_infrastructure_init(err_msg)
Definition: MOM_diag_mediator.F90:2963
mom_diag_mediator::pss
integer pss
x:point,y:sum,z:point
Definition: MOM_diag_mediator.F90:157
mom_diag_mediator::register_diag_field_expand_axes
integer function register_diag_field_expand_axes(module_name, field_name, axes, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count)
Returns an FMS id from register_diag_field_fms (the diag_manager routine) after expanding axes (axes-...
Definition: MOM_diag_mediator.F90:2286
mom_diag_mediator::downsample_mask_2d
subroutine downsample_mask_2d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_d, jsc_d, jec_d, isd_d, ied_d, jsd_d, jed_d)
Allocate and compute the 2d down sampled mask The masks are down sampled based on a minority rule,...
Definition: MOM_diag_mediator.F90:4151
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_diag_mediator::post_data_0d
subroutine post_data_0d(diag_field_id, field, diag_cs, is_static)
Make a real scalar diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:1199
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::set_masks_for_axes
subroutine, public set_masks_for_axes(G, diag_cs)
set_masks_for_axes sets up the 2d and 3d masks for diagnostics using the current grid recorded after ...
Definition: MOM_diag_mediator.F90:706
mom_diag_mediator::diag_copy_storage_to_diag
subroutine, public diag_copy_storage_to_diag(diag, grid_storage)
Copy from the stored diagnostic arrays to the main diagnostic grids.
Definition: MOM_diag_mediator.F90:3530
mom_diag_mediator::id_clock_diag_remap
integer id_clock_diag_remap
Sets up diagnostics axes.
Definition: MOM_diag_mediator.F90:338
mom_diag_mediator::diag_grid_storage_init
subroutine, public diag_grid_storage_init(grid_storage, G, diag)
Allocates fields necessary to store diagnostic remapping fields.
Definition: MOM_diag_mediator.F90:3486
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_diag_mediator::set_axes_info_dsamp
subroutine set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_native)
Definition: MOM_diag_mediator.F90:533