MOM_diag_mediator.F90

! This file is part of MOM6, the Modular Ocean Model version 6.

! See the LICENSE file for licensing information.

! SPDX-License-Identifier: Apache-2.0

!> The subroutines here provide convenient wrappers to the FMS diag_manager

!! interfaces with additional diagnostic capabilities.

module mom_diag_mediator

use mom_checksums,        only : chksum0, zchksum, hchksum, uchksum, vchksum, bchksum

use mom_coms,             only : pe_here

use mom_cpu_clock,        only : cpu_clock_id, cpu_clock_begin, cpu_clock_end

use mom_cpu_clock,        only : clock_module, clock_routine

use mom_diag_buffers,     only : diag_buffer_2d, diag_buffer_3d

use mom_diag_manager_infra, only : mom_diag_manager_init, mom_diag_manager_end

use mom_diag_manager_infra, only : diag_axis_init=>mom_diag_axis_init, get_mom_diag_axis_name

use mom_diag_manager_infra, only : send_data_infra, mom_diag_field_add_attribute, east, north

use mom_diag_manager_infra, only : register_diag_field_infra, register_static_field_infra

use mom_diag_manager_infra, only : get_mom_diag_field_id, diag_field_not_found

use mom_diag_remap,       only : diag_remap_ctrl, diag_remap_update, diag_remap_calc_hmask

use mom_diag_remap,       only : diag_remap_init, diag_remap_end, diag_remap_do_remap

use mom_diag_remap,       only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field

use mom_diag_remap,       only : horizontally_average_diag_field, diag_remap_get_axes_info

use mom_diag_remap,       only : diag_remap_configure_axes, diag_remap_axes_configured

use mom_diag_remap,       only : diag_remap_diag_registration_closed, diag_remap_set_active

use mom_eos,              only : eos_type

use mom_error_handler,    only : mom_error, fatal, warning, is_root_pe, assert, calltree_showquery

use mom_error_handler,    only : calltree_enter, calltree_leave, calltree_waypoint

use mom_file_parser,      only : get_param, log_version, param_file_type

use mom_grid,             only : ocean_grid_type

use mom_interface_heights, only : thickness_to_dz

use mom_io,               only : vardesc, query_vardesc

use mom_io,               only : get_filename_appendix

use mom_safe_alloc,       only : safe_alloc_ptr, safe_alloc_alloc

use mom_string_functions, only : lowercase, slasher, ints_to_string, trim_trailing_commas

use mom_time_manager,     only : time_type, get_time

use mom_unit_scaling,     only : unit_scale_type

use mom_variables,        only : thermo_var_ptrs

use mom_verticalgrid,     only : verticalgrid_type

implicit none ; private

#undef __DO_SAFETY_CHECKS__

#define IMPLIES(A, B) ((.not. (A)) .or. (B))

#define MAX_DSAMP_LEV 2

public set_axes_info, post_data, register_diag_field, time_type

public post_data_3d_by_column, post_data_3d_final

public post_product_u, post_product_sum_u, post_product_v, post_product_sum_v

public set_masks_for_axes

! post_data_1d_k is a deprecated interface that can be replaced by a call to post_data, but

! it is being retained for backward compatibility to older versions of the ocean_BGC code.

public post_data_1d_k

public safe_alloc_ptr, safe_alloc_alloc

public enable_averaging, enable_averages, disable_averaging, query_averaging_enabled

public diag_mediator_init, diag_mediator_end, set_diag_mediator_grid

public diag_mediator_infrastructure_init

public diag_mediator_close_registration, get_diag_time_end

public diag_axis_init, ocean_register_diag, register_static_field

public register_scalar_field

public define_axes_group, diag_masks_set

public set_piecemeal_extents

public diag_register_area_ids

public register_cell_measure, diag_associate_volume_cell_measure

public diag_get_volume_cell_measure_dm_id

public diag_set_state_ptrs, diag_update_remap_grids

public diag_grid_storage_init, diag_grid_storage_end

public diag_copy_diag_to_storage, diag_copy_storage_to_diag

public diag_save_grids, diag_restore_grids

public found_in_diagtable

!> Make a diagnostic available for averaging or output.

interface post_data

  module procedure post_data_3d, post_data_2d, post_data_1d_k, post_data_0d

end interface post_data

!> Down sample a field

interface downsample_field

  module procedure downsample_field_2d, downsample_field_3d

end interface downsample_field

!> Down sample the mask of a field

interface downsample_mask

  module procedure downsample_mask_2d, downsample_mask_3d

end interface downsample_mask

!> Down sample a diagnostic field

interface downsample_diag_field

  module procedure downsample_diag_field_2d, downsample_diag_field_3d

end interface downsample_diag_field

!> Contained for down sampled masks

type, private :: diag_dsamp

  real, pointer, dimension(:,:)   :: mask2d => null() !< Mask for 2d (x-y) axes [nondim]

  real, pointer, dimension(:,:,:) :: mask3d => null() !< Mask for 3d axes [nondim]

end type diag_dsamp

!> A group of 1D axes that comprise a 1D/2D/3D mesh

type, public :: axes_grp

  character(len=15) :: id   !< The id string for this particular combination of handles.

  integer           :: rank !< Number of dimensions in the list of axes.

  integer, dimension(:), allocatable :: handles !< Handles to 1D axes.

  type(diag_ctrl), pointer :: diag_cs => null() !< Circular link back to the main diagnostics control structure

                                                !! (Used to avoid passing said structure into every possible call).

  ! ID's for cell_methods

  character(len=9) :: x_cell_method = '' !< Default nature of data representation, if axes group

                                         !! includes x-direction.

  character(len=9) :: y_cell_method = '' !< Default nature of data representation, if axes group

                                         !! includes y-direction.

  character(len=9) :: v_cell_method = '' !< Default nature of data representation, if axes group

                                         !! includes vertical direction.

  ! For remapping

  integer :: nz = 0 !< Vertical dimension of diagnostic

  integer :: vertical_coordinate_number = 0 !< Index of the corresponding diag_remap_ctrl for this axis group

  ! For detecting position on the grid

  logical :: is_h_point = .false. !< If true, indicates that this axes group is for an h-point located field.

  logical :: is_q_point = .false. !< If true, indicates that this axes group is for a q-point located field.

  logical :: is_u_point = .false. !< If true, indicates that this axes group is for a u-point located field.

  logical :: is_v_point = .false. !< If true, indicates that this axes group is for a v-point located field.

  logical :: is_layer = .false. !< If true, indicates that this axes group is for a layer vertically-located field.

  logical :: is_interface = .false. !< If true, indicates that this axes group is for an interface

                                    !! vertically-located field.

  logical :: is_native = .true. !< If true, indicates that this axes group is for a native model grid.

                                !! False for any other grid. Used for rank>2.

  logical :: needs_remapping = .false. !< If true, indicates that this axes group is for a intensive layer-located

                                       !! field that must be remapped to these axes. Used for rank>2.

  logical :: needs_interpolating = .false. !< If true, indicates that this axes group is for a sampled

                                         !! interface-located field that must be interpolated to

                                         !! these axes. Used for rank>2.

  integer :: downsample_level = 1 !< If greater than 1, the factor by which this diagnostic/axes/masks be downsampled

  ! For horizontally averaged diagnostics (applies to 2d and 3d fields only)

  type(axes_grp), pointer :: xyave_axes => null() !< The associated 1d axes for horizontally area-averaged diagnostics

  ! ID's for cell_measures

  integer :: id_area = -1 !< The diag_manager id for area to be used for cell_measure of variables with this axes_grp.

  integer :: id_volume = -1 !< The diag_manager id for volume to be used for cell_measure of variables

                            !! with this axes_grp.

  ! For masking

  real, pointer, dimension(:,:)   :: mask2d => null() !< Mask for 2d (x-y) axes [nondim]

  real, pointer, dimension(:,:,:) :: mask3d => null() !< Mask for 3d axes [nondim]

  type(diag_dsamp), dimension(2:MAX_DSAMP_LEV) :: dsamp !< Downsample container

  ! For diagnostics posted piecemeal

  type(diag_buffer_2d) :: piecemeal_2d !< A dynamically reallocated buffer for 2d piecemeal diagnostics

  type(diag_buffer_3d) :: piecemeal_3d !< A dynamically reallocated buffer for 3d piecemeal diagnostics

end type axes_grp

!> Contains an array to store a diagnostic target grid

type, private :: diag_grids_type

  real, dimension(:,:,:), allocatable :: h  !< Target grid for remapped coordinate [H ~> m or kg m-2] or [Z ~> m]

end type diag_grids_type

!> Stores all the remapping grids and the model's native space thicknesses

type, public :: diag_grid_storage

  integer                                          :: num_diag_coords !< Number of target coordinates

  real, dimension(:,:,:), allocatable              :: h_state         !< Layer thicknesses in native

                                                                      !! space [H ~> m or kg m-2]

  type(diag_grids_type), dimension(:), allocatable :: diag_grids      !< Primarily empty, except h field

end type diag_grid_storage

! Integers to encode the total cell methods

!integer :: PPP=111  ! x:point,y:point,z:point, this kind of diagnostic is not currently present in diag_table.MOM6

!integer :: PPS=112  ! x:point,y:point,z:sum  , this kind of diagnostic is not currently present in diag_table.MOM6

!integer :: PPM=113  ! x:point,y:point,z:mean , this kind of diagnostic is not currently present in diag_table.MOM6

integer :: psp=121  !< x:point,y:sum,z:point

integer :: pss=122  !< x:point,y:sum,z:point

integer :: psm=123  !< x:point,y:sum,z:mean

integer :: pmp=131  !< x:point,y:mean,z:point

integer :: pmm=133  !< x:point,y:mean,z:mean

integer :: spp=211  !< x:sum,y:point,z:point

integer :: sps=212  !< x:sum,y:point,z:sum

integer :: ssp=221  !< x:sum,y:sum,z:point

integer :: mpp=311  !< x:mean,y:point,z:point

integer :: mpm=313  !< x:mean,y:point,z:mean

integer :: mmp=331  !< x:mean,y:mean,z:point

integer :: mms=332  !< x:mean,y:mean,z:sum

integer :: sss=222  !< x:sum,y:sum,z:sum

integer :: mmm=333  !< x:mean,y:mean,z:mean

integer :: msk=-1   !< Use the downsample method of a mask

!> This type is used to represent a diagnostic at the diag_mediator level.

!!

!! There can be both 'primary' and 'secondary' diagnostics. The primaries

!! reside in the diag_cs%diags array. They have an id which is an index

!! into this array. The secondaries are 'variations' on the primary diagnostic.

!! For example the CMOR diagnostics are secondary. The secondary diagnostics

!! are kept in a list with the primary diagnostic as the head.

type, private :: diag_type

  logical :: in_use !< True if this entry is being used.

  integer :: fms_diag_id !< Underlying FMS diag_manager id.

  integer :: fms_xyave_diag_id = -1 !< For a horizontally area-averaged diagnostic.

  integer :: downsample_diag_id = -1 !< For a horizontally area-downsampled diagnostic.

  character(len=64) :: debug_str = '' !< The diagnostic name and module for FATAL errors and debugging.

  type(axes_grp), pointer :: axes => null() !< The axis group for this diagnostic

  type(diag_type), pointer :: next => null() !< Pointer to the next diagnostic

  real :: conversion_factor = 0. !< If non-zero, a factor to multiply data by before posting to FMS,

                                 !! often including factors to undo internal scaling in units of [a A-1 ~> 1]

  logical :: v_extensive = .false. !< True for vertically extensive fields (vertically integrated).

                                   !! False for intensive (concentrations).

  integer :: xyz_method = 0 !< A 3 digit integer encoding the diagnostics cell method

                            !! It can be used to determine the downsample algorithm

end type diag_type

!> Container for down sampling information

type diagcs_dsamp

  integer :: isc !< The start i-index of cell centers within the computational domain

  integer :: iec !< The end i-index of cell centers within the computational domain

  integer :: jsc !< The start j-index of cell centers within the computational domain

  integer :: jec !< The end j-index of cell centers within the computational domain

  integer :: isd !< The start i-index of cell centers within the data domain

  integer :: ied !< The end i-index of cell centers within the data domain

  integer :: jsd !< The start j-index of cell centers within the data domain

  integer :: jed !< The end j-index of cell centers within the data domain

  integer :: isg !< The start i-index of cell centers within the global domain

  integer :: ieg !< The end i-index of cell centers within the global domain

  integer :: jsg !< The start j-index of cell centers within the global domain

  integer :: jeg !< The end j-index of cell centers within the global domain

  integer :: isgb !< The start i-index of cell corners within the global domain

  integer :: iegb !< The end i-index of cell corners within the global domain

  integer :: jsgb !< The start j-index of cell corners within the global domain

  integer :: jegb !< The end j-index of cell corners within the global domain

  !>@{ Axes for each location on a diagnostic grid

  type(axes_grp)  :: axesbl, axestl, axescul, axescvl

  type(axes_grp)  :: axesbi, axesti, axescui, axescvi

  type(axes_grp)  :: axesb1, axest1, axescu1, axescv1

  type(axes_grp), dimension(:), allocatable :: remap_axestl, remap_axesbl, remap_axescul, remap_axescvl

  type(axes_grp), dimension(:), allocatable :: remap_axesti, remap_axesbi, remap_axescui, remap_axescvi

  !>@}

  real, dimension(:,:),   pointer :: mask2dt   => null() !< 2D mask array for cell-center points [nondim]

  real, dimension(:,:),   pointer :: mask2dbu  => null() !< 2D mask array for cell-corner points [nondim]

  real, dimension(:,:),   pointer :: mask2dcu  => null() !< 2D mask array for east-face points [nondim]

  real, dimension(:,:),   pointer :: mask2dcv  => null() !< 2D mask array for north-face points [nondim]

  !>@{ 3D mask arrays for diagnostics at layers (mask...L) and interfaces (mask...i), all [nondim]

  real, dimension(:,:,:), pointer :: mask3dtl  => null()

  real, dimension(:,:,:), pointer :: mask3dbl  => null()

  real, dimension(:,:,:), pointer :: mask3dcul => null()

  real, dimension(:,:,:), pointer :: mask3dcvl => null()

  real, dimension(:,:,:), pointer :: mask3dti  => null()

  real, dimension(:,:,:), pointer :: mask3dbi  => null()

  real, dimension(:,:,:), pointer :: mask3dcui => null()

  real, dimension(:,:,:), pointer :: mask3dcvi => null()

  !>@}

end type diagcs_dsamp

!> The following data type a list of diagnostic fields an their variants,

!! as well as variables that control the handling of model output.

type, public :: diag_ctrl

  integer :: available_diag_doc_unit = -1 !< The unit number of a diagnostic documentation file.

                                          !! This file is open if available_diag_doc_unit is > 0.

  integer :: chksum_iounit = -1           !< The unit number of a diagnostic documentation file.

                                          !! This file is open if available_diag_doc_unit is > 0.

  logical :: diag_as_chksum  !< If true, log chksums in a text file instead of posting diagnostics

  logical :: show_call_tree  !< Display the call tree while running. Set by VERBOSITY level.

  logical :: index_space_axes !< If true, diagnostic horizontal coordinates axes are in index space.

  ! The following fields are used for the output of the data.

  ! These give the computational-domain sizes, and are relative to a start value

  ! of 1 in memory for the tracer-point arrays.

  integer :: is  !< The start i-index of cell centers within the computational domain

  integer :: ie  !< The end i-index of cell centers within the computational domain

  integer :: js  !< The start j-index of cell centers within the computational domain

  integer :: je  !< The end j-index of cell centers within the computational domain

  ! These give the memory-domain sizes, and can be start at any value on each PE.

  integer :: isd !< The start i-index of cell centers within the data domain

  integer :: ied !< The end i-index of cell centers within the data domain

  integer :: jsd !< The start j-index of cell centers within the data domain

  integer :: jed !< The end j-index of cell centers within the data domain

  real :: time_int              !< The time interval for any fields

                                !! that are offered for averaging [s].

  type(time_type) :: time_end   !< The end time of the valid interval for any offered field.

  logical :: ave_enabled = .false. !< True if averaging is enabled.

  !>@{ The following are 3D and 2D axis groups defined for output.  The names

  !! indicate the horizontal (B, T, Cu, or Cv) and vertical (L, i, or 1) locations.

  type(axes_grp) :: axesbl, axestl, axescul, axescvl

  type(axes_grp) :: axesbi, axesti, axescui, axescvi

  type(axes_grp) :: axesb1, axest1, axescu1, axescv1

  !>@}

  type(axes_grp) :: axeszi !< A 1-D z-space axis at interfaces

  type(axes_grp) :: axeszl !< A 1-D z-space axis at layer centers

  type(axes_grp) :: axesnull !< An axis group for scalars

  ! Mask arrays for 2D diagnostics

  real, dimension(:,:),   pointer :: mask2dt   => null() !< 2D mask array for cell-center points [nondim]

  real, dimension(:,:),   pointer :: mask2dbu  => null() !< 2D mask array for cell-corner points [nondim]

  real, dimension(:,:),   pointer :: mask2dcu  => null() !< 2D mask array for east-face points [nondim]

  real, dimension(:,:),   pointer :: mask2dcv  => null() !< 2D mask array for north-face points [nondim]

  !>@{ 3D mask arrays for diagnostics at layers (mask...L) and interfaces (mask...i) all [nondim]

  real, dimension(:,:,:), pointer :: mask3dtl  => null()

  real, dimension(:,:,:), pointer :: mask3dbl  => null()

  real, dimension(:,:,:), pointer :: mask3dcul => null()

  real, dimension(:,:,:), pointer :: mask3dcvl => null()

  real, dimension(:,:,:), pointer :: mask3dti  => null()

  real, dimension(:,:,:), pointer :: mask3dbi  => null()

  real, dimension(:,:,:), pointer :: mask3dcui => null()

  real, dimension(:,:,:), pointer :: mask3dcvi => null()

  type(diagcs_dsamp), dimension(2:MAX_DSAMP_LEV) :: dsamp !< Downsample control container

  !>@}

! Space for diagnostics is dynamically allocated as it is needed.

! The chunk size is how much the array should grow on each new allocation.

#define DIAG_ALLOC_CHUNK_SIZE 100

  type(diag_type), dimension(:), allocatable :: diags !< The list of diagnostics

  integer :: next_free_diag_id !< The next unused diagnostic ID

  !> default missing value to be sent to ALL diagnostics registrations [various]

  real :: missing_value = -1.0e34

  !> Number of diagnostic vertical coordinates (remapped)

  integer :: num_diag_coords

  !> Control structure for each possible coordinate

  type(diag_remap_ctrl), dimension(:), allocatable :: diag_remap_cs

  type(diag_grid_storage) :: diag_grid_temp !< Stores the remapped diagnostic grid

  logical :: diag_grid_overridden = .false. !< True if the diagnostic grids have been overriden

  type(axes_grp), dimension(:), allocatable :: &

    remap_axeszl, &  !< The 1-D z-space cell-centered axis for remapping

    remap_axeszi     !< The 1-D z-space interface axis for remapping

  !>@{ Axes used for remapping

  type(axes_grp), dimension(:), allocatable :: remap_axestl, remap_axesbl, remap_axescul, remap_axescvl

  type(axes_grp), dimension(:), allocatable :: remap_axesti, remap_axesbi, remap_axescui, remap_axescvi

  !>@}

  ! Pointer to H, G and T&S needed for remapping

  real, dimension(:,:,:), pointer :: h => null() !< The thicknesses needed for remapping [H ~> m or kg m-2]

  real, dimension(:,:,:), pointer :: t => null() !< The temperatures needed for remapping [C ~> degC]

  real, dimension(:,:,:), pointer :: s => null() !< The salinities needed for remapping [S ~> ppt]

  type(eos_type),         pointer :: eqn_of_state => null() !< The equation of state type

  type(thermo_var_ptrs),  pointer :: tv => null()   !< A structure with thermodynamic variables that are

                                                    !! are used to convert thicknesses to vertical extents

  type(ocean_grid_type), pointer :: g => null()  !< The ocean grid type

  type(verticalgrid_type), pointer :: gv => null()  !< The model's vertical ocean grid

  type(unit_scale_type), pointer :: us => null() !< A dimensional unit scaling type

  !> The volume cell measure (special diagnostic) manager id

  integer :: volume_cell_measure_dm_id = -1

#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)

  ! Keep a copy of h so that we know whether it has changed [H ~> m or kg m-2]. If it has then

  ! need the target grid for vertical remapping needs to have been updated.

  real, dimension(:,:,:), allocatable :: h_old

#endif

  !> Number of checksum-only diagnostics

  integer :: num_chksum_diags

  real, dimension(:,:,:), allocatable :: h_begin !< Layer thicknesses at the beginning of the timestep used

                                                 !! for remapping of extensive variables [H ~> m or kg m-2]

  real, dimension(:,:,:), allocatable :: dz_begin !< Layer vertical extents at the beginning of the timestep used

                                                 !! for remapping of extensive variables [Z ~> m]

end type diag_ctrl

!>@{ CPU clocks

integer :: id_clock_diag_mediator, id_clock_diag_remap, id_clock_diag_grid_updates

!>@}

contains

!> Sets up diagnostics axes

subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical)

  type(ocean_grid_type),   intent(inout) :: g  !< Ocean grid structure

  type(verticalgrid_type), intent(in)    :: gv !< ocean vertical grid structure

  type(unit_scale_type),   intent(in)    :: us !< A dimensional unit scaling type

  type(param_file_type),   intent(in)    :: param_file !< Parameter file structure

  type(diag_ctrl),         intent(inout) :: diag_cs !< Diagnostics control structure

  logical,       optional, intent(in)    :: set_vertical !< If true or missing, set up

                                                       !! vertical axes

  ! Local variables

  integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, id_null

  integer :: id_zl_native, id_zi_native

  integer :: i, j, nz

  real :: zlev(gv%ke)     ! Numerical values for layer vertical coordinates, in unscaled units

                          ! that might be [m], [kg m-3] or [nondim], depending on the coordinate.

  real :: zinter(gv%ke+1) ! Numerical values for interface vertical coordinates, in unscaled units

                          ! that might be [m], [kg m-3] or [nondim], depending on the coordinate.

  logical :: set_vert

  real, allocatable, dimension(:) :: iaxb, iax ! Index-based integer and half-integer i-axis labels [nondim]

  real, allocatable, dimension(:) :: jaxb, jax ! Index-based integer and half-integer j-axis labels [nondim]

  set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical

  if (diag_cs%index_space_axes) then

    allocate(iaxb(g%IsgB:g%IegB))

    do i=g%IsgB, g%IegB

      iaxb(i)=real(i)

    enddo

    allocate(iax(g%isg:g%ieg))

    do i=g%isg, g%ieg

      iax(i)=real(i)-0.5

    enddo

    allocate(jaxb(g%JsgB:g%JegB))

    do j=g%JsgB, g%JegB

      jaxb(j)=real(j)

    enddo

    allocate(jax(g%jsg:g%jeg))

    do j=g%jsg, g%jeg

      jax(j)=real(j)-0.5

    enddo

  endif

  ! Horizontal axes for the native grids

  if (g%symmetric) then

    if (diag_cs%index_space_axes) then

      id_xq = diag_axis_init('iq', iaxb(g%isgB:g%iegB), 'none', 'x', &

          'q point grid-space longitude', g%Domain, position=east)

      id_yq = diag_axis_init('jq', jaxb(g%jsgB:g%jegB), 'none', 'y', &

          'q point grid space latitude', g%Domain, position=north)

    else

      id_xq = diag_axis_init('xq', g%gridLonB(g%isgB:g%iegB), g%x_axis_units, 'x', &

          'q point nominal longitude', g%Domain, position=east)

      id_yq = diag_axis_init('yq', g%gridLatB(g%jsgB:g%jegB), g%y_axis_units, 'y', &

          'q point nominal latitude', g%Domain, position=north)

    endif

  else

    if (diag_cs%index_space_axes) then

      id_xq = diag_axis_init('Iq', iaxb(g%isg:g%ieg), 'none', 'x', &

          'q point grid-space longitude', g%Domain, position=east)

      id_yq = diag_axis_init('Jq', jaxb(g%jsg:g%jeg), 'none', 'y', &

          'q point grid space latitude', g%Domain, position=north)

    else

      id_xq = diag_axis_init('xq', g%gridLonB(g%isg:g%ieg), g%x_axis_units, 'x', &

          'q point nominal longitude', g%Domain, position=east)

      id_yq = diag_axis_init('yq', g%gridLatB(g%jsg:g%jeg), g%y_axis_units, 'y', &

          'q point nominal latitude', g%Domain, position=north)

    endif

  endif

  if (diag_cs%index_space_axes) then

    id_xh = diag_axis_init('ih', iax(g%isg:g%ieg), 'none', 'x', &

        'h point grid-space longitude', g%Domain)

    id_yh = diag_axis_init('jh', jax(g%jsg:g%jeg), 'none', 'y', &

        'h point grid space latitude', g%Domain)

  else

    id_xh = diag_axis_init('xh', g%gridLonT(g%isg:g%ieg), g%x_axis_units, 'x', &

        'h point nominal longitude', g%Domain)

    id_yh = diag_axis_init('yh', g%gridLatT(g%jsg:g%jeg), g%y_axis_units, 'y', &

        'h point nominal latitude', g%Domain)

  endif

  if (set_vert) then

    nz = gv%ke

    zinter(1:nz+1) = gv%sInterface(1:nz+1)

    zlev(1:nz) = gv%sLayer(1:nz)

    id_zl = diag_axis_init('zl', zlev, trim(gv%zAxisUnits), 'z', &

                           'Layer '//trim(gv%zAxisLongName), direction=gv%direction)

    id_zi = diag_axis_init('zi', zinter, trim(gv%zAxisUnits), 'z', &

                           'Interface '//trim(gv%zAxisLongName), direction=gv%direction)

  else

    id_zl = -1 ; id_zi = -1

  endif

  id_zl_native = id_zl ; id_zi_native = id_zi

  ! Vertical axes for the interfaces and layers

  call define_axes_group(diag_cs, (/ id_zi /), diag_cs%axesZi, &

       v_cell_method='point', is_interface=.true.)

  call define_axes_group(diag_cs, (/ id_zl /), diag_cs%axesZL, &

       v_cell_method='mean', is_layer=.true.)

  ! Axis groupings for the model layers

  call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zl /), diag_cs%axesTL, &

       x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &

       is_h_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)

  call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zl /), diag_cs%axesBL, &

       x_cell_method='point', y_cell_method='point', v_cell_method='mean', &

       is_q_point=.true., is_layer=.true.)

  call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zl /), diag_cs%axesCuL, &

       x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &

       is_u_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)

  call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zl /), diag_cs%axesCvL, &

       x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &

       is_v_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)

  ! Axis groupings for the model interfaces

  call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%axesTi, &

       x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &

       is_h_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)

  call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%axesBi, &

       x_cell_method='point', y_cell_method='point', v_cell_method='point', &

       is_q_point=.true., is_interface=.true.)

  call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%axesCui, &

       x_cell_method='point', y_cell_method='mean', v_cell_method='point', &

       is_u_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)

  call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%axesCvi, &

       x_cell_method='mean', y_cell_method='point', v_cell_method='point', &

       is_v_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)

  ! Axis groupings for 2-D arrays

  call define_axes_group(diag_cs, (/ id_xh, id_yh /), diag_cs%axesT1, &

       x_cell_method='mean', y_cell_method='mean', is_h_point=.true.)

  call define_axes_group(diag_cs, (/ id_xq, id_yq /), diag_cs%axesB1, &

       x_cell_method='point', y_cell_method='point', is_q_point=.true.)

  call define_axes_group(diag_cs, (/ id_xq, id_yh /), diag_cs%axesCu1, &

       x_cell_method='point', y_cell_method='mean', is_u_point=.true.)

  call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, &

       x_cell_method='mean', y_cell_method='point', is_v_point=.true.)

  ! Define array extents for all piecemeal buffers

  call set_piecemeal_extents(diag_cs)

  ! Axis group for special null axis for scalars from diag manager.

  id_null = diag_axis_init('scalar_axis', (/0./), 'none', 'N', 'none', null_axis=.true.)

  call define_axes_group(diag_cs, (/ id_null /), diag_cs%axesNull)

  ! Set axis groups for non-native, non-downsampled grids

  if (diag_cs%num_diag_coords>0) then

    allocate(diag_cs%remap_axesZL(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesTL(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesBL(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesCuL(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesCvL(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesZi(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesTi(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesBi(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesCui(diag_cs%num_diag_coords))

    allocate(diag_cs%remap_axesCvi(diag_cs%num_diag_coords))

  endif

  do i=1, diag_cs%num_diag_coords

    ! For each possible diagnostic coordinate

    call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), g, gv, us, param_file)

    ! Allocate these arrays since the size of the diagnostic array is now known

    allocate(diag_cs%diag_remap_cs(i)%h(g%isd:g%ied,g%jsd:g%jed, diag_cs%diag_remap_cs(i)%nz))

    allocate(diag_cs%diag_remap_cs(i)%h_extensive(g%isd:g%ied,g%jsd:g%jed, diag_cs%diag_remap_cs(i)%nz))

    ! This vertical coordinate has been configured so can be used.

    if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then

      ! This fetches the 1D-axis id for layers and interfaces and overwrite

      ! id_zl and id_zi from above. It also returns the number of layers.

      call diag_remap_get_axes_info(diag_cs%diag_remap_cs(i), nz, id_zl, id_zi)

      ! Axes for z layers

      call define_axes_group(diag_cs, (/ id_zl /), diag_cs%remap_axesZL(i), &

           nz=nz, vertical_coordinate_number=i, &

           v_cell_method='mean', &

           is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true.)

      call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zl /), diag_cs%remap_axesTL(i), &

           nz=nz, vertical_coordinate_number=i, &

           x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &

           is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &

           xyave_axes=diag_cs%remap_axesZL(i))

       !! \note Remapping for B points is not yet implemented so needs_remapping is not

       !! provided for remap_axesBL

      call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zl /), diag_cs%remap_axesBL(i), &

           nz=nz, vertical_coordinate_number=i, &

           x_cell_method='point', y_cell_method='point', v_cell_method='mean', &

           is_q_point=.true., is_layer=.true., is_native=.false.)

      call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zl /), diag_cs%remap_axesCuL(i), &

           nz=nz, vertical_coordinate_number=i, &

           x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &

           is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &

           xyave_axes=diag_cs%remap_axesZL(i))

      call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zl /), diag_cs%remap_axesCvL(i), &

           nz=nz, vertical_coordinate_number=i, &

           x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &

           is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &

           xyave_axes=diag_cs%remap_axesZL(i))

      ! Axes for z interfaces

      call define_axes_group(diag_cs, (/ id_zi /), diag_cs%remap_axesZi(i), &

           nz=nz, vertical_coordinate_number=i, &

           v_cell_method='point', &

           is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true.)

      call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%remap_axesTi(i), &

           nz=nz, vertical_coordinate_number=i, &

           x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &

           is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true., &

           xyave_axes=diag_cs%remap_axesZi(i))

      !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBi

      call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%remap_axesBi(i), &

           nz=nz, vertical_coordinate_number=i, &

           x_cell_method='point', y_cell_method='point', v_cell_method='point', &

           is_q_point=.true., is_interface=.true., is_native=.false.)

      call define_axes_group(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%remap_axesCui(i), &

           nz=nz, vertical_coordinate_number=i, &

           x_cell_method='point', y_cell_method='mean', v_cell_method='point', &

           is_u_point=.true., is_interface=.true., is_native=.false., &

           needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))

      call define_axes_group(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%remap_axesCvi(i), &

           nz=nz, vertical_coordinate_number=i, &

           x_cell_method='mean', y_cell_method='point', v_cell_method='point', &

           is_v_point=.true., is_interface=.true., is_native=.false., &

           needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))

    endif

  enddo

  if (diag_cs%index_space_axes) then

    deallocate(iaxb, iax, jaxb, jax)

  endif

  ! Define the downsampled axes

  call set_axes_info_dsamp(g, gv, param_file, diag_cs, id_zl_native, id_zi_native)

  call diag_grid_storage_init(diag_cs%diag_grid_temp, g, gv, diag_cs)

end subroutine set_axes_info

subroutine set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_native)

  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure

  type(verticalgrid_type), intent(in)  :: GV !< ocean vertical grid structure

  type(param_file_type), intent(in)    :: param_file !< Parameter file structure

  type(diag_ctrl),       intent(inout) :: diag_cs !< Diagnostics control structure

  integer,               intent(in)    :: id_zl_native !< ID of native layers

  integer,               intent(in)    :: id_zi_native !< ID of native interfaces

  ! Local variables

  integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh

  integer :: i, j, nz, dl

  real, dimension(:), pointer :: gridLonT_dsamp =>null() ! The longitude of downsampled T points for labeling

                                                         ! the output axes, often in units of [degrees_N] or

                                                         ! [km] or [m] or [gridpoints].

  real, dimension(:), pointer :: gridLatT_dsamp =>null() ! The latitude of downsampled T points for labeling

                                                         ! the output axes, often in units of [degrees_N] or

                                                         ! [km] or [m] or [gridpoints].

  real, dimension(:), pointer :: gridLonB_dsamp =>null() ! The longitude of downsampled B points for labeling

                                                         ! the output axes, often in units of [degrees_N] or

                                                         ! [km] or [m] or [gridpoints].

  real, dimension(:), pointer :: gridLatB_dsamp =>null() ! The latitude of downsampled B points for labeling

                                                         ! the output axes, often in units of [degrees_N] or

                                                         ! [km] or [m] or [gridpoints].

  id_zl = id_zl_native ; id_zi = id_zi_native

  ! Axes group for native downsampled diagnostics

  do dl=2,max_dsamp_lev

    if (dl /= 2) call mom_error(fatal, "set_axes_info_dsamp: Downsample level other than 2 is not supported yet!")

    if (g%symmetric) then

      allocate(gridlonb_dsamp(diag_cs%dsamp(dl)%isgB:diag_cs%dsamp(dl)%iegB))

      allocate(gridlatb_dsamp(diag_cs%dsamp(dl)%jsgB:diag_cs%dsamp(dl)%jegB))

      do i=diag_cs%dsamp(dl)%isgB,diag_cs%dsamp(dl)%iegB ; gridlonb_dsamp(i) = g%gridLonB(g%isgB+dl*i) ; enddo

      do j=diag_cs%dsamp(dl)%jsgB,diag_cs%dsamp(dl)%jegB ; gridlatb_dsamp(j) = g%gridLatB(g%jsgB+dl*j) ; enddo

      id_xq = diag_axis_init('xq', gridlonb_dsamp, g%x_axis_units, 'x', &

            'q point nominal longitude', g%Domain, coarsen=2)

      id_yq = diag_axis_init('yq', gridlatb_dsamp, g%y_axis_units, 'y', &

            'q point nominal latitude', g%Domain, coarsen=2)

      deallocate(gridlonb_dsamp, gridlatb_dsamp)

    else

      allocate(gridlonb_dsamp(diag_cs%dsamp(dl)%isg:diag_cs%dsamp(dl)%ieg))

      allocate(gridlatb_dsamp(diag_cs%dsamp(dl)%jsg:diag_cs%dsamp(dl)%jeg))

      do i=diag_cs%dsamp(dl)%isg,diag_cs%dsamp(dl)%ieg ; gridlonb_dsamp(i) = g%gridLonB(g%isg+dl*i-2) ; enddo

      do j=diag_cs%dsamp(dl)%jsg,diag_cs%dsamp(dl)%jeg ; gridlatb_dsamp(j) = g%gridLatB(g%jsg+dl*j-2) ; enddo

      id_xq = diag_axis_init('xq', gridlonb_dsamp, g%x_axis_units, 'x', &

            'q point nominal longitude', g%Domain, coarsen=2)

      id_yq = diag_axis_init('yq', gridlatb_dsamp, g%y_axis_units, 'y', &

            'q point nominal latitude', g%Domain, coarsen=2)

      deallocate(gridlonb_dsamp, gridlatb_dsamp)

    endif

    allocate(gridlont_dsamp(diag_cs%dsamp(dl)%isg:diag_cs%dsamp(dl)%ieg))

    allocate(gridlatt_dsamp(diag_cs%dsamp(dl)%jsg:diag_cs%dsamp(dl)%jeg))

    do i=diag_cs%dsamp(dl)%isg,diag_cs%dsamp(dl)%ieg ; gridlont_dsamp(i) = g%gridLonT(g%isg+dl*i-2) ; enddo

    do j=diag_cs%dsamp(dl)%jsg,diag_cs%dsamp(dl)%jeg ; gridlatt_dsamp(j) = g%gridLatT(g%jsg+dl*j-2) ; enddo

    id_xh = diag_axis_init('xh', gridlont_dsamp, g%x_axis_units, 'x', &

          'h point nominal longitude', g%Domain, coarsen=2)

    id_yh = diag_axis_init('yh', gridlatt_dsamp, g%y_axis_units, 'y', &

          'h point nominal latitude', g%Domain, coarsen=2)

    deallocate(gridlont_dsamp, gridlatt_dsamp)

    ! Axis groupings for the model layers

    call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zl /), diag_cs%dsamp(dl)%axesTL, dl, &

            x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &

            is_h_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)

    call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zl /), diag_cs%dsamp(dl)%axesBL, dl, &

            x_cell_method='point', y_cell_method='point', v_cell_method='mean', &

            is_q_point=.true., is_layer=.true.)

    call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zl /), diag_cs%dsamp(dl)%axesCuL, dl, &

            x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &

            is_u_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)

    call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zl /), diag_cs%dsamp(dl)%axesCvL, dl, &

            x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &

            is_v_point=.true., is_layer=.true., xyave_axes=diag_cs%axesZL)

    ! Axis groupings for the model interfaces

    call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%dsamp(dl)%axesTi, dl, &

            x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &

            is_h_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)

    call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%dsamp(dl)%axesBi, dl, &

            x_cell_method='point', y_cell_method='point', v_cell_method='point', &

            is_q_point=.true., is_interface=.true.)

    call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%dsamp(dl)%axesCui, dl, &

            x_cell_method='point', y_cell_method='mean', v_cell_method='point', &

            is_u_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)

    call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%dsamp(dl)%axesCvi, dl, &

            x_cell_method='mean', y_cell_method='point', v_cell_method='point', &

            is_v_point=.true., is_interface=.true., xyave_axes=diag_cs%axesZi)

    ! Axis groupings for 2-D arrays

    call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh /), diag_cs%dsamp(dl)%axesT1, dl, &

            x_cell_method='mean', y_cell_method='mean', is_h_point=.true.)

    call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq /), diag_cs%dsamp(dl)%axesB1, dl, &

            x_cell_method='point', y_cell_method='point', is_q_point=.true.)

    call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh /), diag_cs%dsamp(dl)%axesCu1, dl, &

            x_cell_method='point', y_cell_method='mean', is_u_point=.true.)

    call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq /), diag_cs%dsamp(dl)%axesCv1, dl, &

            x_cell_method='mean', y_cell_method='point', is_v_point=.true.)

    ! Axis groupings with a non-native vertical coordinate

    if (diag_cs%num_diag_coords>0) then

      allocate(diag_cs%dsamp(dl)%remap_axesTL(diag_cs%num_diag_coords))

      allocate(diag_cs%dsamp(dl)%remap_axesBL(diag_cs%num_diag_coords))

      allocate(diag_cs%dsamp(dl)%remap_axesCuL(diag_cs%num_diag_coords))

      allocate(diag_cs%dsamp(dl)%remap_axesCvL(diag_cs%num_diag_coords))

      allocate(diag_cs%dsamp(dl)%remap_axesTi(diag_cs%num_diag_coords))

      allocate(diag_cs%dsamp(dl)%remap_axesBi(diag_cs%num_diag_coords))

      allocate(diag_cs%dsamp(dl)%remap_axesCui(diag_cs%num_diag_coords))

      allocate(diag_cs%dsamp(dl)%remap_axesCvi(diag_cs%num_diag_coords))

    endif

    do i=1, diag_cs%num_diag_coords

      ! For each possible diagnostic coordinate

      ! call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), G, GV, param_file)

      ! This vertical coordinate has been configured so can be used.

      if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then

        ! This fetches the 1D-axis id for layers and interfaces and overwrite

        ! id_zl and id_zi from above. It also returns the number of layers.

        call diag_remap_get_axes_info(diag_cs%diag_remap_cs(i), nz, id_zl, id_zi)

        ! Axes for z layers

        call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zl /), diag_cs%dsamp(dl)%remap_axesTL(i), dl, &

                nz=nz, vertical_coordinate_number=i, &

                x_cell_method='mean', y_cell_method='mean', v_cell_method='mean', &

                is_h_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &

                xyave_axes=diag_cs%remap_axesZL(i))

        !! \note Remapping for B points is not yet implemented so needs_remapping is not

        !! provided for remap_axesBL

        call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zl /), diag_cs%dsamp(dl)%remap_axesBL(i), dl, &

                nz=nz, vertical_coordinate_number=i, &

                x_cell_method='point', y_cell_method='point', v_cell_method='mean', &

                is_q_point=.true., is_layer=.true., is_native=.false.)

        call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zl /), diag_cs%dsamp(dl)%remap_axesCuL(i), dl, &

                nz=nz, vertical_coordinate_number=i, &

                x_cell_method='point', y_cell_method='mean', v_cell_method='mean', &

                is_u_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &

                xyave_axes=diag_cs%remap_axesZL(i))

        call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zl /), diag_cs%dsamp(dl)%remap_axesCvL(i), dl, &

                nz=nz, vertical_coordinate_number=i, &

                x_cell_method='mean', y_cell_method='point', v_cell_method='mean', &

                is_v_point=.true., is_layer=.true., is_native=.false., needs_remapping=.true., &

                xyave_axes=diag_cs%remap_axesZL(i))

        ! Axes for z interfaces

        call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%dsamp(dl)%remap_axesTi(i), dl, &

                nz=nz, vertical_coordinate_number=i, &

                x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &

                is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true., &

                xyave_axes=diag_cs%remap_axesZi(i))

        !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBi

        call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%dsamp(dl)%remap_axesBi(i), dl, &

                nz=nz, vertical_coordinate_number=i, &

                x_cell_method='point', y_cell_method='point', v_cell_method='point', &

                is_q_point=.true., is_interface=.true., is_native=.false.)

        call define_axes_group_dsamp(diag_cs, (/ id_xq, id_yh, id_zi /), diag_cs%dsamp(dl)%remap_axesCui(i), dl, &

                nz=nz, vertical_coordinate_number=i, &

                x_cell_method='point', y_cell_method='mean', v_cell_method='point', &

                is_u_point=.true., is_interface=.true., is_native=.false., &

                needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))

        call define_axes_group_dsamp(diag_cs, (/ id_xh, id_yq, id_zi /), diag_cs%dsamp(dl)%remap_axesCvi(i), dl, &

                nz=nz, vertical_coordinate_number=i, &

                x_cell_method='mean', y_cell_method='point', v_cell_method='point', &

                is_v_point=.true., is_interface=.true., is_native=.false., &

                needs_interpolating=.true., xyave_axes=diag_cs%remap_axesZi(i))

      endif

    enddo

  enddo

end subroutine set_axes_info_dsamp

!> set_masks_for_axes sets up the 2d and 3d masks for diagnostics using the current grid

!! recorded after calling diag_update_remap_grids()

subroutine set_masks_for_axes(G, diag_cs)

  type(ocean_grid_type), target, intent(in) :: g !< The ocean grid type.

  type(diag_ctrl),               pointer    :: diag_cs !< A pointer to a type with many variables

                                                       !! used for diagnostics

  ! Local variables

  integer :: c, nk, i, j, k

  type(axes_grp), pointer :: axes => null(), h_axes => null() ! Current axes, for convenience

  do c=1, diag_cs%num_diag_coords

    ! This vertical coordinate has been configured so can be used.

    if (diag_remap_axes_configured(diag_cs%diag_remap_cs(c))) then

      ! Level/layer h-points in diagnostic coordinate

      axes => diag_cs%remap_axesTL(c)

      nk = axes%nz

      allocate( axes%mask3d(g%isd:g%ied,g%jsd:g%jed,nk), source=0. )

      call diag_remap_calc_hmask(diag_cs%diag_remap_cs(c), g, axes%mask3d)

      h_axes => diag_cs%remap_axesTL(c) ! Use the h-point masks to generate the u-, v- and q- masks

      ! Level/layer u-points in diagnostic coordinate

      axes => diag_cs%remap_axesCuL(c)

      call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at u-layers')

      call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')

      allocate( axes%mask3d(g%IsdB:g%IedB,g%jsd:g%jed,nk), source=0. )

      do k = 1, nk ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec

        if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j,k) > 0.) axes%mask3d(i,j,k) = 1.

      enddo ; enddo ; enddo

      ! Level/layer v-points in diagnostic coordinate

      axes => diag_cs%remap_axesCvL(c)

      call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at v-layers')

      call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')

      allocate( axes%mask3d(g%isd:g%ied,g%JsdB:g%JedB,nk), source=0. )

      do k = 1, nk ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec

        if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,j,k) = 1.

      enddo ; enddo ; enddo

      ! Level/layer q-points in diagnostic coordinate

      axes => diag_cs%remap_axesBL(c)

      call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at q-layers')

      call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')

      allocate( axes%mask3d(g%IsdB:g%IedB,g%JsdB:g%JedB,nk), source=0. )

      do k = 1, nk ; do j=g%jsc-1,g%jec ; do i=g%isc-1,g%iec

        if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j+1,k) + &

            h_axes%mask3d(i+1,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,j,k) = 1.

      enddo ; enddo ; enddo

      ! Interface h-points in diagnostic coordinate (w-point)

      axes => diag_cs%remap_axesTi(c)

      call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at h-interfaces')

      call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')

      allocate( axes%mask3d(g%isd:g%ied,g%jsd:g%jed,nk+1), source=0. )

      do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1

        if (h_axes%mask3d(i,j,1) > 0.) axes%mask3d(i,j,1) = 1.

        do k = 2, nk

          if (h_axes%mask3d(i,j,k-1) + h_axes%mask3d(i,j,k) > 0.) axes%mask3d(i,j,k) = 1.

        enddo

        if (h_axes%mask3d(i,j,nk) > 0.) axes%mask3d(i,j,nk+1) = 1.

      enddo ; enddo

      h_axes => diag_cs%remap_axesTi(c) ! Use the w-point masks to generate the u-, v- and q- masks

      ! Interface u-points in diagnostic coordinate

      axes => diag_cs%remap_axesCui(c)

      call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at u-interfaces')

      call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')

      allocate( axes%mask3d(g%IsdB:g%IedB,g%jsd:g%jed,nk+1), source=0. )

      do k = 1, nk+1 ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec

        if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j,k) > 0.) axes%mask3d(i,j,k) = 1.

      enddo ; enddo ; enddo

      ! Interface v-points in diagnostic coordinate

      axes => diag_cs%remap_axesCvi(c)

      call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at v-interfaces')

      call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')

      allocate( axes%mask3d(g%isd:g%ied,g%JsdB:g%JedB,nk+1), source=0. )

      do k = 1, nk+1 ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec

        if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,j,k) = 1.

      enddo ; enddo ; enddo

      ! Interface q-points in diagnostic coordinate

      axes => diag_cs%remap_axesBi(c)

      call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at q-interfaces')

      call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')

      allocate( axes%mask3d(g%IsdB:g%IedB,g%JsdB:g%JedB,nk+1), source=0. )

      do k = 1, nk ; do j=g%jsc-1,g%jec ; do i=g%isc-1,g%iec

        if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j+1,k) + &

            h_axes%mask3d(i+1,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(i,j,k) = 1.

      enddo ; enddo ; enddo

    endif

  enddo

  ! Allocate and initialize the downsampled masks for the axes

  call set_masks_for_axes_dsamp(g, diag_cs)

end subroutine set_masks_for_axes

subroutine set_masks_for_axes_dsamp(G, diag_cs)

  type(ocean_grid_type), target, intent(in) :: G !< The ocean grid type.

  type(diag_ctrl),               pointer    :: diag_cs !< A pointer to a type with many variables

                                                       !! used for diagnostics

  ! Local variables

  integer :: c, dl

  type(axes_grp), pointer :: axes => null() ! Current axes, for convenience

  ! Each downsampled axis needs both downsampled and non-downsampled masks.

  ! The downsampled mask is needed for sending out the diagnostics output via diag_manager.

  ! The non-downsampled mask is needed for downsampling the diagnostics field.

  do dl=2,max_dsamp_lev

    if (dl /= 2) call mom_error(fatal, "set_masks_for_axes_dsamp: Downsample level other than 2 is not supported!")

    do c=1, diag_cs%num_diag_coords

      ! Level/layer h-points in diagnostic coordinate

      axes => diag_cs%remap_axesTL(c)

      call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTL(c)%dsamp(dl)%mask3d, &

              dl, g%isc, g%jsc, g%isd, g%jsd, &

              g%HId2%isc, g%HId2%iec, g%HId2%jsc, g%HId2%jec, g%HId2%isd, g%HId2%ied, g%HId2%jsd, g%HId2%jed)

      diag_cs%dsamp(dl)%remap_axesTL(c)%mask3d => axes%mask3d ! Set a pointer to the non-downsampled mask

      ! Level/layer u-points in diagnostic coordinate

      axes => diag_cs%remap_axesCuL(c)

      call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCuL(c)%dsamp(dl)%mask3d, &

              dl, g%IscB, g%jsc, g%IsdB, g%jsd, &

              g%HId2%IscB, g%HId2%IecB, g%HId2%jsc, g%HId2%jec, g%HId2%IsdB, g%HId2%IedB, g%HId2%jsd, g%HId2%jed)

      diag_cs%dsamp(dl)%remap_axesCul(c)%mask3d => axes%mask3d ! Set a pointer to the non-downsampled mask

      ! Level/layer v-points in diagnostic coordinate

      axes => diag_cs%remap_axesCvL(c)

      call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvL(c)%dsamp(dl)%mask3d, &

              dl, g%isc, g%JscB, g%isd, g%JsdB, &

              g%HId2%isc, g%HId2%iec, g%HId2%JscB, g%HId2%JecB, g%HId2%isd, g%HId2%ied, g%HId2%JsdB, g%HId2%JedB)

      diag_cs%dsamp(dl)%remap_axesCvL(c)%mask3d => axes%mask3d ! Set a pointer to the non-downsampled mask

      ! Level/layer q-points in diagnostic coordinate

      axes => diag_cs%remap_axesBL(c)

      call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBL(c)%dsamp(dl)%mask3d, &

              dl, g%IscB, g%JscB, g%IsdB, g%JsdB, &

              g%HId2%IscB, g%HId2%IecB, g%HId2%JscB, g%HId2%JecB, g%HId2%IsdB, g%HId2%IedB, g%HId2%JsdB, g%HId2%JedB)

      diag_cs%dsamp(dl)%remap_axesBL(c)%mask3d => axes%mask3d ! Set a pointer to the non-downsampled mask

      ! Interface h-points in diagnostic coordinate (w-point)

      axes => diag_cs%remap_axesTi(c)

      call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTi(c)%dsamp(dl)%mask3d,  &

              dl, g%isc, g%jsc, g%isd, g%jsd, &

              g%HId2%isc, g%HId2%iec, g%HId2%jsc, g%HId2%jec, g%HId2%isd, g%HId2%ied, g%HId2%jsd, g%HId2%jed)

      diag_cs%dsamp(dl)%remap_axesTi(c)%mask3d => axes%mask3d ! Set a pointer to the non-downsampled mask

      ! Interface u-points in diagnostic coordinate

      axes => diag_cs%remap_axesCui(c)

      call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCui(c)%dsamp(dl)%mask3d,  &

              dl, g%IscB, g%jsc, g%IsdB, g%jsd, &

              g%HId2%IscB, g%HId2%IecB, g%HId2%jsc, g%HId2%jec, g%HId2%IsdB, g%HId2%IedB, g%HId2%jsd, g%HId2%jed)

      diag_cs%dsamp(dl)%remap_axesCui(c)%mask3d => axes%mask3d ! Set a pointer to the non-downsampled mask

      ! Interface v-points in diagnostic coordinate

      axes => diag_cs%remap_axesCvi(c)

      call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvi(c)%dsamp(dl)%mask3d,  &

              dl, g%isc, g%JscB, g%isd, g%JsdB, &

              g%HId2%isc, g%HId2%iec, g%HId2%JscB, g%HId2%JecB, g%HId2%isd, g%HId2%ied, g%HId2%JsdB, g%HId2%JedB)

      diag_cs%dsamp(dl)%remap_axesCvi(c)%mask3d => axes%mask3d ! Set a pointer to the non-downsampled mask

      ! Interface q-points in diagnostic coordinate

      axes => diag_cs%remap_axesBi(c)

      call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBi(c)%dsamp(dl)%mask3d,  &

              dl, g%IscB, g%JscB, g%IsdB, g%JsdB, &

              g%HId2%IscB, g%HId2%IecB, g%HId2%JscB, g%HId2%JecB, g%HId2%IsdB, g%HId2%IedB, g%HId2%JsdB, g%HId2%JedB)

      diag_cs%dsamp(dl)%remap_axesBi(c)%mask3d => axes%mask3d ! Set a pointer to the non-downsampled mask

    enddo

  enddo

end subroutine set_masks_for_axes_dsamp

!> Attaches the id of cell areas to axes groups for use with cell_measures

subroutine diag_register_area_ids(diag_cs, id_area_t, id_area_q)

  type(diag_ctrl),   intent(inout) :: diag_cs   !< Diagnostics control structure

  integer, optional, intent(in)    :: id_area_t !< Diag_mediator id for area of h-cells

  integer, optional, intent(in)    :: id_area_q !< Diag_mediator id for area of q-cells

  ! Local variables

  integer :: fms_id, i

  if (present(id_area_t)) then

    fms_id = diag_cs%diags(id_area_t)%fms_diag_id

    diag_cs%axesT1%id_area = fms_id

    diag_cs%axesTi%id_area = fms_id

    diag_cs%axesTL%id_area = fms_id

    do i=1, diag_cs%num_diag_coords

      diag_cs%remap_axesTL(i)%id_area = fms_id

      diag_cs%remap_axesTi(i)%id_area = fms_id

    enddo

  endif

  if (present(id_area_q)) then

    fms_id = diag_cs%diags(id_area_q)%fms_diag_id

    diag_cs%axesB1%id_area = fms_id

    diag_cs%axesBi%id_area = fms_id

    diag_cs%axesBL%id_area = fms_id

    do i=1, diag_cs%num_diag_coords

      diag_cs%remap_axesBL(i)%id_area = fms_id

      diag_cs%remap_axesBi(i)%id_area = fms_id

    enddo

  endif

end subroutine diag_register_area_ids

!> Sets a handle inside diagnostics mediator to associate 3d cell measures

subroutine register_cell_measure(G, diag, Time)

  type(ocean_grid_type),   intent(in)    :: g    !< Ocean grid structure

  type(diag_ctrl), target, intent(inout) :: diag !< Regulates diagnostic output

  type(time_type),         intent(in)    :: time !< Model time

  ! Local variables

  integer :: id

  id = register_diag_field('ocean_model', 'volcello', diag%axesTL, &

                           time, 'Ocean grid-cell volume', units='m3', conversion=1.0, &

                           standard_name='ocean_volume', v_extensive=.true., &

                           x_cell_method='sum', y_cell_method='sum')

  call diag_associate_volume_cell_measure(diag, id)

end subroutine register_cell_measure

!> Attaches the id of cell volumes to axes groups for use with cell_measures

subroutine diag_associate_volume_cell_measure(diag_cs, id_h_volume)

  type(diag_ctrl),   intent(inout) :: diag_cs     !< Diagnostics control structure

  integer,           intent(in)    :: id_h_volume !< Diag_manager id for volume of h-cells

  ! Local variables

  type(diag_type), pointer :: tmp => null()

  if (id_h_volume<=0) return ! Do nothing

  diag_cs%volume_cell_measure_dm_id = id_h_volume ! Record for diag_get_volume_cell_measure_dm_id()

  ! Set the cell measure for this axes group to the FMS id in this coordinate system

  diag_cs%diags(id_h_volume)%axes%id_volume = diag_cs%diags(id_h_volume)%fms_diag_id

  tmp => diag_cs%diags(id_h_volume)%next ! First item in the list, if any

  do while (associated(tmp))

    ! Set the cell measure for this axes group to the FMS id in this coordinate system

    tmp%axes%id_volume = tmp%fms_diag_id

    tmp => tmp%next ! Move to next axes group for this field

  enddo

end subroutine diag_associate_volume_cell_measure

!> Returns diag_manager id for cell measure of h-cells

integer function diag_get_volume_cell_measure_dm_id(diag_cs)

  type(diag_ctrl),   intent(in) :: diag_cs   !< Diagnostics control structure

  diag_get_volume_cell_measure_dm_id = diag_cs%volume_cell_measure_dm_id

end function diag_get_volume_cell_measure_dm_id

!> Define a group of "axes" from a list of handles and associate a mask with it

subroutine 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)

  type(diag_ctrl), target,    intent(in)  :: diag_cs !< Diagnostics control structure

  integer, dimension(:),      intent(in)  :: handles !< A list of 1D axis handles

  type(axes_grp),             intent(out) :: axes    !< The group of 1D axes

  integer,          optional, intent(in)  :: nz      !< Number of layers in this diagnostic grid

  integer,          optional, intent(in)  :: vertical_coordinate_number !< Index number for vertical coordinate

  character(len=*), optional, intent(in)  :: x_cell_method !< A x-direction cell method used to construct the

                                                           !! "cell_methods" attribute in CF convention

  character(len=*), optional, intent(in)  :: y_cell_method !< A y-direction cell method used to construct the

                                                           !! "cell_methods" attribute in CF convention

  character(len=*), optional, intent(in)  :: v_cell_method !< A vertical direction cell method used to construct

                                                        !! the "cell_methods" attribute in CF convention

  logical,          optional, intent(in)  :: is_h_point !< If true, indicates this axes group for h-point

                                                        !! located fields

  logical,          optional, intent(in)  :: is_q_point !< If true, indicates this axes group for q-point

                                                        !! located fields

  logical,          optional, intent(in)  :: is_u_point !< If true, indicates this axes group for

                                                        !! u-point located fields

  logical,          optional, intent(in)  :: is_v_point !< If true, indicates this axes group for

                                                        !! v-point located fields

  logical,          optional, intent(in)  :: is_layer   !< If true, indicates that this axes group is

                                                        !! for a layer vertically-located field.

  logical,          optional, intent(in)  :: is_interface !< If true, indicates that this axes group

                                                        !! is for an interface vertically-located field.

  logical,          optional, intent(in)  :: is_native  !< If true, indicates that this axes group is

                                                        !! for a native model grid. False for any other grid.

  logical,          optional, intent(in)  :: needs_remapping !< If true, indicates that this axes group is

                                                        !! for a intensive layer-located field that must

                                                        !! be remapped to these axes. Used for rank>2.

  logical,          optional, intent(in)  :: needs_interpolating !< If true, indicates that this axes group

                                                        !! is for a sampled interface-located field that must

                                                        !! be interpolated to these axes. Used for rank>2.

  type(axes_grp),   optional, target      :: xyave_axes !< The corresponding axes group for horizontally

                                                        !! area-average diagnostics

  ! Local variables

  integer :: n

  n = size(handles)

  if (n<1 .or. n>3) call mom_error(fatal, "define_axes_group: wrong size for list of handles!")

  allocate( axes%handles(n) )

  axes%id = ints_to_string(handles, max(n,3)) ! Identifying string

  axes%rank = n

  axes%handles(:) = handles(:)

  axes%diag_cs => diag_cs ! A (circular) link back to the diag_cs structure

  if (present(x_cell_method)) then

    if (axes%rank<2) call mom_error(fatal, 'define_axes_group: ' // &

                                           'Can not set x_cell_method for rank<2.')

    axes%x_cell_method = trim(x_cell_method)

  else

    axes%x_cell_method = ''

  endif

  if (present(y_cell_method)) then

    if (axes%rank<2) call mom_error(fatal, 'define_axes_group: ' // &

                                           'Can not set y_cell_method for rank<2.')

    axes%y_cell_method = trim(y_cell_method)

  else

    axes%y_cell_method = ''

  endif

  if (present(v_cell_method)) then

    if (axes%rank/=1 .and. axes%rank/=3) call mom_error(fatal, 'define_axes_group: ' // &

                                           'Can not set v_cell_method for rank<>1 or 3.')

    axes%v_cell_method = trim(v_cell_method)

  else

    axes%v_cell_method = ''

  endif

  if (present(nz)) axes%nz = nz

  if (present(vertical_coordinate_number)) axes%vertical_coordinate_number = vertical_coordinate_number

  if (present(is_h_point)) axes%is_h_point = is_h_point

  if (present(is_q_point)) axes%is_q_point = is_q_point

  if (present(is_u_point)) axes%is_u_point = is_u_point

  if (present(is_v_point)) axes%is_v_point = is_v_point

  if (present(is_layer)) axes%is_layer = is_layer

  if (present(is_interface)) axes%is_interface = is_interface

  if (present(is_native)) axes%is_native = is_native

  if (present(needs_remapping)) axes%needs_remapping = needs_remapping

  if (present(needs_interpolating)) axes%needs_interpolating = needs_interpolating

  if (present(xyave_axes)) axes%xyave_axes => xyave_axes

  ! Setup masks for this axes group

  axes%mask2d => null()

  if (axes%rank==2) then

    if (axes%is_h_point) axes%mask2d => diag_cs%mask2dT

    if (axes%is_u_point) axes%mask2d => diag_cs%mask2dCu

    if (axes%is_v_point) axes%mask2d => diag_cs%mask2dCv

    if (axes%is_q_point) axes%mask2d => diag_cs%mask2dBu

  endif

  ! A static 3d mask for non-native coordinates can only be setup when a grid is available

  axes%mask3d => null()

  if (axes%rank==3 .and. axes%is_native) then

    ! Native variables can/should use the native masks copied into diag_cs

    if (axes%is_layer) then

      if (axes%is_h_point) axes%mask3d => diag_cs%mask3dTL

      if (axes%is_u_point) axes%mask3d => diag_cs%mask3dCuL

      if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvL

      if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBL

    elseif (axes%is_interface) then

      if (axes%is_h_point) axes%mask3d => diag_cs%mask3dTi

      if (axes%is_u_point) axes%mask3d => diag_cs%mask3dCui

      if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvi

      if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBi

    endif

  endif

end subroutine define_axes_group

!> Defines a group of downsampled "axes" from list of handles

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)

  type(diag_ctrl), target,    intent(in)  :: diag_cs !< Diagnostics control structure

  integer, dimension(:),      intent(in)  :: handles !< A list of 1D axis handles

  type(axes_grp),             intent(out) :: axes    !< The group of 1D axes

  integer,                    intent(in)  :: dl      !< Downsample level

  integer,          optional, intent(in)  :: nz      !< Number of layers in this diagnostic grid

  integer,          optional, intent(in)  :: vertical_coordinate_number !< Index number for vertical coordinate

  character(len=*), optional, intent(in)  :: x_cell_method !< A x-direction cell method used to construct the

                                                           !! "cell_methods" attribute in CF convention

  character(len=*), optional, intent(in)  :: y_cell_method !< A y-direction cell method used to construct the

                                                           !! "cell_methods" attribute in CF convention

  character(len=*), optional, intent(in)  :: v_cell_method !< A vertical direction cell method used to construct

                                                        !! the "cell_methods" attribute in CF convention

  logical,          optional, intent(in)  :: is_h_point !< If true, indicates this axes group for h-point

                                                        !! located fields

  logical,          optional, intent(in)  :: is_q_point !< If true, indicates this axes group for q-point

                                                        !! located fields

  logical,          optional, intent(in)  :: is_u_point !< If true, indicates this axes group for

                                                        !! u-point located fields

  logical,          optional, intent(in)  :: is_v_point !< If true, indicates this axes group for

                                                        !! v-point located fields

  logical,          optional, intent(in)  :: is_layer   !< If true, indicates that this axes group is

                                                        !! for a layer vertically-located field.

  logical,          optional, intent(in)  :: is_interface !< If true, indicates that this axes group

                                                        !! is for an interface vertically-located field.

  logical,          optional, intent(in)  :: is_native  !< If true, indicates that this axes group is

                                                        !! for a native model grid. False for any other grid.

  logical,          optional, intent(in)  :: needs_remapping !< If true, indicates that this axes group is

                                                        !! for a intensive layer-located field that must

                                                        !! be remapped to these axes. Used for rank>2.

  logical,          optional, intent(in)  :: needs_interpolating !< If true, indicates that this axes group

                                                        !! is for a sampled interface-located field that must

                                                        !! be interpolated to these axes. Used for rank>2.

  type(axes_grp),   optional, target      :: xyave_axes !< The corresponding axes group for horizontally

                                                        !! area-average diagnostics

  ! Local variables

  integer :: n

  n = size(handles)

  if (n<1 .or. n>3) call mom_error(fatal, "define_axes_group: wrong size for list of handles!")

  allocate( axes%handles(n) )

  axes%id = ints_to_string(handles, max(n,3)) ! Identifying string

  axes%rank = n

  axes%handles(:) = handles(:)

  axes%diag_cs => diag_cs ! A (circular) link back to the diag_cs structure

  if (present(x_cell_method)) then

    if (axes%rank<2) call mom_error(fatal, 'define_axes_group: ' // &

                                           'Can not set x_cell_method for rank<2.')

    axes%x_cell_method = trim(x_cell_method)

  else

    axes%x_cell_method = ''

  endif

  if (present(y_cell_method)) then

    if (axes%rank<2) call mom_error(fatal, 'define_axes_group: ' // &

                                           'Can not set y_cell_method for rank<2.')

    axes%y_cell_method = trim(y_cell_method)

  else

    axes%y_cell_method = ''

  endif

  if (present(v_cell_method)) then

    if (axes%rank/=1 .and. axes%rank/=3) call mom_error(fatal, 'define_axes_group: ' // &

                                           'Can not set v_cell_method for rank<>1 or 3.')

    axes%v_cell_method = trim(v_cell_method)

  else

    axes%v_cell_method = ''

  endif

  axes%downsample_level = dl

  if (present(nz)) axes%nz = nz

  if (present(vertical_coordinate_number)) axes%vertical_coordinate_number = vertical_coordinate_number

  if (present(is_h_point)) axes%is_h_point = is_h_point

  if (present(is_q_point)) axes%is_q_point = is_q_point

  if (present(is_u_point)) axes%is_u_point = is_u_point

  if (present(is_v_point)) axes%is_v_point = is_v_point

  if (present(is_layer)) axes%is_layer = is_layer

  if (present(is_interface)) axes%is_interface = is_interface

  if (present(is_native)) axes%is_native = is_native

  if (present(needs_remapping)) axes%needs_remapping = needs_remapping

  if (present(needs_interpolating)) axes%needs_interpolating = needs_interpolating

  if (present(xyave_axes)) axes%xyave_axes => xyave_axes

  ! Setup masks for this axes group

  axes%mask2d => null()

  if (axes%rank==2) then

    if (axes%is_h_point) axes%mask2d => diag_cs%mask2dT

    if (axes%is_u_point) axes%mask2d => diag_cs%mask2dCu

    if (axes%is_v_point) axes%mask2d => diag_cs%mask2dCv

    if (axes%is_q_point) axes%mask2d => diag_cs%mask2dBu

  endif

  ! A static 3d mask for non-native coordinates can only be setup when a grid is available

  axes%mask3d => null()

  if (axes%rank==3 .and. axes%is_native) then

    ! Native variables can/should use the native masks copied into diag_cs

    if (axes%is_layer) then

      if (axes%is_h_point) axes%mask3d => diag_cs%mask3dTL

      if (axes%is_u_point) axes%mask3d => diag_cs%mask3dCuL

      if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvL

      if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBL

    elseif (axes%is_interface) then

      if (axes%is_h_point) axes%mask3d => diag_cs%mask3dTi

      if (axes%is_u_point) axes%mask3d => diag_cs%mask3dCui

      if (axes%is_v_point) axes%mask3d => diag_cs%mask3dCvi

      if (axes%is_q_point) axes%mask3d => diag_cs%mask3dBi

    endif

  endif

  axes%dsamp(dl)%mask2d => null()

  if (axes%rank==2) then

    if (axes%is_h_point) axes%dsamp(dl)%mask2d => diag_cs%dsamp(dl)%mask2dT

    if (axes%is_u_point) axes%dsamp(dl)%mask2d => diag_cs%dsamp(dl)%mask2dCu

    if (axes%is_v_point) axes%dsamp(dl)%mask2d => diag_cs%dsamp(dl)%mask2dCv

    if (axes%is_q_point) axes%dsamp(dl)%mask2d => diag_cs%dsamp(dl)%mask2dBu

  endif

  ! A static 3d mask for non-native coordinates can only be setup when a grid is available

  axes%dsamp(dl)%mask3d => null()

  if (axes%rank==3 .and. axes%is_native) then

    ! Native variables can/should use the native masks copied into diag_cs

    if (axes%is_layer) then

      if (axes%is_h_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dTL

      if (axes%is_u_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dCuL

      if (axes%is_v_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dCvL

      if (axes%is_q_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dBL

    elseif (axes%is_interface) then

      if (axes%is_h_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dTi

      if (axes%is_u_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dCui

      if (axes%is_v_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dCvi

      if (axes%is_q_point) axes%dsamp(dl)%mask3d => diag_cs%dsamp(dl)%mask3dBi

    endif

  endif

end subroutine define_axes_group_dsamp

!> Set up the array extents for doing diagnostics

subroutine set_diag_mediator_grid(G, diag_cs)

  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure

  type(diag_ctrl),  intent(inout) :: diag_cs !< Structure used to regulate diagnostic output

  diag_cs%is = g%isc - (g%isd-1) ; diag_cs%ie = g%iec - (g%isd-1)

  diag_cs%js = g%jsc - (g%jsd-1) ; diag_cs%je = g%jec - (g%jsd-1)

  diag_cs%isd = g%isd ; diag_cs%ied = g%ied

  diag_cs%jsd = g%jsd ; diag_cs%jed = g%jed

end subroutine set_diag_mediator_grid

!> Make a real scalar diagnostic available for averaging or output

subroutine post_data_0d(diag_field_id, field, diag_cs, is_static)

  integer,           intent(in) :: diag_field_id !< The id for an output variable returned by a

                                                 !! previous call to register_diag_field.

  real,              intent(in) :: field         !< real value being offered for output or averaging

                                                 !! in internally scaled arbitrary units [A ~> a]

  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.

  ! Local variables

  real :: locfield ! The field being offered in arbitrary unscaled units [a]

  logical :: used, is_stat

  type(diag_type), pointer :: diag => null()

  integer :: time_days

  integer :: time_seconds

  character(len=300) :: debug_mesg

  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)

  is_stat = .false. ; if (present(is_static)) is_stat = is_static

  ! Iterate over list of diag 'variants', e.g. CMOR aliases, call send_data

  ! for each one.

  call assert(diag_field_id < diag_cs%next_free_diag_id, &

              'post_data_0d: Unregistered diagnostic id')

  diag => diag_cs%diags(diag_field_id)

  do while (associated(diag))

    locfield = field

    if (diag%conversion_factor /= 0.) &

      locfield = locfield * diag%conversion_factor

    if (diag_cs%diag_as_chksum) then

      ! Append timestep to mesg

      call get_time(diag_cs%time_end, time_seconds, days=time_days)

      write(debug_mesg, '(a, 1x, i0, 1x, i0)') &

          trim(diag%debug_str), time_days, time_seconds

      call chksum0(locfield, debug_mesg, logunit=diag_cs%chksum_iounit)

    elseif (is_stat) then

      used = send_data_infra(diag%fms_diag_id, locfield)

    elseif (diag_cs%ave_enabled) then

      used = send_data_infra(diag%fms_diag_id, locfield, diag_cs%time_end)

    endif

    diag => diag%next

  enddo

  if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator)

end subroutine post_data_0d

!> Make a real 1-d array diagnostic available for averaging or output

subroutine post_data_1d_k(diag_field_id, field, diag_cs, is_static)

  integer,           intent(in) :: diag_field_id !< The id for an output variable returned by a

                                                 !! previous call to register_diag_field.

  real, target,      intent(in) :: field(:)      !< 1-d array being offered for output or averaging

                                                 !! in internally scaled arbitrary units [A ~> a]

  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.

  ! Local variables

  logical :: used  ! The return value of send_data is not used for anything.

  real, dimension(:), pointer :: locfield => null() ! The field being offered in arbitrary unscaled units [a]

  logical :: is_stat

  integer :: k, ks, ke

  type(diag_type), pointer :: diag => null()

  integer :: time_days

  integer :: time_seconds

  character(len=300) :: debug_mesg

  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)

  is_stat = .false. ; if (present(is_static)) is_stat = is_static

  ! Iterate over list of diag 'variants', e.g. CMOR aliases.

  call assert(diag_field_id < diag_cs%next_free_diag_id, &

              'post_data_1d_k: Unregistered diagnostic id')

  diag => diag_cs%diags(diag_field_id)

  do while (associated(diag))

    if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) then

      ks = lbound(field,1) ; ke = ubound(field,1)

      allocate( locfield( ks:ke ) )

      do k=ks,ke

        locfield(k) = field(k) * diag%conversion_factor

      enddo

    else

      locfield => field

    endif

    if (diag_cs%diag_as_chksum) then

      ! Append timestep to mesg

      call get_time(diag_cs%time_end, time_seconds, days=time_days)

      write(debug_mesg, '(a, 1x, i0, 1x, i0)') &

          trim(diag%debug_str), time_days, time_seconds

      call zchksum(locfield, debug_mesg, logunit=diag_cs%chksum_iounit)

    elseif (is_stat) then

      used = send_data_infra(diag%fms_diag_id, locfield)

    elseif (diag_cs%ave_enabled) then

      used = send_data_infra(diag%fms_diag_id, locfield, time=diag_cs%time_end, weight=diag_cs%time_int)

    endif

    if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield )

    diag => diag%next

  enddo

  if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator)

end subroutine post_data_1d_k

!> Make a real 2-d array diagnostic available for averaging or output

subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask)

  integer,           intent(in) :: diag_field_id !< The id for an output variable returned by a

                                                 !! previous call to register_diag_field.

  real,              intent(in) :: field(:,:)    !< 2-d array being offered for output or averaging

                                                 !! in internally scaled arbitrary units [A ~> a]

  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.

  real,    optional, intent(in) :: mask(:,:) !< If present, use this real array as the data mask [nondim]

  ! Local variables

  type(diag_type), pointer :: diag => null()

  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)

  ! Iterate over list of diag 'variants' (e.g. CMOR aliases) and post each.

  call assert(diag_field_id < diag_cs%next_free_diag_id, &

              'post_data_2d: Unregistered diagnostic id')

  diag => diag_cs%diags(diag_field_id)

  do while (associated(diag))

    call post_data_2d_low(diag, field, diag_cs, is_static, mask)

    diag => diag%next

  enddo

  if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator)

end subroutine post_data_2d

!> Make a real 2-d array diagnostic available for averaging or output

!! using a diag_type instead of an integer id.

subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask)

  type(diag_type),   intent(in) :: diag       !< A structure describing the diagnostic to post

  real,    target,   intent(in) :: field(:,:) !< 2-d array being offered for output or averaging

                                              !! in internally scaled arbitrary units [A ~> a]

  type(diag_ctrl),   intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.

  real, optional, target, intent(in) :: mask(:,:) !< If present, use this real array as the data mask [nondim]

  ! Local variables

  real, dimension(:,:), pointer :: locfield ! The field being offered in arbitrary unscaled units [a]

  real, dimension(:,:), pointer :: locmask  ! A pointer to the data mask to use [nondim]

  logical :: used  ! The return value of send_data is not used for anything.

  logical :: is_stat

  integer :: cszi, cszj, dszi, dszj

  integer :: isv, iev, jsv, jev, i, j, isv_o, jsv_o

  real, dimension(:,:), allocatable, target :: locfield_dsamp ! A downsampled version of locfield [a]

  real, dimension(:,:), allocatable, target :: locmask_dsamp  ! A downsampled version of locmask [nondim]

  integer :: dl

  integer :: time_days

  integer :: time_seconds

  character(len=300) :: mesg

  character(len=300) :: debug_mesg

  locfield => null()

  locmask => null()

  is_stat = .false. ; if (present(is_static)) is_stat = is_static

  ! Determine the proper array indices, noting that because of the (:,:)

  ! declaration of field, symmetric arrays are using a SW-grid indexing,

  ! but non-symmetric arrays are using a NE-grid indexing.  Send_data

  ! actually only uses the difference between ie and is to determine

  ! the output data size and assumes that halos are symmetric.

  isv = diag_cs%is ; iev = diag_cs%ie ; jsv = diag_cs%js ; jev = diag_cs%je

  cszi = diag_cs%ie-diag_cs%is +1 ; dszi = diag_cs%ied-diag_cs%isd +1

  cszj = diag_cs%je-diag_cs%js +1 ; dszj = diag_cs%jed-diag_cs%jsd +1

  if ( size(field,1) == dszi ) then

    isv = diag_cs%is ; iev = diag_cs%ie     ! Data domain

  elseif ( size(field,1) == dszi + 1 ) then

    isv = diag_cs%is ; iev = diag_cs%ie+1   ! Symmetric data domain

  elseif ( size(field,1) == cszi) then

    isv = 1 ; iev = cszi                    ! Computational domain

  elseif ( size(field,1) == cszi + 1 ) then

    isv = 1 ; iev = cszi+1                  ! Symmetric computational domain

  else

    write (mesg,*) " peculiar size ",size(field,1)," in i-direction\n"//&

       "does not match one of ", cszi, cszi+1, dszi, dszi+1

    call mom_error(fatal,"post_data_2d_low: "//trim(diag%debug_str)//trim(mesg))

  endif

  if ( size(field,2) == dszj ) then

    jsv = diag_cs%js ; jev = diag_cs%je     ! Data domain

  elseif ( size(field,2) == dszj + 1 ) then

    jsv = diag_cs%js ; jev = diag_cs%je+1   ! Symmetric data domain

  elseif ( size(field,2) == cszj ) then

    jsv = 1 ; jev = cszj                    ! Computational domain

  elseif ( size(field,2) == cszj+1 ) then

    jsv = 1 ; jev = cszj+1                  ! Symmetric computational domain

  else

    write (mesg,*) " peculiar size ",size(field,2)," in j-direction\n"//&

       "does not match one of ", cszj, cszj+1, dszj, dszj+1

    call mom_error(fatal,"post_data_2d_low: "//trim(diag%debug_str)//trim(mesg))

  endif

  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) then

    allocate( locfield( lbound(field,1):ubound(field,1), lbound(field,2):ubound(field,2) ) )

    do j=jsv,jev ; do i=isv,iev

      locfield(i,j) = field(i,j) * diag%conversion_factor

    enddo ; enddo

  else

    locfield => field

  endif

  if (present(mask)) then

    locmask => mask

  elseif (.NOT. is_stat .and. associated(diag%axes)) then

    if (associated(diag%axes%mask2d)) locmask => diag%axes%mask2d

  endif

  dl = 1

  if (.NOT. is_stat .and. associated(diag%axes)) &

    dl = diag%axes%downsample_level ! Static field downsampling is not supported yet.

  ! Downsample the diag field and mask as appropriate.

  if (dl > 1) then

    isv_o = isv ; jsv_o = jsv

    call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag, isv, iev, jsv, jev, mask)

    if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield )

    locfield => locfield_dsamp

    if (present(mask)) then

      call downsample_field_2d(locmask, locmask_dsamp, dl, msk, locmask, diag_cs, diag, &

                               isv_o, jsv_o, isv, iev, jsv, jev)

      locmask => locmask_dsamp

    elseif (associated(diag%axes%dsamp(dl)%mask2d)) then

      locmask => diag%axes%dsamp(dl)%mask2d

    endif

  endif

  if (diag_cs%diag_as_chksum) then

    ! Append timestep to mesg

    call get_time(diag_cs%time_end, time_seconds, days=time_days)

    write(debug_mesg, '(a, 1x, i0, 1x, i0)') &

        trim(diag%debug_str), time_days, time_seconds

    if (diag%axes%is_h_point) then

      call hchksum(locfield, debug_mesg, diag_cs%G%HI, &

                   logunit=diag_cs%chksum_iounit)

    elseif (diag%axes%is_u_point) then

      call uchksum(locfield, debug_mesg, diag_cs%G%HI, &

                   logunit=diag_cs%chksum_iounit)

    elseif (diag%axes%is_v_point) then

      call vchksum(locfield, debug_mesg, diag_cs%G%HI, &

                   logunit=diag_cs%chksum_iounit)

    elseif (diag%axes%is_q_point) then

      call bchksum(locfield, debug_mesg, diag_cs%G%HI, &

                   logunit=diag_cs%chksum_iounit)

    else

      call mom_error(fatal, "post_data_2d_low: unknown axis type.")

    endif

  else

    if (is_stat) then

      if (present(mask)) then

        call assert(size(locfield) == size(locmask), &

            'post_data_2d_low is_stat: mask size mismatch: '//trim(diag%debug_str))

        used = send_data_infra(diag%fms_diag_id, locfield, &

                         is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, rmask=locmask)

     !elseif (associated(diag%axes%mask2d)) then

     !  used = send_data(diag%fms_diag_id, locfield, &

     !                   is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, rmask=diag%axes%mask2d)

      else

        used = send_data_infra(diag%fms_diag_id, locfield, &

                         is_in=isv, ie_in=iev, js_in=jsv, je_in=jev)

      endif

    elseif (diag_cs%ave_enabled) then

      if (associated(locmask)) then

        call assert(size(locfield) == size(locmask), &

            'post_data_2d_low: mask size mismatch: '//trim(diag%debug_str))

        used = send_data_infra(diag%fms_diag_id, locfield, &

                         is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, &

                         time=diag_cs%time_end, weight=diag_cs%time_int, rmask=locmask)

      else

        used = send_data_infra(diag%fms_diag_id, locfield, &

                         is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, &

                         time=diag_cs%time_end, weight=diag_cs%time_int)

      endif

    endif

  endif

  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.) .and. dl<2) &

    deallocate( locfield )

end subroutine post_data_2d_low

!> Make a real 3-d array diagnostic available for averaging or output.

subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h)

  integer,           intent(in) :: diag_field_id !< The id for an output variable returned by a

                                                 !! previous call to register_diag_field.

  real,              intent(in) :: field(:,:,:)  !< 3-d array being offered for output or averaging

                                                 !! in internally scaled arbitrary units [A ~> a]

  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.

  real,    optional, intent(in) :: mask(:,:,:) !< If present, use this real array as the data mask [nondim]

  real, dimension(:,:,:), &

         target, optional, intent(in) :: alt_h  !< An alternate thickness to use for vertically

                                                !! remapping this diagnostic [H ~> m or kg m-2].

  ! Local variables

  type(diag_type), pointer :: diag => null()

  real, dimension(:,:,:), allocatable :: remapped_field !< The vertically remapped diagnostic [A ~> a]

  logical :: staggered_in_x, staggered_in_y, dz_diag_needed, dz_begin_needed

  real, dimension(:,:,:), pointer :: h_diag => null() !< A pointer to the thickness to use for vertically

                                                !! remapping this diagnostic [H ~> m or kg m-2].

  real, dimension(diag_cs%G%isd:diag_cS%G%ied, diag_cs%G%jsd:diag_cS%G%jed, diag_cs%GV%ke) :: &

    dz_diag     ! Layer vertical extents for remapping [Z ~> m]

  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)

  ! For intensive variables only, we can choose to use a different diagnostic grid to map to

  if (present(alt_h)) then

    h_diag => alt_h

  else

    h_diag => diag_cs%h

  endif

  ! Iterate over list of diag 'variants', e.g. CMOR aliases, different vertical

  ! grids, and post each.

  call assert(diag_field_id < diag_cs%next_free_diag_id, &

              'post_data_3d: Unregistered diagnostic id')

  if (diag_cs%show_call_tree) &

    call calltree_enter("post_data_3d("//trim(diag_cs%diags(diag_field_id)%debug_str)//")")

  ! Find out whether there are any z-based diagnostics

  diag => diag_cs%diags(diag_field_id)

  dz_diag_needed = .false.

  do while (associated(diag))

    if (diag%axes%needs_remapping .or. diag%axes%needs_interpolating) then

      if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) &

        dz_diag_needed = .true.

    endif

    diag => diag%next

  enddo

  ! Determine the diagnostic grid spacing in height units, if it is needed.

  if (dz_diag_needed) then

    call thickness_to_dz(h_diag, diag_cs%tv, dz_diag, diag_cs%G, diag_cs%GV, diag_cs%US, halo_size=1)

  endif

  diag => diag_cs%diags(diag_field_id)

  do while (associated(diag))

    call assert(associated(diag%axes), 'post_data_3d: axes is not associated')

    staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point

    staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point

    if (diag%v_extensive .and. .not.diag%axes%is_native) then

      ! The field is vertically integrated and needs to be re-gridded

      if (present(mask)) then

        call mom_error(fatal,"post_data_3d: no mask for regridded field.")

      endif

      if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)

      allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz))

      if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) then

        call vertically_reintegrate_diag_field(                                    &

                diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, &

                diag_cs%dz_begin, diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, &

                staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field)

      else

        call vertically_reintegrate_diag_field(                                    &

                diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, &

                diag_cs%h_begin, diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, &

                staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field)

      endif

      if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)

      if (associated(diag%axes%mask3d)) then

        ! Since 3d masks do not vary in the vertical, just use as much as is

        ! needed.

        call post_data_3d_low(diag, remapped_field, diag_cs, is_static, &

                              mask=diag%axes%mask3d)

      else

        call post_data_3d_low(diag, remapped_field, diag_cs, is_static)

      endif

      if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)

      deallocate(remapped_field)

      if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)

    elseif (diag%axes%needs_remapping) then

      ! Remap this field to another vertical coordinate.

      if (present(mask)) then

        call mom_error(fatal,"post_data_3d: no mask for regridded field.")

      endif

      if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)

      allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz))

      if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) then

        call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), &

                diag_cs%G, diag_cs%GV, diag_cs%US, dz_diag, staggered_in_x, staggered_in_y, &

                diag%axes%mask3d, field, remapped_field)

      else

        call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), &

                diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, staggered_in_x, staggered_in_y, &

                diag%axes%mask3d, field, remapped_field)

      endif

      if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)

      if (associated(diag%axes%mask3d)) then

        ! Since 3d masks do not vary in the vertical, just use as much as is

        ! needed.

        call post_data_3d_low(diag, remapped_field, diag_cs, is_static, &

                              mask=diag%axes%mask3d)

      else

        call post_data_3d_low(diag, remapped_field, diag_cs, is_static)

      endif

      if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)

      deallocate(remapped_field)

      if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)

    elseif (diag%axes%needs_interpolating) then

      ! Interpolate this field to another vertical coordinate.

      if (present(mask)) then

        call mom_error(fatal,"post_data_3d: no mask for regridded field.")

      endif

      if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)

      allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz+1))

      if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) then

        call vertically_interpolate_diag_field(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), &

                diag_cs%G, dz_diag, staggered_in_x, staggered_in_y, &

                diag%axes%mask3d, field, remapped_field)

      else

        call vertically_interpolate_diag_field(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), &

                diag_cs%G, h_diag, staggered_in_x, staggered_in_y, &

                diag%axes%mask3d, field, remapped_field)

      endif

      if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)

      if (associated(diag%axes%mask3d)) then

        ! Since 3d masks do not vary in the vertical, just use as much as is needed.

        call post_data_3d_low(diag, remapped_field, diag_cs, is_static, &

                              mask=diag%axes%mask3d)

      else

        call post_data_3d_low(diag, remapped_field, diag_cs, is_static)

      endif

      if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)

      deallocate(remapped_field)

      if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)

    else

      call post_data_3d_low(diag, field, diag_cs, is_static, mask)

    endif

    diag => diag%next

  enddo

  if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator)

  if (diag_cs%show_call_tree) &

    call calltree_leave("post_data_3d("//trim(diag_cs%diags(diag_field_id)%debug_str)//")")

end subroutine post_data_3d

!> Make a real 3-d array diagnostic available for averaging or output

!! using a diag_type instead of an integer id.

subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask)

  type(diag_type),   intent(in) :: diag       !< A structure describing the diagnostic to post

  real,    target,   intent(in) :: field(:,:,:) !< 3-d array being offered for output or averaging

                                                !! in internally scaled arbitrary units [A ~> a]

  type(diag_ctrl),   intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered.

  real,    optional,target, intent(in) :: mask(:,:,:) !< If present, use this real array as the data mask [nondim]

  ! Local variables

  real, dimension(:,:,:), pointer :: locfield ! The field being offered in arbitrary unscaled units [a]

  real, dimension(:,:,:), pointer :: locmask  ! A pointer to the data mask to use [nondim]

  character(len=300) :: mesg

  logical :: used  ! The return value of send_data is not used for anything.

  logical :: staggered_in_x, staggered_in_y

  logical :: is_stat

  integer :: cszi, cszj, dszi, dszj

  integer :: isv, iev, jsv, jev, ks, ke, i, j, k, isv_c, jsv_c, isv_o, jsv_o

  real, dimension(:,:,:), allocatable, target :: locfield_dsamp ! A downsampled version of locfield [a]

  real, dimension(:,:,:), allocatable, target :: locmask_dsamp  ! A downsampled version of locmask [nondim]

  integer :: dl

  integer :: time_days

  integer :: time_seconds

  character(len=300) :: debug_mesg

  locfield => null()

  locmask => null()

  is_stat = .false. ; if (present(is_static)) is_stat = is_static

  ! Determine the proper array indices, noting that because of the (:,:)

  ! declaration of field, symmetric arrays are using a SW-grid indexing,

  ! but non-symmetric arrays are using a NE-grid indexing.  Send_data

  ! actually only uses the difference between ie and is to determine

  ! the output data size and assumes that halos are symmetric.

  isv = diag_cs%is ; iev = diag_cs%ie ; jsv = diag_cs%js ; jev = diag_cs%je

  cszi = (diag_cs%ie-diag_cs%is) +1 ; dszi = (diag_cs%ied-diag_cs%isd) +1

  cszj = (diag_cs%je-diag_cs%js) +1 ; dszj = (diag_cs%jed-diag_cs%jsd) +1

  if ( size(field,1) == dszi ) then

    isv = diag_cs%is ; iev = diag_cs%ie     ! Data domain

  elseif ( size(field,1) == dszi + 1 ) then

    isv = diag_cs%is ; iev = diag_cs%ie+1   ! Symmetric data domain

  elseif ( size(field,1) == cszi) then

    isv = 1 ; iev = cszi                    ! Computational domain

  elseif ( size(field,1) == cszi + 1 ) then

    isv = 1 ; iev = cszi+1                  ! Symmetric computational domain

  else

    write (mesg,*) " peculiar size ",size(field,1)," in i-direction\n"//&

       "does not match one of ", cszi, cszi+1, dszi, dszi+1

    call mom_error(fatal,"post_data_3d_low: "//trim(diag%debug_str)//trim(mesg))

  endif

  if ( size(field,2) == dszj ) then

    jsv = diag_cs%js ; jev = diag_cs%je     ! Data domain

  elseif ( size(field,2) == dszj + 1 ) then

    jsv = diag_cs%js ; jev = diag_cs%je+1   ! Symmetric data domain

  elseif ( size(field,2) == cszj ) then

    jsv = 1 ; jev = cszj                    ! Computational domain

  elseif ( size(field,2) == cszj+1 ) then

    jsv = 1 ; jev = cszj+1                  ! Symmetric computational domain

  else

    write (mesg,*) " peculiar size ",size(field,2)," in j-direction\n"//&

       "does not match one of ", cszj, cszj+1, dszj, dszj+1

    call mom_error(fatal,"post_data_3d_low: "//trim(diag%debug_str)//trim(mesg))

  endif

  ks = lbound(field,3) ; ke = ubound(field,3)

  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) then

    allocate( locfield( lbound(field,1):ubound(field,1), lbound(field,2):ubound(field,2), ks:ke ) )

    ! locfield(:,:,:) = 0.0  ! Zeroing out this array would be a good idea, but it appears

                             ! not to be necessary.

    isv_c = isv ; jsv_c = jsv

    if (diag%fms_xyave_diag_id>0) then

      staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point

      staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point

      ! When averaging a staggered field, edge points are always required.

      if (staggered_in_x) isv_c = iev - (diag_cs%ie - diag_cs%is) - 1

      if (staggered_in_y) jsv_c = jev - (diag_cs%je - diag_cs%js) - 1

      if (isv_c < lbound(locfield,1)) call mom_error(fatal, &

        "It is an error to average a staggered diagnostic field that does not "//&

        "have i-direction space to represent the symmetric computational domain.")

      if (jsv_c < lbound(locfield,2)) call mom_error(fatal, &

        "It is an error to average a staggered diagnostic field that does not "//&

        "have j-direction space to represent the symmetric computational domain.")

    endif

    do k=ks,ke ; do j=jsv,jev ; do i=isv,iev

      locfield(i,j,k) = field(i,j,k) * diag%conversion_factor

    enddo ; enddo ; enddo

  else

    locfield => field

  endif

  if (present(mask)) then

    locmask => mask

  elseif (associated(diag%axes%mask3d)) then

    locmask => diag%axes%mask3d

  endif

  dl = 1

  if (.NOT. is_stat) &

    dl = diag%axes%downsample_level ! Static field downsampling is not supported yet.

  ! Downsample the diag field and mask as appropriate.

  if (dl > 1) then

    isv_o = isv ; jsv_o = jsv

    call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag, isv, iev, jsv, jev, mask)

    if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield )

    locfield => locfield_dsamp

    if (present(mask)) then

      call downsample_field_3d(locmask, locmask_dsamp, dl, msk, locmask, diag_cs, diag, &

                               isv_o, jsv_o, isv, iev, jsv, jev)

      locmask => locmask_dsamp

    elseif (associated(diag%axes%dsamp(dl)%mask3d)) then

      locmask => diag%axes%dsamp(dl)%mask3d

    endif

  endif

  if (diag%fms_diag_id>0) then

    if (diag_cs%diag_as_chksum) then

      ! Append timestep to mesg

      call get_time(diag_cs%time_end, time_seconds, days=time_days)

      write(debug_mesg, '(a, 1x, i0, 1x, i0)') &

          trim(diag%debug_str), time_days, time_seconds

      if (diag%axes%is_h_point) then

        call hchksum(locfield, debug_mesg, diag_cs%G%HI, &

                     logunit=diag_cs%chksum_iounit)

      elseif (diag%axes%is_u_point) then

        call uchksum(locfield, debug_mesg, diag_cs%G%HI, &

                     logunit=diag_cs%chksum_iounit)

      elseif (diag%axes%is_v_point) then

        call vchksum(locfield, debug_mesg, diag_cs%G%HI, &

                     logunit=diag_cs%chksum_iounit)

      elseif (diag%axes%is_q_point) then

        call bchksum(locfield, debug_mesg, diag_cs%G%HI, &

                     logunit=diag_cs%chksum_iounit)

      else

        call mom_error(fatal, "post_data_3d_low: unknown axis type.")

      endif

    else

      if (is_stat) then

        if (present(mask)) then

          call assert(size(locfield) == size(locmask), &

              'post_data_3d_low is_stat: mask size mismatch: '//diag%debug_str)

          used = send_data_infra(diag%fms_diag_id, locfield, &

                         is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, rmask=locmask)

       !elseif (associated(diag%axes%mask2d)) then

       !  used = send_data(diag%fms_diag_id, locfield, &

       !                   is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, rmask=diag%axes%mask2d)

        else

          used = send_data_infra(diag%fms_diag_id, locfield, &

                           is_in=isv, ie_in=iev, js_in=jsv, je_in=jev)

        endif

      elseif (diag_cs%ave_enabled) then

        if (associated(locmask)) then

          call assert(size(locfield) == size(locmask), &

              'post_data_3d_low: mask size mismatch: '//diag%debug_str)

          used = send_data_infra(diag%fms_diag_id, locfield, &

                           is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, &

                           time=diag_cs%time_end, weight=diag_cs%time_int, rmask=locmask)

        else

          used = send_data_infra(diag%fms_diag_id, locfield, &

                           is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, &

                           time=diag_cs%time_end, weight=diag_cs%time_int)

        endif

      endif

    endif

  endif

  if (diag%fms_xyave_diag_id>0) then

    call post_xy_average(diag_cs, diag, locfield)

  endif

  if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.) .and. dl<2) &

    deallocate( locfield )

end subroutine post_data_3d_low

!> Put data into the buffer for a diagnostic one column at a time

subroutine post_data_3d_by_column(diag_field_id, field, diag_cs, i, j)

  integer,            intent(in) :: diag_field_id !< The id for an output variable returned by a

                                                  !! previous call to register_diag_field.

  real, dimension(:), intent(in) :: field         !< 3-d array being offered for output or averaging

                                                  !! in internally scaled arbitrary units [A ~> a]

  type(diag_ctrl), target, intent(in) :: diag_cs  !< Structure used to regulate diagnostic output

  integer,            intent(in) :: i             !< The i-index to post the data in the buffer

  integer,            intent(in) :: j             !< The j-index to post the data in the buffer

  type(diag_type), pointer :: diag => null()

  integer :: buffer_slot

  diag => diag_cs%diags(diag_field_id)

  buffer_slot = diag%axes%piecemeal_3d%check_capacity_by_id(diag_field_id)

  diag%axes%piecemeal_3d%buffer(buffer_slot)%field(i,j,:) = field(:)

end subroutine post_data_3d_by_column

!> Put data into the buffer for a diagnostic one point at a time

subroutine post_data_3d_by_point(diag_field_id, field, diag_cs, i, j, k)

  integer,           intent(in) :: diag_field_id !< The id for an output variable returned by a

                                                 !! previous call to register_diag_field.

  real,              intent(in) :: field         !< 3-d array being offered for output or averaging

                                                 !! in internally scaled arbitrary units [A ~> a]

  type(diag_ctrl), target, intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  integer,           intent(in) :: i             !< The i-index to post the data in the buffer

  integer,           intent(in) :: j             !< The j-index to post the data in the buffer

  integer,           intent(in) :: k             !< The k-index to post the data in the buffer

  type(diag_type), pointer :: diag => null()

  integer :: buffer_slot

  diag => diag_cs%diags(diag_field_id)

  buffer_slot = diag%axes%piecemeal_3d%check_capacity_by_id(diag_field_id)

  diag%axes%piecemeal_3d%buffer(buffer_slot)%field(i,j,k) = field

end subroutine post_data_3d_by_point

!> Post the final buffer using the standard post_data interface

subroutine post_data_3d_final(diag_field_id, diag_cs)

  integer,           intent(in) :: diag_field_id !< The id for an output variable returned by a

                                                 !! previous call to register_diag_field.

  type(diag_ctrl), target, intent(in) :: diag_cs !< Structure used to regulate diagnostic output

  type(diag_type), pointer :: diag => null()

  integer :: buffer_slot

  diag => diag_cs%diags(diag_field_id)

  buffer_slot = diag%axes%piecemeal_3d%find_buffer_slot(diag_field_id)

  ! Only perform an action if the buffer slot was actually used

  if (buffer_slot>0) then

    call post_data(diag_field_id, diag%axes%piecemeal_3d%buffer(buffer_slot)%field(:,:,:), diag_cs)

    call diag%axes%piecemeal_3d%mark_available(diag_field_id)

  endif

end subroutine post_data_3d_final

!> Calculate and write out diagnostics that are the product of two 3-d arrays at u-points

subroutine post_product_u(id, u_a, u_b, G, nz, diag, mask, alt_h)

  integer,                  intent(in) :: id   !< The ID for this diagnostic

  type(ocean_grid_type),    intent(in) :: g    !< ocean grid structure

  integer,                  intent(in) :: nz   !< The size of the arrays in the vertical

  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz), &

                            intent(in) :: u_a  !< The first u-point array in arbitrary units [A]

  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz), &

                            intent(in) :: u_b  !< The second u-point array in arbitrary units [B]

  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output

  real,           optional, intent(in) :: mask(:,:,:)  !< If present, use this real array as the data mask [nondim]

  real,   target, optional, intent(in) :: alt_h(:,:,:) !< An alternate thickness to use for vertically

                                               !! remapping this diagnostic [H ~> m or kg m-2]

  ! Local variables

  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz) :: u_prod ! The product of u_a and u_b [A B]

  integer :: i, j, k

  if (id <= 0) return

  do k=1,nz ; do j=g%jsc,g%jec ; do i=g%IscB,g%IecB

    u_prod(i,j,k) = u_a(i,j,k) * u_b(i,j,k)

  enddo ; enddo ; enddo

  call post_data(id, u_prod, diag, mask=mask, alt_h=alt_h)

end subroutine post_product_u

!> Calculate and write out diagnostics that are the vertical sum of the product of two 3-d arrays at u-points

subroutine post_product_sum_u(id, u_a, u_b, G, nz, diag)

  integer,                  intent(in) :: id   !< The ID for this diagnostic

  type(ocean_grid_type),    intent(in) :: g    !< ocean grid structure

  integer,                  intent(in) :: nz   !< The size of the arrays in the vertical

  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz), &

                            intent(in) :: u_a  !< The first u-point array in arbitrary units [A]

  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed, nz), &

                            intent(in) :: u_b  !< The second u-point array in arbitrary units [B]

  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output

  real, dimension(G%IsdB:G%IedB, G%jsd:G%jed) :: u_sum  ! The vertical sum of the product of u_a and u_b [A B]

  integer :: i, j, k

  if (id <= 0) return

  u_sum(:,:) = 0.0

  do k=1,nz ; do j=g%jsc,g%jec ; do i=g%IscB,g%IecB

    u_sum(i,j) = u_sum(i,j) + u_a(i,j,k) * u_b(i,j,k)

  enddo ; enddo ; enddo

  call post_data(id, u_sum, diag)

end subroutine post_product_sum_u

!> Calculate and write out diagnostics that are the product of two 3-d arrays at v-points

subroutine post_product_v(id, v_a, v_b, G, nz, diag, mask, alt_h)

  integer,                  intent(in) :: id   !< The ID for this diagnostic

  type(ocean_grid_type),    intent(in) :: g    !< ocean grid structure

  integer,                  intent(in) :: nz   !< The size of the arrays in the vertical

  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz), &

                            intent(in) :: v_a  !< The first v-point array in arbitrary units [A]

  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz), &

                            intent(in) :: v_b  !< The second v-point array in arbitrary units [B]

  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output

  real,           optional, intent(in) :: mask(:,:,:)  !< If present, use this real array as the data mask [nondim]

  real,   target, optional, intent(in) :: alt_h(:,:,:) !< An alternate thickness to use for vertically

                                               !! remapping this diagnostic [H ~> m or kg m-2]

  ! Local variables

  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz) :: v_prod ! The product of v_a and v_b [A B]

  integer :: i, j, k

  if (id <= 0) return

  do k=1,nz ; do j=g%JscB,g%JecB ; do i=g%isc,g%iec

    v_prod(i,j,k) = v_a(i,j,k) * v_b(i,j,k)

  enddo ; enddo ; enddo

  call post_data(id, v_prod, diag, mask=mask, alt_h=alt_h)

end subroutine post_product_v

!> Calculate and write out diagnostics that are the vertical sum of the product of two 3-d arrays at v-points

subroutine post_product_sum_v(id, v_a, v_b, G, nz, diag)

  integer,                  intent(in) :: id   !< The ID for this diagnostic

  type(ocean_grid_type),    intent(in) :: g    !< ocean grid structure

  integer,                  intent(in) :: nz   !< The size of the arrays in the vertical

  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz), &

                            intent(in) :: v_a  !< The first v-point array in arbitrary units [A]

  real, dimension(G%isd:G%ied, G%JsdB:G%JedB, nz), &

                            intent(in) :: v_b  !< The second v-point array in arbitrary units [B]

  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output

  real, dimension(G%isd:G%ied, G%JsdB:G%JedB) :: v_sum ! The vertical sum of the product of v_a and v_b [A B]

  integer :: i, j, k

  if (id <= 0) return

  v_sum(:,:) = 0.0

  do k=1,nz ; do j=g%JscB,g%JecB ; do i=g%isc,g%iec

    v_sum(i,j) = v_sum(i,j) + v_a(i,j,k) * v_b(i,j,k)

  enddo ; enddo ; enddo

  call post_data(id, v_sum, diag)

end subroutine post_product_sum_v

!> Post the horizontally area-averaged diagnostic

subroutine post_xy_average(diag_cs, diag, field)

  type(diag_type),   intent(in) :: diag !< This diagnostic

  real,    target,   intent(in) :: field(:,:,:) !< Diagnostic field in arbitrary units [A ~> a]

  type(diag_ctrl),   intent(in) :: diag_cs !< Diagnostics mediator control structure

  ! Local variable

  real, dimension(size(field,3)) :: averaged_field ! The horizontally averaged field [A ~> a]

  logical, dimension(size(field,3)) :: averaged_mask

  logical :: staggered_in_x, staggered_in_y, used

  integer :: nz, remap_nz, coord

  integer :: time_days

  integer :: time_seconds

  character(len=300) :: debug_mesg

  if (.not. diag_cs%ave_enabled) then

    return

  endif

  staggered_in_x = diag%axes%is_u_point .or. diag%axes%is_q_point

  staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point

  if (diag%axes%is_native) then

    call horizontally_average_diag_field(diag_cs%G, diag_cs%GV, diag_cs%h, &

                                         staggered_in_x, staggered_in_y, &

                                         diag%axes%is_layer, diag%v_extensive, &

                                         field, averaged_field, averaged_mask)

  else

    nz = size(field, 3)

    coord = diag%axes%vertical_coordinate_number

    remap_nz = diag_cs%diag_remap_cs(coord)%nz

    call assert(diag_cs%diag_remap_cs(coord)%initialized, &

                'post_xy_average: remap_cs not initialized.')

    call assert(implies(diag%axes%is_layer, nz == remap_nz), &

              'post_xy_average: layer field dimension mismatch.')

    call assert(implies(.not. diag%axes%is_layer, nz == remap_nz+1), &

              'post_xy_average: interface field dimension mismatch.')

    call horizontally_average_diag_field(diag_cs%G, diag_cs%GV, &

                                         diag_cs%diag_remap_cs(coord)%h, &

                                         staggered_in_x, staggered_in_y, &

                                         diag%axes%is_layer, diag%v_extensive, &

                                         field, averaged_field, averaged_mask)

  endif

  if (diag_cs%diag_as_chksum) then

    ! Append timestep to mesg

    call get_time(diag_cs%time_end, time_seconds, days=time_days)

    write(debug_mesg, '(a, 1x, i0, 1x, i0)') &

        trim(diag%debug_str)//'_xyave', time_days, time_seconds

    call zchksum(averaged_field, debug_mesg, logunit=diag_cs%chksum_iounit)

  else

    used = send_data_infra(diag%fms_xyave_diag_id, averaged_field, &

                           time=diag_cs%time_end, weight=diag_cs%time_int, mask=averaged_mask)

  endif

end subroutine post_xy_average

!> This subroutine enables the accumulation of time averages over the specified time interval.

subroutine enable_averaging(time_int_in, time_end_in, diag_cs)

  real,            intent(in)    :: time_int_in !< The time interval [s] over which any

                                                !! values that are offered are valid.

  type(time_type), intent(in)    :: time_end_in !< The end time of the valid interval

  type(diag_ctrl), intent(inout) :: diag_cs     !< Structure used to regulate diagnostic output

! This subroutine enables the accumulation of time averages over the specified time interval.

!  if (num_file==0) return

  diag_cs%time_int = time_int_in

  diag_cs%time_end = time_end_in

  diag_cs%ave_enabled = .true.

end subroutine enable_averaging

!> Enable the accumulation of time averages over the specified time interval in time units.

subroutine enable_averages(time_int, time_end, diag_CS, T_to_s)

  real,            intent(in)    :: time_int !< The time interval over which any values

                                             !! that are offered are valid [T ~> s].

  type(time_type), intent(in)    :: time_end !< The end time of the valid interval.

  type(diag_ctrl), intent(inout) :: diag_cs  !< A structure that is used to regulate diagnostic output

  real,  optional, intent(in)    :: t_to_s   !< A conversion factor for time_int to seconds [s T-1 ~> 1].

  ! This subroutine enables the accumulation of time averages over the specified time interval.

  if (present(t_to_s)) then

    diag_cs%time_int = time_int*t_to_s

  elseif (associated(diag_cs%US)) then

    diag_cs%time_int = time_int*diag_cs%US%T_to_s

  else

    diag_cs%time_int = time_int

  endif

  diag_cs%time_end = time_end

  diag_cs%ave_enabled = .true.

end subroutine enable_averages

!> Call this subroutine to avoid averaging any offered fields.

subroutine disable_averaging(diag_cs)

  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output

  diag_cs%time_int = 0.0

  diag_cs%ave_enabled = .false.

end subroutine disable_averaging

!> Call this subroutine to determine whether the averaging is

!! currently enabled.  .true. is returned if it is.

function query_averaging_enabled(diag_cs, time_int, time_end)

  type(diag_ctrl),           intent(in)  :: diag_cs  !< Structure used to regulate diagnostic output

  real,            optional, intent(out) :: time_int !< Current setting of diag%time_int [s]

  type(time_type), optional, intent(out) :: time_end !< Current setting of diag%time_end

  logical :: query_averaging_enabled

  if (present(time_int)) time_int = diag_cs%time_int

  if (present(time_end)) time_end = diag_cs%time_end

  query_averaging_enabled = diag_cs%ave_enabled

end function query_averaging_enabled

!> This function returns the valid end time for use with diagnostics that are

!! handled outside of the MOM6 diagnostics infrastructure.

function get_diag_time_end(diag_cs)

  type(diag_ctrl), intent(in)  :: diag_cs !< Structure used to regulate diagnostic output

  type(time_type) :: get_diag_time_end

  !   This function returns the valid end time for diagnostics that are handled

  ! outside of the MOM6 infrastructure, such as via the generic tracer code.

  get_diag_time_end = diag_cs%time_end

end function get_diag_time_end

!> Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics

!! derived from one field.

integer function 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)

  character(len=*),           intent(in) :: module_name !< Name of this module, usually "ocean_model"

                                                        !! or "ice_shelf_model"

  character(len=*),           intent(in) :: field_name !< Name of the diagnostic field

  type(axes_grp),     target, intent(in) :: axes_in   !< Container w/ up to 3 integer handles that

                                                      !! indicates axes for this field

  type(time_type),            intent(in) :: init_time !< Time at which a field is first available?

  character(len=*), optional, intent(in) :: long_name !< Long name of a field.

  character(len=*), optional, intent(in) :: units !< Units of a field.

  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field

  real,             optional, intent(in) :: missing_value !< A value that indicates missing values in

                                                          !! output files, in unscaled arbitrary units [a]

  real,             optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)

                                                     !! in arbitrary units [a]

  logical,          optional, intent(in) :: mask_variant !< If true a logical mask must be provided with

                                                         !! post_data calls (not used in MOM?)

  logical,          optional, intent(in) :: verbose !< If true, FMS is verbose (not used in MOM?)

  logical,          optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?)

  character(len=*), optional, intent(out):: err_msg !< String into which an error message might be

                                                         !! placed (not used in MOM?)

  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not

                                                         !! be interpolated as a scalar

  integer,          optional, intent(in) :: tile_count !< no clue (not used in MOM?)

  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field

  character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field

  character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field

  character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field

  character(len=*), optional, intent(in) :: cell_methods !< String to append as cell_methods attribute. Use '' to

                                                         !! have no attribute.  If present, this overrides the

                                                         !! default constructed from the default for

                                                         !! each individual axis direction.

  character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction.

                                                         !! Use '' have no method.

  character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction.

                                                         !! Use '' have no method.

  character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction.

                                                         !! Use '' have no method.

  real,             optional, intent(in) :: conversion !< A value to multiply data by before writing to files,

                                                       !! often including factors to undo internal scaling and

                                                       !! in units of [a A-1 ~> 1]

  logical,          optional, intent(in) :: v_extensive  !< True for vertically extensive fields (vertically

                                                         !! integrated). Default/absent for intensive.

  ! Local variables

  real :: mom_missing_value ! A value used to indicate missing values in output files, in arbitrary units [a]

  type(diag_ctrl), pointer :: diag_cs => null() ! A structure that is used to regulate diagnostic output

  type(axes_grp), pointer :: remap_axes

  type(axes_grp), pointer :: axes

  type(axes_grp), pointer :: axes_d2

  integer :: dm_id, i, dl

  character(len=256) :: msg, cm_string

  character(len=256) :: new_module_name

  character(len=480) :: module_list, var_list

  character(len=16)  :: dimensions

  integer :: num_modnm, num_varnm

  logical :: active

  diag_cs => axes_in%diag_cs

  ! Check if the axes match a standard grid axis.

  ! If not, allocate the new axis and copy the contents.

  if (axes_in%id == diag_cs%axesTL%id) then

    axes => diag_cs%axesTL

  elseif (axes_in%id == diag_cs%axesBL%id) then

    axes => diag_cs%axesBL

  elseif (axes_in%id == diag_cs%axesCuL%id) then

    axes => diag_cs%axesCuL

  elseif (axes_in%id == diag_cs%axesCvL%id) then

    axes => diag_cs%axesCvL

  elseif (axes_in%id == diag_cs%axesTi%id) then

    axes => diag_cs%axesTi

  elseif (axes_in%id == diag_cs%axesBi%id) then

    axes => diag_cs%axesBi

  elseif (axes_in%id == diag_cs%axesCui%id) then

    axes => diag_cs%axesCui

  elseif (axes_in%id == diag_cs%axesCvi%id) then

    axes => diag_cs%axesCvi

  else

    allocate(axes)

    axes = axes_in

  endif

  mom_missing_value = axes%diag_cs%missing_value

  if (present(missing_value)) mom_missing_value = missing_value

  diag_cs => axes%diag_cs

  dm_id = -1

  module_list = "{"//trim(module_name)

  num_modnm = 1

  ! Register the native diagnostic

  active = register_diag_field_expand_cmor(dm_id, module_name, field_name, axes, &

             init_time, long_name=long_name, units=units, missing_value=mom_missing_value, &

             range=range, mask_variant=mask_variant, standard_name=standard_name, &

             verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

             interp_method=interp_method, tile_count=tile_count, &

             cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, &

             cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, &

             cell_methods=cell_methods, x_cell_method=x_cell_method, &

             y_cell_method=y_cell_method, v_cell_method=v_cell_method, &

             conversion=conversion, v_extensive=v_extensive)

  if (associated(axes%xyave_axes)) then

    num_varnm = 2 ; var_list = "{"//trim(field_name)//","//trim(field_name)//"_xyave"

  else

    num_varnm = 1 ; var_list = "{"//trim(field_name)

  endif

  if (present(cmor_field_name)) then

    if (associated(axes%xyave_axes)) then

      num_varnm = num_varnm + 2

      var_list = trim(var_list)//","//trim(cmor_field_name)//","//trim(cmor_field_name)//"_xyave"

    else

      num_varnm = num_varnm + 1

      var_list = trim(var_list)//","//trim(cmor_field_name)

    endif

  endif

  var_list = trim(var_list)//"}"

  ! For each diagnostic coordinate register the diagnostic again under a different module name

  do i=1,diag_cs%num_diag_coords

    new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)

    ! Register diagnostics remapped to z vertical coordinate

    if (axes_in%rank == 3) then

      remap_axes => null()

      if ((axes_in%id == diag_cs%axesTL%id)) then

        remap_axes => diag_cs%remap_axesTL(i)

      elseif (axes_in%id == diag_cs%axesBL%id) then

        remap_axes => diag_cs%remap_axesBL(i)

      elseif (axes_in%id == diag_cs%axesCuL%id ) then

        remap_axes => diag_cs%remap_axesCuL(i)

      elseif (axes_in%id == diag_cs%axesCvL%id) then

        remap_axes => diag_cs%remap_axesCvL(i)

      elseif (axes_in%id == diag_cs%axesTi%id) then

        remap_axes => diag_cs%remap_axesTi(i)

      elseif (axes_in%id == diag_cs%axesBi%id) then

        remap_axes => diag_cs%remap_axesBi(i)

      elseif (axes_in%id == diag_cs%axesCui%id ) then

        remap_axes => diag_cs%remap_axesCui(i)

      elseif (axes_in%id == diag_cs%axesCvi%id) then

        remap_axes => diag_cs%remap_axesCvi(i)

      endif

      ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will

      ! always exist but in the mean-time we have to do this check:

      ! call assert(associated(remap_axes), 'register_diag_field: remap_axes not set')

      if (associated(remap_axes)) then

        if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then

          active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, remap_axes, &

                     init_time, long_name=long_name, units=units, missing_value=mom_missing_value, &

                     range=range, mask_variant=mask_variant, standard_name=standard_name, &

                     verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                     interp_method=interp_method, tile_count=tile_count, &

                     cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, &

                     cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, &

                     cell_methods=cell_methods, x_cell_method=x_cell_method, &

                     y_cell_method=y_cell_method, v_cell_method=v_cell_method, &

                     conversion=conversion, v_extensive=v_extensive)

          if (active) then

            call diag_remap_set_active(diag_cs%diag_remap_cs(i))

          endif

          module_list = trim(module_list)//","//trim(new_module_name)

          num_modnm = num_modnm + 1

        endif ! remap_axes%needs_remapping

      endif ! associated(remap_axes)

    endif ! axes%rank == 3

  enddo ! i

  ! Register downsampled diagnostics

  do dl=2,max_dsamp_lev

    ! Do not attempt to checksum the downsampled diagnostics

    if (diag_cs%diag_as_chksum) cycle

    new_module_name = trim(module_name)//'_d2'

    axes_d2 => null()

    if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then

      if (axes_in%id == diag_cs%axesTL%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesTL

      elseif (axes_in%id == diag_cs%axesBL%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesBL

      elseif (axes_in%id == diag_cs%axesCuL%id ) then

        axes_d2 => diag_cs%dsamp(dl)%axesCuL

      elseif (axes_in%id == diag_cs%axesCvL%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesCvL

      elseif (axes_in%id == diag_cs%axesTi%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesTi

      elseif (axes_in%id == diag_cs%axesBi%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesBi

      elseif (axes_in%id == diag_cs%axesCui%id ) then

        axes_d2 => diag_cs%dsamp(dl)%axesCui

      elseif (axes_in%id == diag_cs%axesCvi%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesCvi

      elseif (axes_in%id == diag_cs%axesT1%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesT1

      elseif (axes_in%id == diag_cs%axesB1%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesB1

      elseif (axes_in%id == diag_cs%axesCu1%id ) then

        axes_d2 => diag_cs%dsamp(dl)%axesCu1

      elseif (axes_in%id == diag_cs%axesCv1%id) then

        axes_d2 => diag_cs%dsamp(dl)%axesCv1

      else

        !Niki: Should we worry about these, e.g., diag_to_Z_CS?

        call mom_error(warning,"register_diag_field: Could not find a proper axes for " &

              //trim(new_module_name)//"-"//trim(field_name))

      endif

    endif

    ! Register the native diagnostic

    if (associated(axes_d2)) then

      active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes_d2, &

                init_time, long_name=long_name, units=units, missing_value=mom_missing_value, &

                range=range, mask_variant=mask_variant, standard_name=standard_name, &

                verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                interp_method=interp_method, tile_count=tile_count, &

                cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, &

                cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, &

                cell_methods=cell_methods, x_cell_method=x_cell_method, &

                y_cell_method=y_cell_method, v_cell_method=v_cell_method, &

                conversion=conversion, v_extensive=v_extensive)

      module_list = trim(module_list)//","//trim(new_module_name)

      num_modnm = num_modnm + 1

    endif

    ! For each diagnostic coordinate register the diagnostic again under a different module name

    do i=1,diag_cs%num_diag_coords

      new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)//'_d2'

      ! Register diagnostics remapped to z vertical coordinate

      if (axes_in%rank == 3) then

        remap_axes => null()

        if ((axes_in%id == diag_cs%axesTL%id)) then

          remap_axes => diag_cs%dsamp(dl)%remap_axesTL(i)

        elseif (axes_in%id == diag_cs%axesBL%id) then

          remap_axes => diag_cs%dsamp(dl)%remap_axesBL(i)

        elseif (axes_in%id == diag_cs%axesCuL%id ) then

          remap_axes => diag_cs%dsamp(dl)%remap_axesCuL(i)

        elseif (axes_in%id == diag_cs%axesCvL%id) then

          remap_axes => diag_cs%dsamp(dl)%remap_axesCvL(i)

        elseif (axes_in%id == diag_cs%axesTi%id) then

          remap_axes => diag_cs%dsamp(dl)%remap_axesTi(i)

        elseif (axes_in%id == diag_cs%axesBi%id) then

          remap_axes => diag_cs%dsamp(dl)%remap_axesBi(i)

        elseif (axes_in%id == diag_cs%axesCui%id ) then

          remap_axes => diag_cs%dsamp(dl)%remap_axesCui(i)

        elseif (axes_in%id == diag_cs%axesCvi%id) then

          remap_axes => diag_cs%dsamp(dl)%remap_axesCvi(i)

        endif

        ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will

        ! always exist but in the mean-time we have to do this check:

        ! call assert(associated(remap_axes), 'register_diag_field: remap_axes not set')

        if (associated(remap_axes)) then

          if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then

            active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, remap_axes, &

                    init_time, long_name=long_name, units=units, missing_value=mom_missing_value, &

                    range=range, mask_variant=mask_variant, standard_name=standard_name, &

                    verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                    interp_method=interp_method, tile_count=tile_count, &

                    cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, &

                    cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, &

                    cell_methods=cell_methods, x_cell_method=x_cell_method, &

                    y_cell_method=y_cell_method, v_cell_method=v_cell_method, &

                    conversion=conversion, v_extensive=v_extensive)

            if (active) then

              call diag_remap_set_active(diag_cs%diag_remap_cs(i))

            endif

            module_list = trim(module_list)//","//trim(new_module_name)

            num_modnm = num_modnm + 1

          endif ! remap_axes%needs_remapping

        endif ! associated(remap_axes)

      endif ! axes%rank == 3

    enddo ! i

  enddo

  dimensions = ""

  if (axes_in%is_h_point)   dimensions = trim(dimensions)//" xh, yh,"

  if (axes_in%is_q_point)   dimensions = trim(dimensions)//" xq, yq,"

  if (axes_in%is_u_point)   dimensions = trim(dimensions)//" xq, yh,"

  if (axes_in%is_v_point)   dimensions = trim(dimensions)//" xh, yq,"

  if (axes_in%is_layer)     dimensions = trim(dimensions)//" zl,"

  if (axes_in%is_interface) dimensions = trim(dimensions)//" zi,"

  if (len_trim(dimensions) > 0) dimensions = trim_trailing_commas(dimensions)

  if (is_root_pe() .and. (diag_cs%available_diag_doc_unit > 0)) then

    msg = ''

    if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"'

    call attach_cell_methods(-1, axes, cm_string, cell_methods, &

                             x_cell_method, y_cell_method, v_cell_method, &

                             v_extensive=v_extensive)

    module_list = trim(module_list)//"}"

    if (num_modnm <= 1) module_list = module_name

    if (num_varnm <= 1) var_list = ''

    call log_available_diag(dm_id>0, module_list, field_name, cm_string, msg, diag_cs, &

                            long_name, units, standard_name, variants=var_list, dimensions=dimensions)

  endif

  register_diag_field = dm_id

end function register_diag_field

!> Returns True if either the native or CMOR version of the diagnostic were registered. Updates 'dm_id'

!! after calling register_diag_field_expand_axes() for both native and CMOR variants of the field.

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)

  integer,          intent(inout) :: dm_id !< The diag_mediator ID for this diagnostic group

  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" or "ice_shelf_model"

  character(len=*), intent(in) :: field_name !< Name of the diagnostic field

  type(axes_grp),   intent(in) :: axes !< Container with up to 3 integer handles that indicates axes

                                             !! for this field

  type(time_type),  intent(in) :: init_time !< Time at which a field is first available?

  character(len=*), optional, intent(in) :: long_name !< Long name of a field.

  character(len=*), optional, intent(in) :: units !< Units of a field.

  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field

  real,             optional, intent(in) :: missing_value !< A value that indicates missing values in

                                                          !! output files, in unscaled arbitrary units [a]

  real,             optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)

                                                     !! in arbitrary units [a]

  logical,          optional, intent(in) :: mask_variant !< If true a logical mask must be provided

                                                         !! with post_data calls (not used in MOM?)

  logical,          optional, intent(in) :: verbose !< If true, FMS is verbose (not used in MOM?)

  logical,          optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?)

  character(len=*), optional, intent(out):: err_msg !< String into which an error message might be

                                                         !! placed (not used in MOM?)

  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should

                                                         !! not be interpolated as a scalar

  integer,          optional, intent(in) :: tile_count !< no clue (not used in MOM?)

  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field

  character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field

  character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field

  character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field

  character(len=*), optional, intent(in) :: cell_methods !< String to append as cell_methods attribute.

                                                         !! Use '' to have no attribute. If present, this

                                                         !! overrides the default constructed from the default

                                                         !! for each individual axis direction.

  character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction.

                                                         !! Use '' have no method.

  character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction.

                                                         !! Use '' have no method.

  character(len=*), optional, intent(in) :: v_cell_method !< Specifies the cell method for the vertical direction.

                                                         !! Use '' have no method.

  real,             optional, intent(in) :: conversion !< A value to multiply data by before writing to files,

                                                       !! often including factors to undo internal scaling and

                                                       !! in units of [a A-1 ~> 1]

  logical,          optional, intent(in) :: v_extensive !< True for vertically extensive fields (vertically

                                                         !! integrated). Default/absent for intensive.

  ! Local variables

  real :: mom_missing_value ! A value used to indicate missing values in output files, in arbitrary units [a]

  type(diag_ctrl), pointer :: diag_cs => null()

  type(diag_type), pointer :: this_diag => null()

  integer :: fms_id, fms_xyave_id

  character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name, cm_string

  mom_missing_value = axes%diag_cs%missing_value

  if (present(missing_value)) mom_missing_value = missing_value

  register_diag_field_expand_cmor = .false.

  diag_cs => axes%diag_cs

  ! Set up the 'primary' diagnostic, first get an underlying FMS id

  fms_id = register_diag_field_expand_axes(module_name, field_name, axes, init_time, &

             long_name=long_name, units=units, missing_value=mom_missing_value, &

             range=range, mask_variant=mask_variant, standard_name=standard_name, &

             verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

             interp_method=interp_method, tile_count=tile_count)

  if (.not. diag_cs%diag_as_chksum) &

    call attach_cell_methods(fms_id, axes, cm_string, cell_methods, &

                             x_cell_method, y_cell_method, v_cell_method, &

                             v_extensive=v_extensive)

  ! Associated horizontally area-averaged diagnostic

  fms_xyave_id = diag_field_not_found

  if (associated(axes%xyave_axes)) then

    fms_xyave_id = register_diag_field_expand_axes(module_name, trim(field_name)//'_xyave', &

             axes%xyave_axes, init_time, &

             long_name=long_name, units=units, missing_value=mom_missing_value, &

             range=range, mask_variant=mask_variant, standard_name=standard_name, &

             verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

             interp_method=interp_method, tile_count=tile_count)

    if (.not. diag_cs%diag_as_chksum) &

      call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, &

                               cell_methods, v_cell_method, v_extensive=v_extensive)

  endif

  this_diag => null()

  if (fms_id /= diag_field_not_found .or. fms_xyave_id /= diag_field_not_found) then

    call add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name)

    this_diag%fms_xyave_diag_id = fms_xyave_id

    ! Encode and save the cell methods for this diagnostic

    call add_xyz_method(this_diag, axes, x_cell_method, y_cell_method, v_cell_method, v_extensive)

    if (present(v_extensive)) this_diag%v_extensive = v_extensive

    if (present(conversion)) this_diag%conversion_factor = conversion

    register_diag_field_expand_cmor = .true.

  endif

  ! For the CMOR variation of the above diagnostic

  if (present(cmor_field_name) .and. .not. diag_cs%diag_as_chksum) then

    ! Fallback values for strings set to "NULL"

    posted_cmor_units = "not provided"         !

    posted_cmor_standard_name = "not provided" ! Values might be able to be replaced with a CS%missing field?

    posted_cmor_long_name = "not provided"     !

    ! If attributes are present for MOM variable names, use them first for the register_diag_field

    ! call for CMOR verison of the variable

    if (present(units)) posted_cmor_units = units

    if (present(standard_name)) posted_cmor_standard_name = standard_name

    if (present(long_name)) posted_cmor_long_name = long_name

    ! If specified in the call to register_diag_field, override attributes with the CMOR versions

    if (present(cmor_units)) posted_cmor_units = cmor_units

    if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name

    if (present(cmor_long_name)) posted_cmor_long_name = cmor_long_name

    fms_id = register_diag_field_expand_axes(module_name, cmor_field_name, axes, init_time,    &

               long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units),                  &

               missing_value=mom_missing_value, range=range, mask_variant=mask_variant,               &

               standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log, &

               err_msg=err_msg, interp_method=interp_method, tile_count=tile_count)

    call attach_cell_methods(fms_id, axes, cm_string, &

                             cell_methods, x_cell_method, y_cell_method, v_cell_method, &

                             v_extensive=v_extensive)

    ! Associated horizontally area-averaged diagnostic

    fms_xyave_id = diag_field_not_found

    if (associated(axes%xyave_axes)) then

      fms_xyave_id = register_diag_field_expand_axes(module_name, trim(cmor_field_name)//'_xyave', &

               axes%xyave_axes, init_time, &

               long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units),                  &

               missing_value=mom_missing_value, range=range, mask_variant=mask_variant,               &

               standard_name=trim(posted_cmor_standard_name), verbose=verbose, do_not_log=do_not_log, &

               err_msg=err_msg, interp_method=interp_method, tile_count=tile_count)

      call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, &

                               cell_methods, v_cell_method, v_extensive=v_extensive)

    endif

    this_diag => null()

    if (fms_id /= diag_field_not_found .or. fms_xyave_id /= diag_field_not_found) then

      call add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name)

      this_diag%fms_xyave_diag_id = fms_xyave_id

      ! Encode and save the cell methods for this diagnostic

      call add_xyz_method(this_diag, axes, x_cell_method, y_cell_method, v_cell_method, v_extensive)

      if (present(v_extensive)) this_diag%v_extensive = v_extensive

      if (present(conversion)) this_diag%conversion_factor = conversion

      register_diag_field_expand_cmor = .true.

    endif

  endif

end function register_diag_field_expand_cmor

!> Returns an FMS id from register_diag_field_fms (the diag_manager routine) after expanding axes

!! (axes-group) into handles and conditionally adding an FMS area_id for cell_measures.

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)

  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model"

                                              !! or "ice_shelf_model"

  character(len=*), intent(in) :: field_name !< Name of the diagnostic field

  type(axes_grp), target, intent(in) :: axes !< Container with up to 3 integer handles that indicates

                                             !! axes for this field

  type(time_type),  intent(in) :: init_time !< Time at which a field is first available?

  character(len=*), optional, intent(in) :: long_name !< Long name of a field.

  character(len=*), optional, intent(in) :: units !< Units of a field.

  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field

  real,             optional, intent(in) :: missing_value !< A value that indicates missing values in

                                                          !! output files, in unscaled arbitrary units [a]

  real,             optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)

                                                     !! in arbitrary units [a]

  logical,          optional, intent(in) :: mask_variant !< If true a logical mask must be provided

                                                         !! with post_data calls (not used in MOM?)

  logical,          optional, intent(in) :: verbose !< If true, FMS is verbose (not used in MOM?)

  logical,          optional, intent(in) :: do_not_log !< If true, do not log something

                                                       !! (not used in MOM?)

  character(len=*), optional, intent(out):: err_msg !< String into which an error message might be

                                                         !! placed (not used in MOM?)

  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should

                                                         !! not be interpolated as a scalar

  integer,          optional, intent(in) :: tile_count !< no clue (not used in MOM?)

  ! Local variables

  integer :: fms_id, area_id, volume_id

  ! This gets the cell area associated with the grid location of this variable

  area_id = axes%id_area

  volume_id = axes%id_volume

  ! Get the FMS diagnostic id

  if (axes%diag_cs%diag_as_chksum) then

    fms_id = axes%diag_cs%num_chksum_diags + 1

    axes%diag_cs%num_chksum_diags = fms_id

  elseif (present(interp_method) .or. axes%is_h_point) then

    ! If interp_method is provided we must use it

    if (area_id>0) then

      if (volume_id>0) then

        fms_id = register_diag_field_infra(module_name, field_name, axes%handles, &

                 init_time, long_name=long_name, units=units, missing_value=missing_value, &

                 range=range, mask_variant=mask_variant, standard_name=standard_name, &

                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                 interp_method=interp_method, tile_count=tile_count, area=area_id, volume=volume_id)

      else

        fms_id = register_diag_field_infra(module_name, field_name, axes%handles, &

                 init_time, long_name=long_name, units=units, missing_value=missing_value, &

                 range=range, mask_variant=mask_variant, standard_name=standard_name, &

                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                 interp_method=interp_method, tile_count=tile_count, area=area_id)

      endif

    else

      if (volume_id>0) then

        fms_id = register_diag_field_infra(module_name, field_name, axes%handles, &

                 init_time, long_name=long_name, units=units, missing_value=missing_value, &

                 range=range, mask_variant=mask_variant, standard_name=standard_name, &

                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                 interp_method=interp_method, tile_count=tile_count, volume=volume_id)

      else

        fms_id = register_diag_field_infra(module_name, field_name, axes%handles, &

                 init_time, long_name=long_name, units=units, missing_value=missing_value, &

                 range=range, mask_variant=mask_variant, standard_name=standard_name, &

                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                 interp_method=interp_method, tile_count=tile_count)

      endif

    endif

  else

    ! If interp_method is not provided and the field is not at an h-point then interp_method='none'

    if (area_id>0) then

      if (volume_id>0) then

        fms_id = register_diag_field_infra(module_name, field_name, axes%handles, &

                 init_time, long_name=long_name, units=units, missing_value=missing_value, &

                 range=range, mask_variant=mask_variant, standard_name=standard_name, &

                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                 interp_method='none', tile_count=tile_count, area=area_id, volume=volume_id)

      else

        fms_id = register_diag_field_infra(module_name, field_name, axes%handles, &

                 init_time, long_name=long_name, units=units, missing_value=missing_value, &

                 range=range, mask_variant=mask_variant, standard_name=standard_name, &

                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                 interp_method='none', tile_count=tile_count, area=area_id)

      endif

    else

      if (volume_id>0) then

        fms_id = register_diag_field_infra(module_name, field_name, axes%handles, &

                 init_time, long_name=long_name, units=units, missing_value=missing_value, &

                 range=range, mask_variant=mask_variant, standard_name=standard_name, &

                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                 interp_method='none', tile_count=tile_count, volume=volume_id)

      else

        fms_id = register_diag_field_infra(module_name, field_name, axes%handles, &

                 init_time, long_name=long_name, units=units, missing_value=missing_value, &

                 range=range, mask_variant=mask_variant, standard_name=standard_name, &

                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &

                 interp_method='none', tile_count=tile_count)

      endif

    endif

  endif

  register_diag_field_expand_axes = fms_id

end function register_diag_field_expand_axes

!> Create a diagnostic type and attached to list

subroutine add_diag_to_list(diag_cs, dm_id, fms_id, this_diag, axes, module_name, field_name)

  type(diag_ctrl),        pointer       :: diag_cs !< Diagnostics mediator control structure

  integer,                intent(inout) :: dm_id !< The diag_mediator ID for this diagnostic group

  integer,                intent(in)    :: fms_id !< The FMS diag_manager ID for this diagnostic

  type(diag_type),        pointer       :: this_diag !< This diagnostic

  type(axes_grp), target, intent(in)    :: axes !< Container with up to 3 integer handles that

                                                !! indicates axes for this field

  character(len=*),       intent(in)    :: module_name !< Name of this module, usually

                                                       !! "ocean_model" or "ice_shelf_model"

  character(len=*),       intent(in)    :: field_name !< Name of diagnostic

  ! If the diagnostic is needed obtain a diag_mediator ID (if needed)

  if (dm_id == -1) dm_id = get_new_diag_id(diag_cs)

  ! Create a new diag_type to store links in

  call alloc_diag_with_id(dm_id, diag_cs, this_diag)

  call assert(associated(this_diag), 'add_diag_to_list: allocation failed for '//trim(field_name))

  ! Record FMS id, masks and conversion factor, in diag_type

  this_diag%fms_diag_id = fms_id

  this_diag%debug_str = trim(module_name)//"-"//trim(field_name)

  this_diag%axes => axes

end subroutine add_diag_to_list

!> Adds the encoded "cell_methods" for a diagnostics as a diag% property

!! This allows access to the cell_method for a given diagnostics at the time of sending

subroutine add_xyz_method(diag, axes, x_cell_method, y_cell_method, v_cell_method, v_extensive)

  type(diag_type),          pointer       :: diag !< This diagnostic

  type(axes_grp),             intent(in)  :: axes !< Container w/ up to 3 integer handles that indicates

                                                  !! axes for this field

  character(len=*), optional, intent(in)  :: x_cell_method !< Specifies the cell method for the x-direction.

                                                         !! Use '' have no method.

  character(len=*), optional, intent(in)  :: y_cell_method !< Specifies the cell method for the y-direction.

                                                         !! Use '' have no method.

  character(len=*), optional, intent(in)  :: v_cell_method !< Specifies the cell method for the vertical direction.

                                                         !! Use '' have no method.

  logical,          optional, intent(in)  :: v_extensive !< True for vertically extensive fields

                                                         !! (vertically integrated). Default/absent for intensive.

  integer :: xyz_method

  character(len=9) :: mstr

  ! This is a simple way to encode the cell method information made from 3 strings

  ! (x_cell_method,y_cell_method,v_cell_method) in a 3 digit integer xyz

  ! x_cell_method,y_cell_method,v_cell_method can each be 'point' or 'sum' or 'mean'

  ! We can encode these with setting  1 for 'point', 2 for 'sum, 3 for 'mean' in

  ! the 100s position for x, 10s position for y, 1s position for z

  ! E.g., x:sum,y:point,z:mean is 213

  xyz_method = 111

  mstr = diag%axes%v_cell_method

  if (present(v_extensive)) then

    if (present(v_cell_method)) call mom_error(fatal, "attach_cell_methods: " // &

       'Vertical cell method was specified along with the vertically extensive flag.')

    if (v_extensive) then

      mstr='sum'

    else

      mstr='mean'

    endif

  elseif (present(v_cell_method)) then

    mstr = v_cell_method

  endif

  if (trim(mstr)=='sum') then

    xyz_method = xyz_method + 1

  elseif (trim(mstr)=='mean') then

    xyz_method = xyz_method + 2

  endif

  mstr = diag%axes%y_cell_method

  if (present(y_cell_method)) mstr = y_cell_method

  if (trim(mstr)=='sum') then

    xyz_method = xyz_method + 10

  elseif (trim(mstr)=='mean') then

    xyz_method = xyz_method + 20

  endif

  mstr = diag%axes%x_cell_method

  if (present(x_cell_method)) mstr = x_cell_method

  if (trim(mstr)=='sum') then

    xyz_method = xyz_method + 100

  elseif (trim(mstr)=='mean') then

    xyz_method = xyz_method + 200

  endif

  diag%xyz_method = xyz_method

end subroutine add_xyz_method

!> Attaches "cell_methods" attribute to a variable based on defaults for axes_grp or optional arguments.

subroutine attach_cell_methods(id, axes, ostring, cell_methods, &

                               x_cell_method, y_cell_method, v_cell_method, v_extensive)

  integer,                    intent(in)  :: id !< Handle to diagnostic

  type(axes_grp),             intent(in)  :: axes !< Container with up to 3 integer handles that indicates

                                                  !! axes for this field

  character(len=*),           intent(out) :: ostring !< The cell_methods strings that would appear in the file

  character(len=*), optional, intent(in)  :: cell_methods !< String to append as cell_methods attribute.

                                                         !! Use '' to have no attribute. If present, this

                                                         !! overrides the default constructed from the default

                                                         !! for each individual axis direction.

  character(len=*), optional, intent(in)  :: x_cell_method !< Specifies the cell method for the x-direction.

                                                         !! Use '' have no method.

  character(len=*), optional, intent(in)  :: y_cell_method !< Specifies the cell method for the y-direction.

                                                         !! Use '' have no method.

  character(len=*), optional, intent(in)  :: v_cell_method !< Specifies the cell method for the vertical direction.

                                                         !! Use '' have no method.

  logical,          optional, intent(in)  :: v_extensive !< True for vertically extensive fields

                                                         !! (vertically integrated). Default/absent for intensive.

  ! Local variables

  character(len=9) :: axis_name

  logical :: x_mean, y_mean, x_sum, y_sum

  x_mean = .false.

  y_mean = .false.

  x_sum = .false.

  y_sum = .false.

  ostring = ''

  if (present(cell_methods)) then

    if (present(x_cell_method) .or. present(y_cell_method) .or. present(v_cell_method) &

        .or. present(v_extensive)) then

      call mom_error(fatal, "attach_cell_methods: " // &

           'Individual direction cell method was specified along with a "cell_methods" string.')

    endif

    if (len(trim(cell_methods))>0) then

      call mom_diag_field_add_attribute(id, 'cell_methods', trim(cell_methods))

      ostring = trim(cell_methods)

    endif

  else

    if (present(x_cell_method)) then

      if (len(trim(x_cell_method))>0) then

        call get_mom_diag_axis_name(axes%handles(1), axis_name)

        call mom_diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(x_cell_method))

        ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(x_cell_method)

        if (trim(x_cell_method)=='mean') x_mean=.true.

        if (trim(x_cell_method)=='sum') x_sum=.true.

      endif

    else

      if (len(trim(axes%x_cell_method))>0) then

        call get_mom_diag_axis_name(axes%handles(1), axis_name)

        call mom_diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%x_cell_method))

        ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%x_cell_method)

        if (trim(axes%x_cell_method)=='mean') x_mean=.true.

        if (trim(axes%x_cell_method)=='sum') x_sum=.true.

      endif

    endif

    if (present(y_cell_method)) then

      if (len(trim(y_cell_method))>0) then

        call get_mom_diag_axis_name(axes%handles(2), axis_name)

        call mom_diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(y_cell_method))

        ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(y_cell_method)

        if (trim(y_cell_method)=='mean') y_mean=.true.

        if (trim(y_cell_method)=='sum') y_sum=.true.

      endif

    else

      if (len(trim(axes%y_cell_method))>0) then

        call get_mom_diag_axis_name(axes%handles(2), axis_name)

        call mom_diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%y_cell_method))

        ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%y_cell_method)

        if (trim(axes%y_cell_method)=='mean') y_mean=.true.

        if (trim(axes%y_cell_method)=='sum') y_sum=.true.

      endif

    endif

    if (present(v_cell_method)) then

      if (present(v_extensive)) call mom_error(fatal, "attach_cell_methods: " // &

           'Vertical cell method was specified along with the vertically extensive flag.')

      if (len(trim(v_cell_method))>0) then

        if (axes%rank==1) then

          call get_mom_diag_axis_name(axes%handles(1), axis_name)

        elseif (axes%rank==3) then

          call get_mom_diag_axis_name(axes%handles(3), axis_name)

        endif

        call mom_diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(v_cell_method))

        ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(v_cell_method)

      endif

    elseif (present(v_extensive)) then

      if (v_extensive) then

        if (axes%rank==1) then

          call get_mom_diag_axis_name(axes%handles(1), axis_name)

        elseif (axes%rank==3) then

          call get_mom_diag_axis_name(axes%handles(3), axis_name)

        endif

        call mom_diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':sum')

        ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':sum'

      endif

    else

      if (len(trim(axes%v_cell_method))>0) then

        if (axes%rank==1) then

          call get_mom_diag_axis_name(axes%handles(1), axis_name)

        elseif (axes%rank==3) then

          call get_mom_diag_axis_name(axes%handles(3), axis_name)

        endif

        call mom_diag_field_add_attribute(id, 'cell_methods', trim(axis_name)//':'//trim(axes%v_cell_method))

        ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(axes%v_cell_method)

      endif

    endif

    if (x_mean .and. y_mean) then

      call mom_diag_field_add_attribute(id, 'cell_methods', 'area:mean')

      ostring = trim(adjustl(ostring))//' area:mean'

    elseif (x_sum .and. y_sum) then

      call mom_diag_field_add_attribute(id, 'cell_methods', 'area:sum')

      ostring = trim(adjustl(ostring))//' area:sum'

    endif

  endif

  ostring = adjustl(ostring)

end subroutine attach_cell_methods

!> Registers a scalar diagnostic, returning an integer handle

function 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, conversion)

  integer :: register_scalar_field !< An integer handle for a diagnostic array.

  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model"

                                              !! or "ice_shelf_model"

  character(len=*), intent(in) :: field_name !< Name of the diagnostic field

  type(time_type),  intent(in) :: init_time !< Time at which a field is first available?

  type(diag_ctrl),  intent(inout) :: diag_cs !< Structure used to regulate diagnostic output

  character(len=*), optional, intent(in) :: long_name !< Long name of a field.

  character(len=*), optional, intent(in) :: units !< Units of a field.

  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field

  real,             optional, intent(in) :: missing_value !< A value that indicates missing values in

                                                          !! output files, in unscaled arbitrary units [a]

  real,             optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)

                                                     !! in arbitrary units [a]

  logical,          optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?)

  character(len=*), optional, intent(out):: err_msg !< String into which an error message might be

                                                         !! placed (not used in MOM?)

  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not

                                                         !! be interpolated as a scalar

  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field

  character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field

  character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field

  character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field

  real,             optional, intent(in) :: conversion !< A value to multiply data by before writing to files,

                                                       !! often including factors to undo internal scaling and

                                                       !! in units of [a A-1 ~> 1]

  ! Local variables

  real :: mom_missing_value ! A value used to indicate missing values in output files, in arbitrary units [a]

  integer :: dm_id, fms_id

  type(diag_type), pointer :: diag => null(), cmor_diag => null()

  character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name

  character(len=16)  :: dimensions

  mom_missing_value = diag_cs%missing_value

  if (present(missing_value)) mom_missing_value = missing_value

  dm_id = -1

  diag => null()

  cmor_diag => null()

  if (diag_cs%diag_as_chksum) then

    fms_id = diag_cs%num_chksum_diags + 1

    diag_cs%num_chksum_diags = fms_id

  else

    fms_id = register_diag_field_infra(module_name, field_name, init_time, &

                long_name=long_name, units=units, missing_value=mom_missing_value, &

                range=range, standard_name=standard_name, do_not_log=do_not_log, &

                err_msg=err_msg)

  endif

  if (fms_id /= diag_field_not_found) then

    dm_id = get_new_diag_id(diag_cs)

    call alloc_diag_with_id(dm_id, diag_cs, diag)

    call assert(associated(diag), 'register_scalar_field: diag allocation failed')

    diag%fms_diag_id = fms_id

    diag%debug_str = trim(module_name)//"-"//trim(field_name)

    if (present(conversion)) diag%conversion_factor = conversion

  endif

  if (present(cmor_field_name)) then

    ! Fallback values for strings set to "not provided"

    posted_cmor_units = "not provided"

    posted_cmor_standard_name = "not provided"

    posted_cmor_long_name = "not provided"

    ! If attributes are present for MOM variable names, use them as defaults for the

    ! register_diag_field_infra call for CMOR verison of the variable

    if (present(units)) posted_cmor_units = units

    if (present(standard_name)) posted_cmor_standard_name = standard_name

    if (present(long_name)) posted_cmor_long_name = long_name

    ! If specified in the call to register_scalar_field, override attributes with the CMOR versions

    if (present(cmor_units)) posted_cmor_units = cmor_units

    if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name

    if (present(cmor_long_name)) posted_cmor_long_name = cmor_long_name

    fms_id = register_diag_field_infra(module_name, cmor_field_name, init_time, &

           long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), &

           missing_value=mom_missing_value, range=range, &

           standard_name=trim(posted_cmor_standard_name), do_not_log=do_not_log, err_msg=err_msg)

    if (fms_id /= diag_field_not_found) then

      if (dm_id == -1) then

        dm_id = get_new_diag_id(diag_cs)

      endif

      call alloc_diag_with_id(dm_id, diag_cs, cmor_diag)

      cmor_diag%fms_diag_id = fms_id

      cmor_diag%debug_str = trim(module_name)//"-"//trim(cmor_field_name)

      if (present(conversion)) cmor_diag%conversion_factor = conversion

    endif

  endif

  dimensions = "scalar"

  ! Document diagnostics in list of available diagnostics

  if (is_root_pe() .and. diag_cs%available_diag_doc_unit > 0) then

    if (present(cmor_field_name)) then

      call log_available_diag(associated(diag), module_name, field_name, '', '', diag_cs, &

                              long_name, units, standard_name, &

                              variants="{"//trim(field_name)//","//trim(cmor_field_name)//"}", &

                              dimensions=dimensions)

    else

      call log_available_diag(associated(diag), module_name, field_name, '', '', diag_cs, &

                              long_name, units, standard_name, dimensions=dimensions)

    endif

  endif

  register_scalar_field = dm_id

end function register_scalar_field

!> Registers a static diagnostic, returning an integer handle

function 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)

  integer :: register_static_field !< An integer handle for a diagnostic array.

  character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model"

                                              !! or "ice_shelf_model"

  character(len=*), intent(in) :: field_name !< Name of the diagnostic field

  type(axes_grp), target, intent(in) :: axes !< Container with up to 3 integer handles that

                                             !! indicates axes for this field

  character(len=*), optional, intent(in) :: long_name !< Long name of a field.

  character(len=*), optional, intent(in) :: units !< Units of a field.

  character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field

  real,             optional, intent(in) :: missing_value !< A value that indicates missing values in

                                                          !! output files, in unscaled arbitrary units [a]

  real,             optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?)

                                                     !! in arbitrary units [a]

  logical,          optional, intent(in) :: mask_variant !< If true a logical mask must be provided with

                                                         !! post_data calls (not used in MOM?)

  logical,          optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?)

  character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not

                                                         !! be interpolated as a scalar

  integer,          optional, intent(in) :: tile_count !< no clue (not used in MOM?)

  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field

  character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field

  character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field

  character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field

  integer,          optional, intent(in) :: area !< fms_id for area_t

  character(len=*), optional, intent(in) :: x_cell_method !< Specifies the cell method for the x-direction.

  character(len=*), optional, intent(in) :: y_cell_method !< Specifies the cell method for the y-direction.

  character(len=*), optional, intent(in) :: area_cell_method !< Specifies the cell method for area

  real,             optional, intent(in) :: conversion !< A value to multiply data by before writing to files,

                                                       !! often including factors to undo internal scaling and

                                                       !! in units of [a A-1 ~> 1]

  ! Local variables

  real :: mom_missing_value ! A value used to indicate missing values in output files, in arbitrary units [a]

  type(diag_ctrl), pointer :: diag_cs => null() !< A structure that is used to regulate diagnostic output

  type(diag_type), pointer :: diag => null(), cmor_diag => null()

  integer :: dm_id, fms_id

  character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name

  character(len=9) :: axis_name

  character(len=16) :: dimensions

  mom_missing_value = axes%diag_cs%missing_value

  if (present(missing_value)) mom_missing_value = missing_value

  diag_cs => axes%diag_cs

  dm_id = -1

  diag => null()

  cmor_diag => null()

  if (diag_cs%diag_as_chksum) then

    fms_id = diag_cs%num_chksum_diags + 1

    diag_cs%num_chksum_diags = fms_id

  else

    fms_id = register_static_field_infra(module_name, field_name, axes%handles, &

           long_name=long_name, units=units, missing_value=mom_missing_value, &

           range=range, mask_variant=mask_variant, standard_name=standard_name, &

           do_not_log=do_not_log, &

           interp_method=interp_method, tile_count=tile_count, area=area)

  endif

  if (fms_id /= diag_field_not_found) then

    dm_id = get_new_diag_id(diag_cs)

    call alloc_diag_with_id(dm_id, diag_cs, diag)

    call assert(associated(diag), 'register_static_field: diag allocation failed')

    diag%fms_diag_id = fms_id

    diag%debug_str = trim(module_name)//"-"//trim(field_name)

    if (present(conversion)) diag%conversion_factor = conversion

    if (diag_cs%diag_as_chksum) then

      diag%axes => axes

    else

      if (present(x_cell_method)) then

        call get_mom_diag_axis_name(axes%handles(1), axis_name)

        call mom_diag_field_add_attribute(fms_id, 'cell_methods', &

            trim(axis_name)//':'//trim(x_cell_method))

      endif

      if (present(y_cell_method)) then

        call get_mom_diag_axis_name(axes%handles(2), axis_name)

        call mom_diag_field_add_attribute(fms_id, 'cell_methods', &

            trim(axis_name)//':'//trim(y_cell_method))

      endif

      if (present(area_cell_method)) then

        call mom_diag_field_add_attribute(fms_id, 'cell_methods', &

            'area:'//trim(area_cell_method))

      endif

    endif

  endif

  if (present(cmor_field_name) .and. .not. diag_cs%diag_as_chksum) then

    ! Fallback values for strings set to "not provided"

    posted_cmor_units = "not provided"

    posted_cmor_standard_name = "not provided"

    posted_cmor_long_name = "not provided"

    ! If attributes are present for MOM variable names, use them first for the register_static_field

    ! call for CMOR verison of the variable

    if (present(units)) posted_cmor_units = units

    if (present(standard_name)) posted_cmor_standard_name = standard_name

    if (present(long_name)) posted_cmor_long_name = long_name

    ! If specified in the call to register_static_field, override attributes with the CMOR versions

    if (present(cmor_units)) posted_cmor_units = cmor_units

    if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name

    if (present(cmor_long_name)) posted_cmor_long_name = cmor_long_name

    fms_id = register_static_field_infra(module_name, cmor_field_name, axes%handles, &

                long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), &

                missing_value=mom_missing_value, range=range, mask_variant=mask_variant, &

                standard_name=trim(posted_cmor_standard_name), do_not_log=do_not_log, &

                interp_method=interp_method, tile_count=tile_count, area=area)

    if (fms_id /= diag_field_not_found) then

      if (dm_id == -1) then

        dm_id = get_new_diag_id(diag_cs)

      endif

      call alloc_diag_with_id(dm_id, diag_cs, cmor_diag)

      cmor_diag%fms_diag_id = fms_id

      cmor_diag%debug_str = trim(module_name)//"-"//trim(cmor_field_name)

      if (present(conversion)) cmor_diag%conversion_factor = conversion

      if (present(x_cell_method)) then

        call get_mom_diag_axis_name(axes%handles(1), axis_name)

        call mom_diag_field_add_attribute(fms_id, 'cell_methods', trim(axis_name)//':'//trim(x_cell_method))

      endif

      if (present(y_cell_method)) then

        call get_mom_diag_axis_name(axes%handles(2), axis_name)

        call mom_diag_field_add_attribute(fms_id, 'cell_methods', trim(axis_name)//':'//trim(y_cell_method))

      endif

      if (present(area_cell_method)) then

        call mom_diag_field_add_attribute(fms_id, 'cell_methods', 'area:'//trim(area_cell_method))

      endif

    endif

  endif

  dimensions = ""

  if (axes%is_h_point)   dimensions = trim(dimensions)//" xh, yh,"

  if (axes%is_q_point)   dimensions = trim(dimensions)//" xq, yq,"

  if (axes%is_u_point)   dimensions = trim(dimensions)//" xq, yh,"

  if (axes%is_v_point)   dimensions = trim(dimensions)//" xh, yq,"

  if (axes%is_layer)     dimensions = trim(dimensions)//" zl,"

  if (axes%is_interface) dimensions = trim(dimensions)//" zi,"

  if (len_trim(dimensions) > 0) dimensions = trim_trailing_commas(dimensions)

  ! Document diagnostics in list of available diagnostics

  if (is_root_pe() .and. diag_cs%available_diag_doc_unit > 0) then

    if (present(cmor_field_name)) then

      call log_available_diag(associated(diag), module_name, field_name, '', '', diag_cs, &

                              long_name, units, standard_name, &

                              variants="{"//trim(field_name)//","//trim(cmor_field_name)//"}", &

                              dimensions=dimensions)

    else

      call log_available_diag(associated(diag), module_name, field_name, '', '', diag_cs, &

                              long_name, units, standard_name, dimensions=dimensions)

    endif

  endif

  register_static_field = dm_id

end function register_static_field

!> Describe an option setting in the diagnostic files.

subroutine describe_option(opt_name, value, diag_CS)

  character(len=*), intent(in) :: opt_name !< The name of the option

  character(len=*), intent(in) :: value   !< A character string with the setting of the option.

  type(diag_ctrl),  intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  character(len=480) :: mesg

  integer :: len_ind

  len_ind = len_trim(value)  ! Add error handling for long values?

  mesg = "    ! "//trim(opt_name)//": "//trim(value)

  write(diag_cs%available_diag_doc_unit, '(a)') trim(mesg)

end subroutine describe_option

!> Registers a diagnostic using the information encapsulated in the vardesc

!! type argument and returns an integer handle to this diagnostic.  That

!! integer handle is negative if the diagnostic is unused.

function ocean_register_diag(var_desc, G, diag_CS, day)

  integer :: ocean_register_diag  !< An integer handle to this diagnostic.

  type(vardesc),         intent(in) :: var_desc !< The vardesc type describing the diagnostic

  type(ocean_grid_type), intent(in) :: g        !< The ocean's grid type

  type(diag_ctrl), intent(in), target :: diag_cs  !< The diagnostic control structure

  type(time_type),       intent(in) :: day      !< The current model time

  character(len=64) :: var_name         ! A variable's name.

  character(len=48) :: units            ! A variable's units.

  character(len=240) :: longname        ! A variable's longname.

  character(len=8) :: hor_grid, z_grid  ! Variable grid info.

  real :: conversion ! A multiplicative factor for unit conversions for output,

                     ! as might be needed to convert from intensive to extensive

                     ! or for dimensional consistency testing [various] or [a A-1 ~> 1]

  type(axes_grp), pointer :: axes => null()

  call query_vardesc(var_desc, units=units, longname=longname, hor_grid=hor_grid, &

                     z_grid=z_grid, conversion=conversion, caller="ocean_register_diag")

  ! Use the hor_grid and z_grid components of vardesc to determine the

  ! desired axes to register the diagnostic field for.

  select case (z_grid)

    case ("L")

      select case (hor_grid)

        case ("q")  ; axes => diag_cs%axesBL

        case ("h")  ; axes => diag_cs%axesTL

        case ("u")  ; axes => diag_cs%axesCuL

        case ("v")  ; axes => diag_cs%axesCvL

        case ("Bu") ; axes => diag_cs%axesBL

        case ("T")  ; axes => diag_cs%axesTL

        case ("Cu") ; axes => diag_cs%axesCuL

        case ("Cv") ; axes => diag_cs%axesCvL

        case ("z")  ; axes => diag_cs%axeszL

        case default ; call mom_error(fatal, "ocean_register_diag: " // &

                                      "unknown hor_grid component "//trim(hor_grid))

      end select

    case ("i")

      select case (hor_grid)

        case ("q")  ; axes => diag_cs%axesBi

        case ("h")  ; axes => diag_cs%axesTi

        case ("u")  ; axes => diag_cs%axesCui

        case ("v")  ; axes => diag_cs%axesCvi

        case ("Bu") ; axes => diag_cs%axesBi

        case ("T")  ; axes => diag_cs%axesTi

        case ("Cu") ; axes => diag_cs%axesCui

        case ("Cv") ; axes => diag_cs%axesCvi

        case ("z")  ; axes => diag_cs%axeszi

        case default ; call mom_error(fatal, "ocean_register_diag: " // &

                                      "unknown hor_grid component "//trim(hor_grid))

      end select

    case ("1")

      select case (hor_grid)

        case ("q")  ; axes => diag_cs%axesB1

        case ("h")  ; axes => diag_cs%axesT1

        case ("u")  ; axes => diag_cs%axesCu1

        case ("v")  ; axes => diag_cs%axesCv1

        case ("Bu") ; axes => diag_cs%axesB1

        case ("T")  ; axes => diag_cs%axesT1

        case ("Cu") ; axes => diag_cs%axesCu1

        case ("Cv") ; axes => diag_cs%axesCv1

        case default ; call mom_error(fatal, "ocean_register_diag: " // &

                                      "unknown hor_grid component "//trim(hor_grid))

      end select

    case default

      call mom_error(fatal,&

        "ocean_register_diag: unknown z_grid component "//trim(z_grid))

  end select

  ocean_register_diag = register_diag_field("ocean_model", trim(var_name), axes, day, &

          trim(longname), units=trim(units), conversion=conversion, missing_value=-1.0e+34)

end function ocean_register_diag

subroutine diag_mediator_infrastructure_init(err_msg)

  ! This subroutine initializes the FMS diag_manager.

  character(len=*), optional, intent(out)   :: err_msg !< An error message

  call mom_diag_manager_init(err_msg=err_msg)

end subroutine diag_mediator_infrastructure_init

!> diag_mediator_init initializes the MOM diag_mediator and opens the available

!! diagnostics file, if appropriate.

subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)

  type(ocean_grid_type), target, intent(inout) :: g  !< The ocean grid type.

  type(verticalgrid_type), target, intent(in)  :: gv !< The ocean vertical grid structure

  type(unit_scale_type),   target, intent(in)  :: us !< A dimensional unit scaling type

  integer,                    intent(in)    :: nz    !< The number of layers in the model's native grid.

  type(param_file_type),      intent(in)    :: param_file !< Parameter file structure

  type(diag_ctrl),            intent(inout) :: diag_cs !< A pointer to a type with many variables

                                                     !! used for diagnostics

  character(len=*), optional, intent(in)    :: doc_file_dir !< A directory in which to create the

                                                     !! file

  ! This subroutine initializes the diag_mediator and the diag_manager.

  ! The grid type should have its dimensions set by this point, but it

  ! is not necessary that the metrics and axis labels be set up yet.

  ! Local variables

  integer :: ios, i, new_unit

  logical :: opened, new_file

  integer :: remap_answer_date    ! The vintage of the order of arithmetic and expressions to use

                                  ! for remapping.  Values below 20190101 recover the remapping

                                  ! answers from 2018, while higher values use more robust

                                  ! forms of the same remapping expressions.

  integer :: default_answer_date  ! The default setting for the various ANSWER_DATE flags.

  logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for diagnostics

  logical :: dz_diag_needed       ! Logical set True if we need to store dz_begin for reintegrating

  character(len=8)   :: this_pe

  character(len=240) :: doc_file, doc_file_dflt, doc_path

  character(len=240), allocatable :: diag_coords(:)

  ! This include declares and sets the variable "version".

# include "version_variable.h"

  character(len=40) :: mdl = "MOM_diag_mediator" ! This module's name.

  character(len=32) :: filename_appendix = '' ! FMS appendix to filename for ensemble runs

  id_clock_diag_mediator = cpu_clock_id('(Ocean diagnostics framework)', grain=clock_module)

  id_clock_diag_remap = cpu_clock_id('(Ocean diagnostics remapping)', grain=clock_routine)

  id_clock_diag_grid_updates = cpu_clock_id('(Ocean diagnostics grid updates)', grain=clock_routine)

  ! Allocate and initialize list of all diagnostics (and variants)

  allocate(diag_cs%diags(diag_alloc_chunk_size))

  diag_cs%next_free_diag_id = 1

  do i=1, diag_alloc_chunk_size

    call initialize_diag_type(diag_cs%diags(i))

  enddo

  diag_cs%show_call_tree = calltree_showquery()

  ! Read all relevant parameters and write them to the model log.

  call log_version(param_file, mdl, version, "")

  call get_param(param_file, mdl, 'NUM_DIAG_COORDS', diag_cs%num_diag_coords, &

                 'The number of diagnostic vertical coordinates to use. '//&

                 'For each coordinate, an entry in DIAG_COORDS must be provided.', &

                 default=1)

  call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &

                 "This sets the default value for the various _ANSWER_DATE parameters.", &

                 default=99991231)

  call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &

                 do_not_log=.true., default=.true.)

  call get_param(param_file, mdl, "DIAG_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &

                 "If true, use the OM4 remapping-via-subcells algorithm for diagnostics. "//&

                 "See REMAPPING_USE_OM4_SUBCELLS for details. "//&

                 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)

  call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &

                 "The vintage of the expressions and order of arithmetic to use for remapping.  "//&

                 "Values below 20190101 result in the use of older, less accurate expressions "//&

                 "that were in use at the end of 2018.  Higher values result in the use of more "//&

                 "robust and accurate forms of mathematically equivalent expressions.", &

                 default=default_answer_date, do_not_log=.not.gv%Boussinesq)

  if (.not.gv%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)

  call get_param(param_file, mdl, 'USE_INDEX_DIAGNOSTIC_AXES', diag_cs%index_space_axes, &

                 'If true, use a grid index coordinate convention for diagnostic axes. ',&

                 default=.false.)

  dz_diag_needed = .false.

  if (diag_cs%num_diag_coords>0) then

    allocate(diag_coords(diag_cs%num_diag_coords))

    if (diag_cs%num_diag_coords==1) then ! The default is to provide just one instance of Z*

      call get_param(param_file, mdl, 'DIAG_COORDS', diag_coords, &

                 'A list of string tuples associating diag_table modules to '//&

                 'a coordinate definition used for diagnostics. Each string '//&

                 'is of the form "MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME".', &

                 default='z Z ZSTAR')

    else ! If using more than 1 diagnostic coordinate, all must be explicitly defined

      call get_param(param_file, mdl, 'DIAG_COORDS', diag_coords, &

                 'A list of string tuples associating diag_table modules to '//&

                 'a coordinate definition used for diagnostics. Each string '//&

                 'is of the form "MODULE_SUFFIX,PARAMETER_SUFFIX,COORDINATE_NAME".', &

                 fail_if_missing=.true.)

    endif

    allocate(diag_cs%diag_remap_cs(diag_cs%num_diag_coords))

    ! Initialize each diagnostic vertical coordinate

    do i=1, diag_cs%num_diag_coords

      call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i), om4_remap_via_sub_cells, remap_answer_date, gv)

      if (diag_cs%diag_remap_cs(i)%Z_based_coord) dz_diag_needed = .true.

    enddo

    deallocate(diag_coords)

  endif

  call get_param(param_file, mdl, 'DIAG_MISVAL', diag_cs%missing_value, &

                 'Set the default missing value to use for diagnostics.', &

                 units="various", default=1.e20)

  call get_param(param_file, mdl, 'DIAG_AS_CHKSUM', diag_cs%diag_as_chksum, &

                 'Instead of writing diagnostics to the diag manager, write '//&

                 'a text file containing the checksum (bitcount) of the array.',  &

                 default=.false.)

  if (diag_cs%diag_as_chksum) &

    diag_cs%num_chksum_diags = 0

  ! Keep pointers to the grid, h, T, S needed for diagnostic remapping

  diag_cs%G => g

  diag_cs%GV => gv

  diag_cs%US => us

  diag_cs%h => null()

  diag_cs%T => null()

  diag_cs%S => null()

  diag_cs%eqn_of_state => null()

  diag_cs%tv => null()

  allocate(diag_cs%h_begin(g%isd:g%ied,g%jsd:g%jed,nz))

  if (dz_diag_needed) allocate(diag_cs%dz_begin(g%isd:g%ied,g%jsd:g%jed,nz))

#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)

  allocate(diag_cs%h_old(g%isd:g%ied,g%jsd:g%jed,nz))

  diag_cs%h_old(:,:,:) = 0.0

#endif

  diag_cs%is = g%isc - (g%isd-1) ; diag_cs%ie = g%iec - (g%isd-1)

  diag_cs%js = g%jsc - (g%jsd-1) ; diag_cs%je = g%jec - (g%jsd-1)

  diag_cs%isd = g%isd ; diag_cs%ied = g%ied

  diag_cs%jsd = g%jsd ; diag_cs%jed = g%jed

  ! Downsample indices for dl=2 (should be generalized to arbitrary dl, perhaps via a G array)

  diag_cs%dsamp(2)%isc = g%HId2%isc - (g%HId2%isd-1) ; diag_cs%dsamp(2)%iec = g%HId2%iec - (g%HId2%isd-1)

  diag_cs%dsamp(2)%jsc = g%HId2%jsc - (g%HId2%jsd-1) ; diag_cs%dsamp(2)%jec = g%HId2%jec - (g%HId2%jsd-1)

  diag_cs%dsamp(2)%isd = g%HId2%isd ; diag_cs%dsamp(2)%ied = g%HId2%ied

  diag_cs%dsamp(2)%jsd = g%HId2%jsd ; diag_cs%dsamp(2)%jed = g%HId2%jed

  diag_cs%dsamp(2)%isg = g%HId2%isg ; diag_cs%dsamp(2)%ieg = g%HId2%ieg

  diag_cs%dsamp(2)%jsg = g%HId2%jsg ; diag_cs%dsamp(2)%jeg = g%HId2%jeg

  diag_cs%dsamp(2)%isgB = g%HId2%isgB ; diag_cs%dsamp(2)%iegB = g%HId2%iegB

  diag_cs%dsamp(2)%jsgB = g%HId2%jsgB ; diag_cs%dsamp(2)%jegB = g%HId2%jegB

  ! Initialze available diagnostic log file

  if (is_root_pe() .and. (diag_cs%available_diag_doc_unit < 0)) then

    write(this_pe,'(i6.6)') pe_here()

    doc_file_dflt = "available_diags."//this_pe

    call get_param(param_file, mdl, "AVAILABLE_DIAGS_FILE", doc_file, &

                 "A file into which to write a list of all available "//&

                 "ocean diagnostics that can be included in a diag_table.", &

                 default=doc_file_dflt, do_not_log=(diag_cs%available_diag_doc_unit/=-1))

    if (len_trim(doc_file) > 0) then

      new_file = .true. ; if (diag_cs%available_diag_doc_unit /= -1) new_file = .false.

    ! Find an unused unit number.

      do new_unit=512,42,-1

        inquire( new_unit, opened=opened)

        if (.not.opened) exit

      enddo

      if (opened) call mom_error(fatal, &

          "diag_mediator_init failed to find an unused unit number.")

      doc_path = doc_file

      if (present(doc_file_dir)) then ; if (len_trim(doc_file_dir) > 0) then

        doc_path = trim(slasher(doc_file_dir))//trim(doc_file)

      endif ; endif

      diag_cs%available_diag_doc_unit = new_unit

      if (new_file) then

        open(diag_cs%available_diag_doc_unit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', &

             action='WRITE', status='REPLACE', iostat=ios)

      else ! This file is being reopened, and should be appended.

        open(diag_cs%available_diag_doc_unit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', &

             action='WRITE', status='OLD', position='APPEND', iostat=ios)

      endif

      inquire(diag_cs%available_diag_doc_unit, opened=opened)

      if ((.not.opened) .or. (ios /= 0)) then

        call mom_error(fatal, "Failed to open available diags file "//trim(doc_path)//".")

      endif

    endif

  endif

  if (is_root_pe() .and. (diag_cs%chksum_iounit < 0) .and. diag_cs%diag_as_chksum) then

    ! write(this_pe,'(i6.6)') PE_here()

    ! doc_file_dflt = "chksum_diag."//this_pe

    doc_file_dflt = "chksum_diag"

    call get_param(param_file, mdl, "CHKSUM_DIAG_FILE", doc_file, &

                 "A file into which to write all checksums of the "//&

                 "diagnostics listed in the diag_table.", &

                 default=doc_file_dflt, do_not_log=(diag_cs%chksum_iounit/=-1))

    call get_filename_appendix(filename_appendix)

    if (len_trim(filename_appendix) > 0) then

      doc_file = trim(doc_file) //'.'//trim(filename_appendix)

    endif

#ifdef STATSLABEL

    doc_file = trim(doc_file)//"."//trim(adjustl(statslabel))

#endif

    if (len_trim(doc_file) > 0) then

      new_file = .true. ; if (diag_cs%chksum_iounit /= -1) new_file = .false.

    ! Find an unused unit number.

      do new_unit=512,42,-1

        inquire( new_unit, opened=opened)

        if (.not.opened) exit

      enddo

      if (opened) call mom_error(fatal, &

          "diag_mediator_init failed to find an unused unit number.")

      doc_path = doc_file

      if (present(doc_file_dir)) then ; if (len_trim(doc_file_dir) > 0) then

        doc_path = trim(slasher(doc_file_dir))//trim(doc_file)

      endif ; endif

      diag_cs%chksum_iounit = new_unit

      if (new_file) then

        open(diag_cs%chksum_iounit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', &

             action='WRITE', status='REPLACE', iostat=ios)

      else ! This file is being reopened, and should be appended.

        open(diag_cs%chksum_iounit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', &

             action='WRITE', status='OLD', position='APPEND', iostat=ios)

      endif

      inquire(diag_cs%chksum_iounit, opened=opened)

      if ((.not.opened) .or. (ios /= 0)) then

        call mom_error(fatal, "Failed to open checksum diags file "//trim(doc_path)//".")

      endif

    endif

  endif

end subroutine diag_mediator_init

!> Set pointers to the default state fields used to remap diagnostics.

subroutine diag_set_state_ptrs(h, tv, diag_cs)

  real, dimension(:,:,:), target, intent(in   ) :: h !< the model thickness array [H ~> m or kg m-2]

  type(thermo_var_ptrs),  target, intent(in   ) :: tv !< A structure with thermodynamic variables that are

                                                      !! are used to convert thicknesses to vertical extents

  type(diag_ctrl),                intent(inout) :: diag_cs !< diag mediator control structure

  ! Keep pointers to h, T, S needed for the diagnostic remapping

  diag_cs%h => h

  diag_cs%T => tv%T

  diag_cs%S => tv%S

  diag_cs%eqn_of_state => tv%eqn_of_state

  diag_cs%tv => tv

end subroutine

!> Build/update vertical grids for diagnostic remapping.

!! \note The target grids need to be updated whenever sea surface

!! height changes.

subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensive, update_extensive )

  type(diag_ctrl),        intent(inout) :: diag_cs      !< Diagnostics control structure

  real, target, optional, intent(in   ) :: alt_h(:,:,:) !< Used if remapped grids should be something other than

                                                        !! the current thicknesses [H ~> m or kg m-2]

  real, target, optional, intent(in   ) :: alt_t(:,:,:) !< Used if remapped grids should be something other than

                                                        !! the current temperatures [C ~> degC]

  real, target, optional, intent(in   ) :: alt_s(:,:,:) !< Used if remapped grids should be something other than

                                                        !! the current salinity [S ~> ppt]

  logical, optional,      intent(in   ) :: update_intensive !< If true (default), update the grids used for

                                                            !! intensive diagnostics

  logical, optional,      intent(in   ) :: update_extensive !< If true (not default), update the grids used for

                                                            !! intensive diagnostics

  ! Local variables

  integer :: m

  real, dimension(:,:,:), pointer :: h_diag => null() ! The layer thicknesses for diagnostics [H ~> m or kg m-2]

  real, dimension(:,:,:), pointer :: t_diag => null() ! The layer temperatures for diagnostics [C ~> degC]

  real, dimension(:,:,:), pointer :: s_diag => null() ! The layer salinities for diagnostics [S ~> ppt]

  real, dimension(diag_cs%G%isd:diag_cS%G%ied, diag_cs%G%jsd:diag_cS%G%jed, diag_cs%GV%ke) :: &

    dz_diag     ! Layer vertical extents for remapping [Z ~> m]

  logical :: update_intensive_local, update_extensive_local, dz_diag_needed

  if (diag_cs%show_call_tree) call calltree_enter("diag_update_remap_grids()")

  ! Set values based on optional input arguments

  if (present(alt_h)) then

    h_diag => alt_h

  else

    h_diag => diag_cs%h

  endif

  if (present(alt_t)) then

    t_diag => alt_t

  else

    t_diag => diag_cs%T

  endif

  if (present(alt_s)) then

    s_diag => alt_s

  else

    s_diag => diag_cs%S

  endif

  ! Defaults here are based on wanting to update intensive quantities frequently as soon as the model state changes.

  ! Conversely, for extensive quantities, in an effort to close budgets and to be consistent with the total time

  ! tendency, we construct the diagnostic grid at the beginning of the baroclinic timestep and remap all extensive

  ! quantities to the same grid

  update_intensive_local = .true.

  if (present(update_intensive)) update_intensive_local = update_intensive

  update_extensive_local = .false.

  if (present(update_extensive)) update_extensive_local = update_extensive

  if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates)

  if (diag_cs%diag_grid_overridden) then

    call mom_error(fatal, "diag_update_remap_grids was called, but current grids in "// &

                          "diagnostic structure have been overridden")

  endif

  ! Determine the diagnostic grid spacing in height units, if it is needed.

  dz_diag_needed = .false.

  if (update_intensive_local .or. update_extensive_local) then

    do m=1, diag_cs%num_diag_coords

      if (diag_cs%diag_remap_cs(m)%Z_based_coord) dz_diag_needed = .true.

    enddo

  endif

  if (dz_diag_needed) then

    call thickness_to_dz(h_diag, diag_cs%tv, dz_diag, diag_cs%G, diag_cs%GV, diag_cs%US, halo_size=1)

  endif

  if (update_intensive_local) then

    do m=1, diag_cs%num_diag_coords

      if (diag_cs%diag_remap_cs(m)%Z_based_coord) then

        call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, dz_diag, t_diag, s_diag, &

                               diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h)

      else

        call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, t_diag, s_diag, &

                               diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h)

      endif

    enddo

  endif

  if (update_extensive_local) then

    diag_cs%h_begin(:,:,:) = diag_cs%h(:,:,:)

    if (dz_diag_needed) diag_cs%dz_begin(:,:,:) = dz_diag(:,:,:)

    do m=1, diag_cs%num_diag_coords

      if (diag_cs%diag_remap_cs(m)%Z_based_coord) then

        call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, dz_diag, t_diag, s_diag, &

                               diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h_extensive)

      else

        call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, t_diag, s_diag, &

                               diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h_extensive)

      endif

    enddo

  endif

#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)

  ! Keep a copy of H - used to check whether grids are up-to-date

  ! when doing remapping.

  diag_cs%h_old(:,:,:) = diag_cs%h(:,:,:)

#endif

  if (id_clock_diag_grid_updates>0) call cpu_clock_end(id_clock_diag_grid_updates)

  if (diag_cs%show_call_tree) call calltree_leave("diag_update_remap_grids()")

end subroutine diag_update_remap_grids

!> Sets up the 2d and 3d masks for native diagnostics

subroutine diag_masks_set(G, nz, diag_cs)

  type(ocean_grid_type), target, intent(in) :: g  !< The ocean grid type.

  integer,                       intent(in) :: nz !< The number of layers in the model's native grid.

  type(diag_ctrl),               pointer    :: diag_cs !< A pointer to a type with many variables

                                                       !! used for diagnostics

  ! Local variables

  integer :: k

  ! 2d masks point to the model masks since they are identical

  diag_cs%mask2dT  => g%mask2dT

  diag_cs%mask2dBu => g%mask2dBu

  diag_cs%mask2dCu => g%mask2dCu

  diag_cs%mask2dCv => g%mask2dCv

  ! 3d native masks are needed by diag_manager but the native variables

  ! can only be masked 2d - for ocean points, all layers exists.

  allocate(diag_cs%mask3dTL(g%isd:g%ied,g%jsd:g%jed,1:nz))

  allocate(diag_cs%mask3dBL(g%IsdB:g%IedB,g%JsdB:g%JedB,1:nz))

  allocate(diag_cs%mask3dCuL(g%IsdB:g%IedB,g%jsd:g%jed,1:nz))

  allocate(diag_cs%mask3dCvL(g%isd:g%ied,g%JsdB:g%JedB,1:nz))

  do k=1,nz

    diag_cs%mask3dTL(:,:,k) = diag_cs%mask2dT(:,:)

    diag_cs%mask3dBL(:,:,k) = diag_cs%mask2dBu(:,:)

    diag_cs%mask3dCuL(:,:,k) = diag_cs%mask2dCu(:,:)

    diag_cs%mask3dCvL(:,:,k) = diag_cs%mask2dCv(:,:)

  enddo

  allocate(diag_cs%mask3dTi(g%isd:g%ied,g%jsd:g%jed,1:nz+1))

  allocate(diag_cs%mask3dBi(g%IsdB:g%IedB,g%JsdB:g%JedB,1:nz+1))

  allocate(diag_cs%mask3dCui(g%IsdB:g%IedB,g%jsd:g%jed,1:nz+1))

  allocate(diag_cs%mask3dCvi(g%isd:g%ied,g%JsdB:g%JedB,1:nz+1))

  do k=1,nz+1

    diag_cs%mask3dTi(:,:,k) = diag_cs%mask2dT(:,:)

    diag_cs%mask3dBi(:,:,k) = diag_cs%mask2dBu(:,:)

    diag_cs%mask3dCui(:,:,k) = diag_cs%mask2dCu(:,:)

    diag_cs%mask3dCvi(:,:,k) = diag_cs%mask2dCv(:,:)

  enddo

  ! Allocate and initialize the downsampled masks

  call downsample_diag_masks_set(g, nz, diag_cs)

end subroutine diag_masks_set

!> Set the extents and fill values for the piecemeal buffers for all axes

subroutine set_piecemeal_extents(diag_cs)

  type(diag_ctrl), intent(inout) :: diag_cs !< A pointer to a type with many variables

                                                       !! used for diagnostics

  ! Piecemeal buffers for 2d axes

  call diag_cs%axesT1%piecemeal_2d%set_extents_from_array(diag_cs%mask2dT, diag_cs%missing_value)

  call diag_cs%axesB1%piecemeal_2d%set_extents_from_array(diag_cs%mask2dBu, diag_cs%missing_value)

  call diag_cs%axesCu1%piecemeal_2d%set_extents_from_array(diag_cs%mask2dCu, diag_cs%missing_value)

  call diag_cs%axesCv1%piecemeal_2d%set_extents_from_array(diag_cs%mask2dCv, diag_cs%missing_value)

  ! Piecemeal buffers for 3d axes

  call diag_cs%axesTL%piecemeal_3d%set_extents_from_array(diag_cs%mask3dTL, diag_cs%missing_value)

  call diag_cs%axesBL%piecemeal_3d%set_extents_from_array(diag_cs%mask3dBL, diag_cs%missing_value)

  call diag_cs%axesCuL%piecemeal_3d%set_extents_from_array(diag_cs%mask3dCuL, diag_cs%missing_value)

  call diag_cs%axesCvL%piecemeal_3d%set_extents_from_array(diag_cs%mask3dCvL, diag_cs%missing_value)

  call diag_cs%axesTi%piecemeal_3d%set_extents_from_array(diag_cs%mask3dTi, diag_cs%missing_value)

  call diag_cs%axesBi%piecemeal_3d%set_extents_from_array(diag_cs%mask3dBi, diag_cs%missing_value)

  call diag_cs%axesCui%piecemeal_3d%set_extents_from_array(diag_cs%mask3dCui, diag_cs%missing_value)

  call diag_cs%axesCvi%piecemeal_3d%set_extents_from_array(diag_cs%mask3dCvi, diag_cs%missing_value)

end subroutine set_piecemeal_extents

subroutine diag_mediator_close_registration(diag_CS)

  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output

  integer :: i

  if (diag_cs%available_diag_doc_unit > -1) then

    close(diag_cs%available_diag_doc_unit) ; diag_cs%available_diag_doc_unit = -2

  endif

  do i=1, diag_cs%num_diag_coords

    call diag_remap_diag_registration_closed(diag_cs%diag_remap_cs(i))

  enddo

end subroutine diag_mediator_close_registration

subroutine axes_grp_end(axes)

  type(axes_grp), intent(inout) :: axes   !< Axes group to be destroyed

  deallocate(axes%handles)

  if (associated(axes%mask2d)) deallocate(axes%mask2d)

  if (associated(axes%mask3d)) deallocate(axes%mask3d)

end subroutine axes_grp_end

subroutine diag_mediator_end(time, diag_CS, end_diag_manager)

  type(time_type),   intent(in)  :: time !< The current model time

  type(diag_ctrl), intent(inout) :: diag_cs !< Structure used to regulate diagnostic output

  logical, optional, intent(in)  :: end_diag_manager !< If true, call diag_manager_end()

  ! Local variables

  type(diag_type), pointer :: diag, next_diag

  integer :: i, dl

  if (diag_cs%available_diag_doc_unit > -1) then

    close(diag_cs%available_diag_doc_unit) ; diag_cs%available_diag_doc_unit = -3

  endif

  if (diag_cs%chksum_iounit > -1) then

    close(diag_cs%chksum_iounit) ; diag_cs%chksum_iounit = -3

  endif

  do i=1, diag_cs%next_free_diag_id - 1

    if (associated(diag_cs%diags(i)%next)) then

      next_diag => diag_cs%diags(i)%next

      do while (associated(next_diag))

        diag => next_diag

        next_diag => diag%next

        deallocate(diag)

      enddo

    endif

  enddo

  deallocate(diag_cs%diags)

  do i=1, diag_cs%num_diag_coords

    call diag_remap_end(diag_cs%diag_remap_cs(i))

  enddo

  call diag_grid_storage_end(diag_cs%diag_grid_temp)

  if (associated(diag_cs%mask3dTL))  deallocate(diag_cs%mask3dTL)

  if (associated(diag_cs%mask3dBL))  deallocate(diag_cs%mask3dBL)

  if (associated(diag_cs%mask3dCuL)) deallocate(diag_cs%mask3dCuL)

  if (associated(diag_cs%mask3dCvL)) deallocate(diag_cs%mask3dCvL)

  if (associated(diag_cs%mask3dTi))  deallocate(diag_cs%mask3dTi)

  if (associated(diag_cs%mask3dBi))  deallocate(diag_cs%mask3dBi)

  if (associated(diag_cs%mask3dCui)) deallocate(diag_cs%mask3dCui)

  if (associated(diag_cs%mask3dCvi)) deallocate(diag_cs%mask3dCvi)

  do dl=2,max_dsamp_lev

    if (associated(diag_cs%dsamp(dl)%mask2dT))   deallocate(diag_cs%dsamp(dl)%mask2dT)

    if (associated(diag_cs%dsamp(dl)%mask2dBu))  deallocate(diag_cs%dsamp(dl)%mask2dBu)

    if (associated(diag_cs%dsamp(dl)%mask2dCu))  deallocate(diag_cs%dsamp(dl)%mask2dCu)

    if (associated(diag_cs%dsamp(dl)%mask2dCv))  deallocate(diag_cs%dsamp(dl)%mask2dCv)

    if (associated(diag_cs%dsamp(dl)%mask3dTL))  deallocate(diag_cs%dsamp(dl)%mask3dTL)

    if (associated(diag_cs%dsamp(dl)%mask3dBL))  deallocate(diag_cs%dsamp(dl)%mask3dBL)

    if (associated(diag_cs%dsamp(dl)%mask3dCuL)) deallocate(diag_cs%dsamp(dl)%mask3dCuL)

    if (associated(diag_cs%dsamp(dl)%mask3dCvL)) deallocate(diag_cs%dsamp(dl)%mask3dCvL)

    if (associated(diag_cs%dsamp(dl)%mask3dTi))  deallocate(diag_cs%dsamp(dl)%mask3dTi)

    if (associated(diag_cs%dsamp(dl)%mask3dBi))  deallocate(diag_cs%dsamp(dl)%mask3dBi)

    if (associated(diag_cs%dsamp(dl)%mask3dCui)) deallocate(diag_cs%dsamp(dl)%mask3dCui)

    if (associated(diag_cs%dsamp(dl)%mask3dCvi)) deallocate(diag_cs%dsamp(dl)%mask3dCvi)

    do i=1,diag_cs%num_diag_coords

      if (associated(diag_cs%dsamp(dl)%remap_axesTL(i)%dsamp(dl)%mask3d)) &

        deallocate(diag_cs%dsamp(dl)%remap_axesTL(i)%dsamp(dl)%mask3d)

      if (associated(diag_cs%dsamp(dl)%remap_axesCuL(i)%dsamp(dl)%mask3d)) &

        deallocate(diag_cs%dsamp(dl)%remap_axesCuL(i)%dsamp(dl)%mask3d)

      if (associated(diag_cs%dsamp(dl)%remap_axesCvL(i)%dsamp(dl)%mask3d)) &

        deallocate(diag_cs%dsamp(dl)%remap_axesCvL(i)%dsamp(dl)%mask3d)

      if (associated(diag_cs%dsamp(dl)%remap_axesBL(i)%dsamp(dl)%mask3d)) &

        deallocate(diag_cs%dsamp(dl)%remap_axesBL(i)%dsamp(dl)%mask3d)

      if (associated(diag_cs%dsamp(dl)%remap_axesTi(i)%dsamp(dl)%mask3d)) &

        deallocate(diag_cs%dsamp(dl)%remap_axesTi(i)%dsamp(dl)%mask3d)

      if (associated(diag_cs%dsamp(dl)%remap_axesCui(i)%dsamp(dl)%mask3d)) &

        deallocate(diag_cs%dsamp(dl)%remap_axesCui(i)%dsamp(dl)%mask3d)

      if (associated(diag_cs%dsamp(dl)%remap_axesCvi(i)%dsamp(dl)%mask3d)) &

        deallocate(diag_cs%dsamp(dl)%remap_axesCvi(i)%dsamp(dl)%mask3d)

      if (associated(diag_cs%dsamp(dl)%remap_axesBi(i)%dsamp(dl)%mask3d)) &

        deallocate(diag_cs%dsamp(dl)%remap_axesBi(i)%dsamp(dl)%mask3d)

    enddo

  enddo

  ! axes_grp masks may point to diag_cs masks, so do these after mask dealloc

  do i=1, diag_cs%num_diag_coords

    call axes_grp_end(diag_cs%remap_axesZL(i))

    call axes_grp_end(diag_cs%remap_axesZi(i))

    call axes_grp_end(diag_cs%remap_axesTL(i))

    call axes_grp_end(diag_cs%remap_axesTi(i))

    call axes_grp_end(diag_cs%remap_axesBL(i))

    call axes_grp_end(diag_cs%remap_axesBi(i))

    call axes_grp_end(diag_cs%remap_axesCuL(i))

    call axes_grp_end(diag_cs%remap_axesCui(i))

    call axes_grp_end(diag_cs%remap_axesCvL(i))

    call axes_grp_end(diag_cs%remap_axesCvi(i))

  enddo

  if (diag_cs%num_diag_coords > 0) then

    deallocate(diag_cs%remap_axesZL)

    deallocate(diag_cs%remap_axesZi)

    deallocate(diag_cs%remap_axesTL)

    deallocate(diag_cs%remap_axesTi)

    deallocate(diag_cs%remap_axesBL)

    deallocate(diag_cs%remap_axesBi)

    deallocate(diag_cs%remap_axesCuL)

    deallocate(diag_cs%remap_axesCui)

    deallocate(diag_cs%remap_axesCvL)

    deallocate(diag_cs%remap_axesCvi)

  endif

  do dl=2,max_dsamp_lev

    if (allocated(diag_cs%dsamp(dl)%remap_axesTL)) &

      deallocate(diag_cs%dsamp(dl)%remap_axesTL)

    if (allocated(diag_cs%dsamp(dl)%remap_axesTi)) &

      deallocate(diag_cs%dsamp(dl)%remap_axesTi)

    if (allocated(diag_cs%dsamp(dl)%remap_axesBL)) &

      deallocate(diag_cs%dsamp(dl)%remap_axesBL)

    if (allocated(diag_cs%dsamp(dl)%remap_axesBi)) &

      deallocate(diag_cs%dsamp(dl)%remap_axesBi)

    if (allocated(diag_cs%dsamp(dl)%remap_axesCuL)) &

      deallocate(diag_cs%dsamp(dl)%remap_axesCuL)

    if (allocated(diag_cs%dsamp(dl)%remap_axesCui)) &

      deallocate(diag_cs%dsamp(dl)%remap_axesCui)

    if (allocated(diag_cs%dsamp(dl)%remap_axesCvL)) &

      deallocate(diag_cs%dsamp(dl)%remap_axesCvL)

    if (allocated(diag_cs%dsamp(dl)%remap_axesCvi)) &

      deallocate(diag_cs%dsamp(dl)%remap_axesCvi)

  enddo

#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)

  deallocate(diag_cs%h_old)

#endif

  if (present(end_diag_manager)) then

    if (end_diag_manager) call mom_diag_manager_end(time)

  endif

end subroutine diag_mediator_end

!> Returns a new diagnostic id, it may be necessary to expand the diagnostics array.

integer function get_new_diag_id(diag_cs)

  type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure

  ! Local variables

  type(diag_type), dimension(:), allocatable :: tmp

  integer :: i

  if (diag_cs%next_free_diag_id > size(diag_cs%diags)) then

    call assert(diag_cs%next_free_diag_id - size(diag_cs%diags) == 1, &

                'get_new_diag_id: inconsistent diag id')

    ! Increase the size of diag_cs%diags and copy data over.

    ! Do not use move_alloc() because it is not supported by Fortran 90

    allocate(tmp(size(diag_cs%diags)))

    tmp(:) = diag_cs%diags(:)

    deallocate(diag_cs%diags)

    allocate(diag_cs%diags(size(tmp) + diag_alloc_chunk_size))

    diag_cs%diags(1:size(tmp)) = tmp(:)

    deallocate(tmp)

    ! Initialize new part of the diag array.

    do i=diag_cs%next_free_diag_id, size(diag_cs%diags)

      call initialize_diag_type(diag_cs%diags(i))

    enddo

  endif

  get_new_diag_id = diag_cs%next_free_diag_id

  diag_cs%next_free_diag_id = diag_cs%next_free_diag_id + 1

end function get_new_diag_id

!> Initializes a diag_type (used after allocating new memory)

subroutine initialize_diag_type(diag)

  type(diag_type), intent(inout) :: diag !< diag_type to be initialized

  diag%in_use = .false.

  diag%fms_diag_id = -1

  diag%axes => null()

  diag%next => null()

  diag%conversion_factor = 0.

end subroutine initialize_diag_type

!> Make a new diagnostic. Either use memory which is in the array of 'primary'

!! diagnostics, or if that is in use, insert it to the list of secondary diagnostics.

subroutine alloc_diag_with_id(diag_id, diag_cs, diag)

  integer,                 intent(in   ) :: diag_id !< id for the diagnostic

  type(diag_ctrl), target, intent(inout) :: diag_cs !< structure used to regulate diagnostic output

  type(diag_type),         pointer       :: diag    !< structure representing a diagnostic (inout)

  type(diag_type), pointer :: tmp => null()

  if (.not. diag_cs%diags(diag_id)%in_use) then

    diag => diag_cs%diags(diag_id)

  else

    allocate(diag)

    tmp => diag_cs%diags(diag_id)%next

    diag_cs%diags(diag_id)%next => diag

    diag%next => tmp

  endif

  diag%in_use = .true.

end subroutine alloc_diag_with_id

!> Log a diagnostic to the available diagnostics file.

subroutine log_available_diag(used, module_name, field_name, cell_methods_string, comment, &

                              diag_CS, long_name, units, standard_name, variants, dimensions)

  logical,          intent(in) :: used !< Whether this diagnostic was in the diag_table or not

  character(len=*), intent(in) :: module_name !< Name of the diagnostic module

  character(len=*), intent(in) :: field_name !< Name of this diagnostic field

  character(len=*), intent(in) :: cell_methods_string !< The spatial component of the CF cell_methods attribute

  character(len=*), intent(in) :: comment !< A comment to append after [Used|Unused]

  type(diag_ctrl),  intent(in) :: diag_CS  !< The diagnostics control structure

  character(len=*), optional, intent(in) :: dimensions !< Descriptor of the horizontal and vertical dimensions

  character(len=*), optional, intent(in) :: long_name !< CF long name of diagnostic

  character(len=*), optional, intent(in) :: units !< Units for diagnostic

  character(len=*), optional, intent(in) :: standard_name !< CF standardized name of diagnostic

  character(len=*), optional, intent(in) :: variants !< Alternate modules and variable names for

                                                     !! this diagnostic and derived diagnostics

  ! Local variables

  character(len=240) :: mesg

  if (used) then

    mesg = '"'//trim(field_name)//'"  [Used]'

  else

    mesg = '"'//trim(field_name)//'"  [Unused]'

  endif

  if (len(trim((comment)))>0) then

    write(diag_cs%available_diag_doc_unit, '(a,1x,"(",a,")")') trim(mesg),trim(comment)

  else

    write(diag_cs%available_diag_doc_unit, '(a)') trim(mesg)

  endif

  call describe_option("modules", module_name, diag_cs)

  if (present(dimensions)) then

    if (len(trim(dimensions)) > 0) then

      call describe_option("dimensions", dimensions, diag_cs)

    endif

  endif

  if (present(long_name)) call describe_option("long_name", long_name, diag_cs)

  if (present(units)) call describe_option("units", units, diag_cs)

  if (present(standard_name)) &

    call describe_option("standard_name", standard_name, diag_cs)

  if (len(trim((cell_methods_string)))>0) &

    call describe_option("cell_methods", trim(cell_methods_string), diag_cs)

  if (present(variants)) then ; if (len(trim(variants)) > 0) then

    call describe_option("variants", variants, diag_cs)

  endif ; endif

end subroutine log_available_diag

!> Log the diagnostic chksum to the chksum diag file

subroutine log_chksum_diag(docunit, description, chksum)

  integer,          intent(in) :: docunit     !< Handle of the log file

  character(len=*), intent(in) :: description !< Name of the diagnostic module

  integer,          intent(in) :: chksum      !< chksum of the diagnostic

  write(docunit, '(a,1x,i9.8)') description, chksum

  flush(docunit)

end subroutine log_chksum_diag

!> Allocates fields necessary to store diagnostic remapping fields

subroutine diag_grid_storage_init(grid_storage, G, GV, diag)

  type(diag_grid_storage), intent(inout) :: grid_storage !< Structure containing a snapshot of the target grids

  type(ocean_grid_type),   intent(in)    :: g           !< Horizontal grid

  type(verticalgrid_type), intent(in)    :: gv          !< ocean vertical grid structure

  type(diag_ctrl),         intent(in)    :: diag        !< Diagnostic control structure used as the constructor

                                                        !! template for this routine

  integer :: m, nz

  grid_storage%num_diag_coords = diag%num_diag_coords

  ! Don't do anything else if there are no remapped coordinates

  if (grid_storage%num_diag_coords < 1) return

  ! Allocate memory for the native space

  allocate( grid_storage%h_state(g%isd:g%ied, g%jsd:g%jed, gv%ke))

  ! Allocate diagnostic remapping structures

  allocate(grid_storage%diag_grids(diag%num_diag_coords))

  ! Loop through and allocate memory for the grid on each target coordinate

  do m = 1, diag%num_diag_coords

    nz = diag%diag_remap_cs(m)%nz

    allocate(grid_storage%diag_grids(m)%h(g%isd:g%ied,g%jsd:g%jed, nz))

  enddo

end subroutine diag_grid_storage_init

!> Copy from the main diagnostic arrays to the grid storage as well as the native thicknesses

subroutine diag_copy_diag_to_storage(grid_storage, h_state, diag)

  type(diag_grid_storage), intent(inout) :: grid_storage !< Structure containing a snapshot of the target grids

  real, dimension(:,:,:),  intent(in)    :: h_state     !< Current model thicknesses [H ~> m or kg m-2]

  type(diag_ctrl),         intent(in)    :: diag     !< Diagnostic control structure used as the constructor

  integer :: m

  ! Don't do anything else if there are no remapped coordinates

  if (grid_storage%num_diag_coords < 1) return

  grid_storage%h_state(:,:,:) = h_state(:,:,:)

  do m = 1,grid_storage%num_diag_coords

    if (diag%diag_remap_cs(m)%nz > 0) &

      grid_storage%diag_grids(m)%h(:,:,:) = diag%diag_remap_cs(m)%h(:,:,:)

  enddo

end subroutine diag_copy_diag_to_storage

!> Copy from the stored diagnostic arrays to the main diagnostic grids

subroutine diag_copy_storage_to_diag(diag, grid_storage)

  type(diag_ctrl),         intent(inout) :: diag     !< Diagnostic control structure used as the constructor

  type(diag_grid_storage), intent(in)    :: grid_storage !< Structure containing a snapshot of the target grids

  integer :: m

  ! Don't do anything else if there are no remapped coordinates

  if (grid_storage%num_diag_coords < 1) return

  diag%diag_grid_overridden = .true.

  do m = 1,grid_storage%num_diag_coords

    if (diag%diag_remap_cs(m)%nz > 0) &

      diag%diag_remap_cs(m)%h(:,:,:) = grid_storage%diag_grids(m)%h(:,:,:)

  enddo

end subroutine diag_copy_storage_to_diag

!> Save the current diagnostic grids in the temporary structure within diag

subroutine diag_save_grids(diag)

  type(diag_ctrl),         intent(inout) :: diag     !< Diagnostic control structure used as the constructor

  integer :: m

  ! Don't do anything else if there are no remapped coordinates

  if (diag%num_diag_coords < 1) return

  do m = 1,diag%num_diag_coords

    if (diag%diag_remap_cs(m)%nz > 0) &

      diag%diag_grid_temp%diag_grids(m)%h(:,:,:) = diag%diag_remap_cs(m)%h(:,:,:)

  enddo

end subroutine diag_save_grids

!> Restore the diagnostic grids from the temporary structure within diag

subroutine diag_restore_grids(diag)

  type(diag_ctrl),         intent(inout) :: diag     !< Diagnostic control structure used as the constructor

  integer :: m

  ! Don't do anything else if there are no remapped coordinates

  if (diag%num_diag_coords < 1) return

  diag%diag_grid_overridden = .false.

  do m = 1,diag%num_diag_coords

    if (diag%diag_remap_cs(m)%nz > 0) &

      diag%diag_remap_cs(m)%h(:,:,:) = diag%diag_grid_temp%diag_grids(m)%h(:,:,:)

  enddo

end subroutine diag_restore_grids

!> Deallocates the fields in the remapping fields container

subroutine diag_grid_storage_end(grid_storage)

  type(diag_grid_storage), intent(inout) :: grid_storage !< Structure containing a snapshot of the target grids

  ! Local variables

  integer :: m

  ! Don't do anything else if there are no remapped coordinates

  if (grid_storage%num_diag_coords < 1) return

  ! Deallocate memory for the native space

  deallocate(grid_storage%h_state)

  ! Loop through and deallocate memory for the grid on each target coordinate

  do m = 1, grid_storage%num_diag_coords

    deallocate(grid_storage%diag_grids(m)%h)

  enddo

  ! Deallocate diagnostic remapping structures

  deallocate(grid_storage%diag_grids)

end subroutine diag_grid_storage_end

!< Allocate and initialize the masks for downsampled diagnostics in diag_cs

!! The downsampled masks in the axes would later "point" to these.

subroutine downsample_diag_masks_set(G, nz, diag_cs)

  type(ocean_grid_type), target, intent(in) :: G  !< The ocean grid type.

  integer,                       intent(in) :: nz !< The number of layers in the model's native grid.

  type(diag_ctrl),               pointer    :: diag_cs !< A pointer to a type with many variables

                                                       !! used for diagnostics

  ! Local variables

  integer :: k, dl

!print*,'original c extents ',G%isc,G%iec,G%jsc,G%jec

!print*,'original c extents ',G%iscb,G%iecb,G%jscb,G%jecb

!print*,'coarse   c extents ',G%HId2%isc,G%HId2%iec,G%HId2%jsc,G%HId2%jec

!print*,'original d extents ',G%isd,G%ied,G%jsd,G%jed

!print*,'coarse   d extents ',G%HId2%isd,G%HId2%ied,G%HId2%jsd,G%HId2%jed

! original c  extents           5          52           5          52

! original cB-nonsym extents    5          52           5          52

! original cB-sym    extents    4          52           4          52

! coarse   c extents            3          26           3          26

! original d extents            1          56           1          56

! original dB-nonsym extents    1          56           1          56

! original dB-sym extents       0          56           0          56

! coarse   d extents            1          28           1          28

  do dl=2,max_dsamp_lev

    ! 2d mask

    call downsample_mask(g%mask2dT, diag_cs%dsamp(dl)%mask2dT,  dl, g%isc, g%jsc, g%isd, g%jsd, &

            g%HId2%isc, g%HId2%iec, g%HId2%jsc, g%HId2%jec, g%HId2%isd, g%HId2%ied, g%HId2%jsd, g%HId2%jed)

    call downsample_mask(g%mask2dBu, diag_cs%dsamp(dl)%mask2dBu, dl,g%IscB, g%JscB, g%IsdB, g%JsdB, &

            g%HId2%IscB,g%HId2%IecB, g%HId2%JscB,g%HId2%JecB,g%HId2%IsdB,g%HId2%IedB,g%HId2%JsdB,g%HId2%JedB)

    call downsample_mask(g%mask2dCu, diag_cs%dsamp(dl)%mask2dCu, dl, g%IscB, g%jsc, g%IsdB, g%jsd, &

            g%HId2%IscB,g%HId2%IecB, g%HId2%jsc, g%HId2%jec,g%HId2%IsdB,g%HId2%IedB,g%HId2%jsd, g%HId2%jed)

    call downsample_mask(g%mask2dCv, diag_cs%dsamp(dl)%mask2dCv, dl,g %isc ,g%JscB, g%isd, g%JsdB, &

            g%HId2%isc ,g%HId2%iec, g%HId2%JscB,g%HId2%JecB,g%HId2%isd ,g%HId2%ied, g%HId2%JsdB,g%HId2%JedB)

    ! 3d native masks are needed by diag_manager but the native variables

    ! can only be masked 2d - for ocean points, all layers exists.

    allocate(diag_cs%dsamp(dl)%mask3dTL(g%HId2%isd:g%HId2%ied,g%HId2%jsd:g%HId2%jed,1:nz))

    allocate(diag_cs%dsamp(dl)%mask3dBL(g%HId2%IsdB:g%HId2%IedB,g%HId2%JsdB:g%HId2%JedB,1:nz))

    allocate(diag_cs%dsamp(dl)%mask3dCuL(g%HId2%IsdB:g%HId2%IedB,g%HId2%jsd:g%HId2%jed,1:nz))

    allocate(diag_cs%dsamp(dl)%mask3dCvL(g%HId2%isd:g%HId2%ied,g%HId2%JsdB:g%HId2%JedB,1:nz))

    do k=1,nz

      diag_cs%dsamp(dl)%mask3dTL(:,:,k) = diag_cs%dsamp(dl)%mask2dT(:,:)

      diag_cs%dsamp(dl)%mask3dBL(:,:,k) = diag_cs%dsamp(dl)%mask2dBu(:,:)

      diag_cs%dsamp(dl)%mask3dCuL(:,:,k) = diag_cs%dsamp(dl)%mask2dCu(:,:)

      diag_cs%dsamp(dl)%mask3dCvL(:,:,k) = diag_cs%dsamp(dl)%mask2dCv(:,:)

    enddo

    allocate(diag_cs%dsamp(dl)%mask3dTi(g%HId2%isd:g%HId2%ied,g%HId2%jsd:g%HId2%jed,1:nz+1))

    allocate(diag_cs%dsamp(dl)%mask3dBi(g%HId2%IsdB:g%HId2%IedB,g%HId2%JsdB:g%HId2%JedB,1:nz+1))

    allocate(diag_cs%dsamp(dl)%mask3dCui(g%HId2%IsdB:g%HId2%IedB,g%HId2%jsd:g%HId2%jed,1:nz+1))

    allocate(diag_cs%dsamp(dl)%mask3dCvi(g%HId2%isd:g%HId2%ied,g%HId2%JsdB:g%HId2%JedB,1:nz+1))

    do k=1,nz+1

      diag_cs%dsamp(dl)%mask3dTi(:,:,k) = diag_cs%dsamp(dl)%mask2dT(:,:)

      diag_cs%dsamp(dl)%mask3dBi(:,:,k) = diag_cs%dsamp(dl)%mask2dBu(:,:)

      diag_cs%dsamp(dl)%mask3dCui(:,:,k) = diag_cs%dsamp(dl)%mask2dCu(:,:)

      diag_cs%dsamp(dl)%mask3dCvi(:,:,k) = diag_cs%dsamp(dl)%mask2dCv(:,:)

    enddo

  enddo

end subroutine downsample_diag_masks_set

!> Get the diagnostics-compute indices (to be passed to send_data) based on the shape of

!! the diagnostic field (the same way they are deduced for non-downsampled fields)

subroutine downsample_diag_indices_get(fo1, fo2, dl, diag_cs, isv, iev, jsv, jev)

  integer,           intent(in)  :: fo1     !< The size of the diag field in x

  integer,           intent(in)  :: fo2     !< The size of the diag field in y

  integer,           intent(in)  :: dl      !< Integer downsample level

  type(diag_ctrl),   intent(in)  :: diag_CS !< Structure used to regulate diagnostic output

  integer,           intent(out) :: isv     !< i-start index for diagnostics

  integer,           intent(out) :: iev     !< i-end index for diagnostics

  integer,           intent(out) :: jsv     !< j-start index for diagnostics

  integer,           intent(out) :: jev     !< j-end index for diagnostics

  ! Local variables

  integer :: dszi, cszi, dszj, cszj, f1, f2

  character(len=500) :: mesg

  logical, save :: first_check = .true.

  !   The current implementation of the downsampled diagnostics assumes that the tracer-point

  ! computational domain on each processor can be evenly divided by dL in each direction, which

  ! avoids the need for halo updates or checks that the halo regions are up-to-date.  The following

  ! check that this assumption is true is only relevant if there are in fact downsampled diagnostics,

  ! which is why it occurs during the first call to this routine instead of during initialization.

  if (first_check) then

    if (mod(diag_cs%ie-diag_cs%is+1, dl) /= 0 .OR. mod(diag_cs%je-diag_cs%js+1, dl) /= 0) then

      write (mesg,*) "Non-commensurate downsampled domain is not supported. "//&

             "Please choose a layout such that NIGLOBAL/Layout_X and NJGLOBAL/Layout_Y are both divisible by dl=",dl,&

             " Current domain extents: ", diag_cs%is,diag_cs%ie, diag_cs%js,diag_cs%je

      call mom_error(fatal,"downsample_diag_indices_get: "//trim(mesg))

    endif

    first_check = .false.

  endif

  cszi = diag_cs%dsamp(dl)%iec-diag_cs%dsamp(dl)%isc +1 ; dszi = diag_cs%dsamp(dl)%ied-diag_cs%dsamp(dl)%isd +1

  cszj = diag_cs%dsamp(dl)%jec-diag_cs%dsamp(dl)%jsc +1 ; dszj = diag_cs%dsamp(dl)%jed-diag_cs%dsamp(dl)%jsd +1

  isv = diag_cs%dsamp(dl)%isc ; iev = diag_cs%dsamp(dl)%iec

  jsv = diag_cs%dsamp(dl)%jsc ; jev = diag_cs%dsamp(dl)%jec

  f1 = fo1/dl

  f2 = fo2/dl

  ! Correction for the symmetric case

  if (diag_cs%G%symmetric) then

    f1 = f1 + mod(fo1,dl)

    f2 = f2 + mod(fo2,dl)

  endif

  ! Find the range of indices in the downscaled computational domain.

  if ( f1 == dszi ) then

    isv = diag_cs%dsamp(dl)%isc ; iev = diag_cs%dsamp(dl)%iec   ! Field on Data domain, take compute domain indices

  elseif ( f1 == dszi + 1 ) then

    isv = diag_cs%dsamp(dl)%isc ; iev = diag_cs%dsamp(dl)%iec+1   ! Symmetric data domain

  elseif ( f1 == cszi) then

    isv = 1 ; iev = (diag_cs%dsamp(dl)%iec-diag_cs%dsamp(dl)%isc) +1  ! Computational domain

  elseif ( f1 == cszi + 1 ) then

    isv = 1 ; iev = (diag_cs%dsamp(dl)%iec-diag_cs%dsamp(dl)%isc) +2  ! Symmetric computational domain

  else

    write (mesg,*) " peculiar size ",f1," in i-direction\n"//&

          "does not match one of ", cszi, cszi+1, dszi, dszi+1

    call mom_error(fatal,"downsample_diag_indices_get: "//trim(mesg))

  endif

  if ( f2 == dszj ) then

    jsv = diag_cs%dsamp(dl)%jsc ; jev = diag_cs%dsamp(dl)%jec     ! Data domain

  elseif ( f2 == dszj + 1 ) then

    jsv = diag_cs%dsamp(dl)%jsc ; jev = diag_cs%dsamp(dl)%jec+1   ! Symmetric data domain

  elseif ( f2 == cszj) then

    jsv = 1 ; jev = (diag_cs%dsamp(dl)%jec-diag_cs%dsamp(dl)%jsc) +1  ! Computational domain

  elseif ( f2 == cszj + 1 ) then

    jsv = 1 ; jev = (diag_cs%dsamp(dl)%jec-diag_cs%dsamp(dl)%jsc) +2  ! Symmetric computational domain

  else

    write (mesg,*) " peculiar size ",f2," in j-direction\n"//&

          "does not match one of ", cszj, cszj+1, dszj, dszj+1

    call mom_error(fatal,"downsample_diag_indices_get: "//trim(mesg))

  endif

end subroutine downsample_diag_indices_get

!> This subroutine allocates and computes a downsampled array from an input array.

!! It also determines the diagnostic computational grid indices for the downsampled array.

!! 3d interface

subroutine downsample_diag_field_3d(locfield, locfield_dsamp, dl, diag_cs, diag, isv, iev, jsv, jev, mask)

  real, dimension(:,:,:), pointer :: locfield  !< Input array pointer in arbitrary units [A ~> a]

  real, dimension(:,:,:), allocatable, intent(inout) :: locfield_dsamp !< Output (downsampled) array [A ~> a]

  type(diag_ctrl),   intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  type(diag_type),   intent(in) :: diag    !< A structure describing the diagnostic to post

  integer, intent(in) :: dl                !< Level of down sampling

  integer, intent(inout) :: isv            !< i-start index for diagnostics

  integer, intent(inout) :: iev            !< i-end index for diagnostics

  integer, intent(inout) :: jsv            !< j-start index for diagnostics

  integer, intent(inout) :: jev            !< j-end index for diagnostics

  real,    optional,target, intent(in) :: mask(:,:,:) !< If present, use this real array as the data mask [nondim]

  ! Local variables

  real, dimension(:,:,:), pointer :: locmask ! A pointer to the mask [nondim]

  integer :: f1, f2, isv_o, jsv_o

  locmask => null()

  ! Get the correct indices corresponding to input field based on its shape.

  f1 = size(locfield, 1)

  f2 = size(locfield, 2)

  ! Save the extents of the original (fine) domain

  isv_o = isv ; jsv_o = jsv

  ! Get the shape of the downsampled field and overwrite isv, iev, jsv and jev with them

  call downsample_diag_indices_get(f1, f2, dl, diag_cs, isv, iev, jsv, jev)

  ! Set the pointer to the non-downsampled mask, which must be associated and initialized

  if (present(mask)) then

    locmask => mask

  elseif (associated(diag%axes%mask3d)) then

    locmask => diag%axes%mask3d

  else

    call mom_error(fatal, "downsample_diag_field_3d: Cannot downsample without a mask!!! ")

  endif

  call downsample_field(locfield, locfield_dsamp, dl, diag%xyz_method, locmask, diag_cs, diag, &

                        isv_o, jsv_o, isv, iev, jsv, jev)

end subroutine downsample_diag_field_3d

!> This subroutine allocates and computes a downsampled array from an input array.

!! It also determines the diagnostic computational grid indices for the downsampled array.

!! 2d interface

subroutine downsample_diag_field_2d(locfield, locfield_dsamp, dl, diag_cs, diag, isv, iev, jsv, jev, mask)

  real, dimension(:,:), pointer :: locfield !< Input array pointer in arbitrary units [A ~> a]

  real, dimension(:,:), allocatable, intent(inout) :: locfield_dsamp !< Output (downsampled) array [A ~> a]

  type(diag_ctrl),   intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  type(diag_type),   intent(in) :: diag    !< A structure describing the diagnostic to post

  integer, intent(in) :: dl                !< Level of down sampling

  integer, intent(inout) :: isv            !< i-start index for diagnostics

  integer, intent(inout) :: iev            !< i-end index for diagnostics

  integer, intent(inout) :: jsv            !< j-start index for diagnostics

  integer, intent(inout) :: jev            !< j-end index for diagnostics

  real,    optional,target, intent(in) :: mask(:,:) !< If present, use this real array as the data mask [nondim].

  ! Local variables

  real, dimension(:,:), pointer :: locmask ! A pointer to the mask [nondim]

  integer :: f1, f2, isv_o, jsv_o

  locmask => null()

  ! Get the correct indices corresponding to input field based on its shape.

  f1 = size(locfield,1)

  f2 = size(locfield,2)

  ! Save the extents of the original (fine) domain

  isv_o = isv ; jsv_o = jsv

  ! Get the shape of the downsampled field and overwrite isv, iev, jsv and jev with them

  call downsample_diag_indices_get(f1, f2, dl, diag_cs, isv, iev, jsv, jev)

  ! Set the non-downsampled mask, it must be associated and initialized

  if (present(mask)) then

    locmask => mask

  elseif (associated(diag%axes%mask2d)) then

    locmask => diag%axes%mask2d

  else

    call mom_error(fatal, "downsample_diag_field_2d: Cannot downsample without a mask!!! ")

  endif

  call downsample_field(locfield, locfield_dsamp, dl, diag%xyz_method, locmask, diag_cs,diag, &

                        isv_o, jsv_o, isv, iev, jsv, jev)

end subroutine downsample_diag_field_2d

!> \section downsampling The down sample algorithm

!!

!! The down sample method could be deduced (before send_data call)

!!  from the diag%x_cell_method, diag%y_cell_method and diag%v_cell_method

!!

!! This is the summary of the down sample algorithm for a diagnostic field f:

!!  \f[

!!     f(Id,Jd) = \sum_{i,j} f(Id+i,Jd+j) * weight(Id+i,Jd+j) / [ \sum_{i,j} weight(Id+i,Jd+j)]

!!  \f]

!!  Here, i and j run from 0 to dl-1 (dl being the down sample level).

!!  Id,Jd are the down sampled (coarse grid) indices run over the coarsened compute grid,

!!  if and jf are the original (fine grid) indices.

!!

!! \verbatim

!! Example   x_cell y_cell v_cell algorithm_id    implemented weight(if,jf)

!! ---------------------------------------------------------------------------------------

!! theta     mean   mean   mean   MMM =222        G%areaT(if,jf)*h(if,jf)

!! u         point  mean   mean   PMM =022        dyCu(if,jf)*h(if,jf)*delta(if,Id)

!! v         mean   point  mean   MPM =202        dxCv(if,jf)*h(if,jf)*delta(jf,Jd)

!! ?         point  sum    mean   PSM =012        h(if,jf)*delta(if,Id)

!! volcello  sum    sum    sum    SSS =111        1

!! T_dfxy_co sum    sum    point  SSP =110        1

!! umo       point  sum    sum    PSS =011        1*delta(if,Id)

!! vmo       sum    point  sum    SPS =101        1*delta(jf,Jd)

!! umo_2d    point  sum    point  PSP =010        1*delta(if,Id)

!! vmo_2d    sum    point  point  SPP =100        1*delta(jf,Jd)

!! ?         point  mean   point  PMP =020        dyCu(if,jf)*delta(if,Id)

!! ?         mean   point  point  MPP =200        dxCv(if,jf)*delta(jf,Jd)

!! w         mean   mean   point  MMP =220        G%areaT(if,jf)

!! h*theta   mean   mean   sum    MMS =221        G%areaT(if,jf)

!!

!! delta is the Kronecker delta

!! \endverbatim

!> This subroutine allocates and computes a down sampled 3d array given an input array

!! The down sample method is based on the "cell_methods" for the diagnostics as explained

!! in the above table

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)

  real, dimension(:,:,:), pointer :: field_in      !< Original field to be downsampled in arbitrary units [A ~> a]

  real, dimension(:,:,:), allocatable :: field_out !< Downsampled field in the same arbitrary units [A ~> a]

  integer, intent(in) :: dl                !< Level of down sampling

  integer,  intent(in) :: method           !< Sampling method

  real,  dimension(:,:,:), pointer :: mask !< Mask for field [nondim]

  type(diag_ctrl), intent(in) :: diag_CS   !< Structure used to regulate diagnostic output

  type(diag_type), intent(in) :: diag      !< A structure describing the diagnostic to post

  integer, intent(in) :: isv_o             !< Original i-start index

  integer, intent(in) :: jsv_o             !< Original j-start index

  integer, intent(in) :: isv_d             !< i-start index of down sampled data

  integer, intent(in) :: iev_d             !< i-end index of down sampled data

  integer, intent(in) :: jsv_d             !< j-start index of down sampled data

  integer, intent(in) :: jev_d             !< j-end index of down sampled data

  ! Local variables

  character(len=240) :: mesg

  integer :: i, j, ii, jj, i0, j0, f1, f2, f_in1, f_in2

  integer :: k, ks, ke

  real :: ave       ! The running sum of the average, in [A ~> a], [A L2 ~> a m2],

                    ! [A H L ~> a m2 or a kg m-1] or [A H L2 ~> a m3 or a kg]

  real :: weight    ! The nondimensional, area-, volume- or mass-based weight for an input

                    ! value [nondim], [L2 ~> m2], [H L ~> m2 or kg m-1] or [H L2 ~> m3 or kg]

  real :: total_weight ! The sum of weights contributing to a point [nondim], [L2 ~> m2],

                    ! [H L ~> m2 or kg m-1] or [H L2 ~> m3 or kg]

  real :: eps_vol   ! A negligibly small volume or mass [H L2 ~> m3 or kg]

  real :: eps_area  ! A negligibly small area [L2 ~> m2]

  real :: eps_face  ! A negligibly small face area [H L ~> m2 or kg m-1]

  ks = 1 ; ke = size(field_in,3)

  eps_face = 1.0e-20 * diag_cs%G%US%m_to_L * diag_cs%GV%m_to_H

  eps_area = 1.0e-20 * diag_cs%G%US%m_to_L**2

  eps_vol = 1.0e-20 * diag_cs%G%US%m_to_L**2 * diag_cs%GV%m_to_H

  ! Allocate the down sampled field on the down sampled data domain

!  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))

!  allocate(field_out(1:size(field_in,1)/dl,1:size(field_in,2)/dl,ks:ke))

  f_in1 = size(field_in,1)

  f_in2 = size(field_in,2)

  f1 = f_in1/dl

  f2 = f_in2/dl

  ! Correction for the symmetric case

  if (diag_cs%G%symmetric) then

    f1 = f1 + mod(f_in1,dl)

    f2 = f2 + mod(f_in2,dl)

  endif

  allocate(field_out(1:f1,1:f2,ks:ke))

  ! Fill the down sampled field on the down sampled diagnostics (almost always compute) domain

  !### The averaging used here is not rotationally invariant.

  !  Also, it would be better to use a max with eps_vol instead of adding it in the denominator.

  if (method == mmm) then

    do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      total_weight = 0.0

      do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

        weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj) * diag_cs%h(ii,jj,k)

        total_weight = total_weight + weight

        ave = ave+field_in(ii,jj,k) * weight

      enddo ; enddo

      field_out(i,j,k)  = ave / (total_weight + eps_vol)  ! Eps_vol avoids division by 0.

    enddo ; enddo ; enddo

  elseif (method == sss) then   ! e.g., volcello

    do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

        weight = mask(ii,jj,k)

        ave = ave+field_in(ii,jj,k)*weight

      enddo ; enddo

      field_out(i,j,k)  = ave ! This is a masked sum, and total_weight = 1.

    enddo ; enddo ; enddo

  elseif (method == mmp .or. method == mms) then   ! e.g., T_advection_xy

    do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      total_weight = 0.0

      do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

        weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj)

        total_weight = total_weight + weight

        ave = ave+field_in(ii,jj,k) * weight

      enddo ; enddo

      field_out(i,j,k)  = ave / (total_weight + eps_area)  ! Eps_area avoids division by 0.

    enddo ; enddo ; enddo

  elseif (method == pmm) then

    do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      total_weight = 0.0

      ii=i0

      do jj=j0,j0+dl-1

        weight = mask(ii,jj,k) * diag_cs%G%dyCu(ii,jj) * diag_cs%h(ii,jj,k)

        total_weight = total_weight + weight

        ave = ave+field_in(ii,jj,k) * weight

      enddo

      field_out(i,j,k)  = ave / (total_weight + eps_face)  ! Eps_face avoids division by 0.

    enddo ; enddo ; enddo

  elseif (method == pss) then    ! e.g. umo

    do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      ii=i0

      do jj=j0,j0+dl-1

        weight = mask(ii,jj,k)

        ave = ave+field_in(ii,jj,k)*weight

      enddo

      field_out(i,j,k)  = ave   ! This is a masked sum, and total_weight = 1.

    enddo ; enddo ; enddo

  elseif (method == sps) then   ! e.g. vmo

    do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      jj=j0

      do ii=i0,i0+dl-1

        weight = mask(ii,jj,k)

        ave = ave+field_in(ii,jj,k)*weight

      enddo

      field_out(i,j,k)  = ave  ! This is a masked sum, and total_weight = 1.

    enddo ; enddo ; enddo

  elseif (method == mpm) then

    do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      total_weight = 0.0

      jj=j0

      do ii=i0,i0+dl-1

        weight = mask(ii,jj,k) * diag_cs%G%dxCv(ii,jj) * diag_cs%h(ii,jj,k)

        total_weight = total_weight + weight

        ave = ave+field_in(ii,jj,k) * weight

      enddo

      field_out(i,j,k)  = ave / (total_weight + eps_face)  ! Eps_face avoids division by 0.

    enddo ; enddo ; enddo

  elseif (method == msk) then ! The input field is a mask, so subsample it instead of averaging.

    field_out(:,:,:) = 0.0

    do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

        ave = ave+field_in(ii,jj,k)

      enddo ; enddo

      if (ave > 0.0) field_out(i,j,k)=1.0

    enddo ; enddo ; enddo

  else

    write (mesg,*) " unknown sampling method: ",method

    call mom_error(fatal, "downsample_field_3d: "//trim(mesg)//" "//trim(diag%debug_str))

  endif

end subroutine downsample_field_3d

!> This subroutine allocates and computes a down sampled 2d array given an input array

!! The down sample method is based on the "cell_methods" for the diagnostics as explained

!! in the above table

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)

  real, dimension(:,:), pointer :: field_in      !< Original field to be downsampled in arbitrary units [A ~> a]

  real, dimension(:,:), allocatable :: field_out !< Downsampled field in the same arbitrary units [A ~> a]

  integer, intent(in) :: dl                !< Level of down sampling

  integer,  intent(in) :: method           !< Sampling method

  real, dimension(:,:), pointer :: mask    !< Mask for field [nondim]

  type(diag_ctrl),   intent(in) :: diag_CS !< Structure used to regulate diagnostic output

  type(diag_type),   intent(in) :: diag    !< A structure describing the diagnostic to post

  integer, intent(in) :: isv_o             !< Original i-start index

  integer, intent(in) :: jsv_o             !< Original j-start index

  integer, intent(in) :: isv_d             !< i-start index of down sampled data

  integer, intent(in) :: iev_d             !< i-end index of down sampled data

  integer, intent(in) :: jsv_d             !< j-start index of down sampled data

  integer, intent(in) :: jev_d             !< j-end index of down sampled data

  ! Local variables

  character(len=240) :: mesg

  integer :: i, j, ii, jj, i0, j0, f1, f2, f_in1, f_in2

  real :: ave       ! The running sum of the average, in [A ~> a] or [A L2 ~> a m2]

  real :: weight    ! The nondimensional or area-weighted weight for an input value [nondim] or [L2 ~> m2]

  real :: total_weight ! The sum of weights contributing to a point [nondim] or [L2 ~> m2]

  real :: eps_area  ! A negligibly small area [L2 ~> m2]

  real :: eps_len   ! A negligibly small horizontal length [L ~> m]

  eps_len = 1.0e-20 * diag_cs%G%US%m_to_L

  eps_area = 1.0e-20 * diag_cs%G%US%m_to_L**2

  ! Allocate the down sampled field on the down sampled data domain

!  allocate(field_out(diag_cs%dsamp(dl)%isd:diag_cs%dsamp(dl)%ied,diag_cs%dsamp(dl)%jsd:diag_cs%dsamp(dl)%jed))

!  allocate(field_out(1:size(field_in,1)/dl,1:size(field_in,2)/dl))

  ! Fill the down sampled field on the down sampled diagnostics (almost always compute) domain

  f_in1 = size(field_in,1)

  f_in2 = size(field_in,2)

  f1 = f_in1/dl

  f2 = f_in2/dl

  ! Correction for the symmetric case

  if (diag_cs%G%symmetric) then

    f1 = f1 + mod(f_in1,dl)

    f2 = f2 + mod(f_in2,dl)

  endif

  allocate(field_out(1:f1,1:f2))

  if (method == mmp) then

    do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      total_weight = 0.0

      do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

        weight = mask(ii,jj)*diag_cs%G%areaT(ii,jj)

        total_weight = total_weight + weight

        ave = ave+field_in(ii,jj) * weight

      enddo ; enddo

      field_out(i,j) = ave / (total_weight + eps_area)  ! Eps_area avoids division by 0.

    enddo ; enddo

  elseif (method == ssp) then    ! e.g., T_dfxy_cont_tendency_2d

    do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

        weight = mask(ii,jj)

        ave = ave+field_in(ii,jj)*weight

      enddo ; enddo

      field_out(i,j) = ave  ! This is a masked sum, and total_weight = 1.

    enddo ; enddo

  elseif (method == psp) then   ! e.g., umo_2d

    do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      ii=i0

      do jj=j0,j0+dl-1

        weight = mask(ii,jj)

        ave = ave+field_in(ii,jj)*weight

      enddo

      field_out(i,j) = ave  ! This is a masked sum, and total_weight = 1.

    enddo ; enddo

  elseif (method == spp) then   ! e.g., vmo_2d

    do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      jj=j0

      do ii=i0,i0+dl-1

        weight = mask(ii,jj)

        ave = ave+field_in(ii,jj)*weight

      enddo

      field_out(i,j) = ave  ! This is a masked sum, and total_weight = 1.

    enddo ; enddo

  elseif (method == pmp) then

    do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      total_weight = 0.0

      ii=i0

      do jj=j0,j0+dl-1

        weight = mask(ii,jj) * diag_cs%G%dyCu(ii,jj)!*diag_cs%h(ii,jj,1) !Niki?

        total_weight = total_weight + weight

        ave = ave+field_in(ii,jj) * weight

      enddo

      field_out(i,j) = ave / (total_weight + eps_len)  ! Eps_len avoids division by 0.

    enddo ; enddo

  elseif (method == mpp) then

    do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      total_weight = 0.0

      jj=j0

      do ii=i0,i0+dl-1

        weight = mask(ii,jj)* diag_cs%G%dxCv(ii,jj)!*diag_cs%h(ii,jj,1) !Niki?

        total_weight = total_weight +weight

        ave = ave+field_in(ii,jj)*weight

      enddo

      field_out(i,j) = ave / (total_weight + eps_len)  ! Eps_len avoids division by 0.

    enddo ; enddo

  elseif (method == msk) then ! The input field is a mask, so subsample it instead of averaging.

    field_out(:,:) = 0.0

    do j=jsv_d,jev_d ; do i=isv_d,iev_d

      i0 = isv_o+dl*(i-isv_d)

      j0 = jsv_o+dl*(j-jsv_d)

      ave = 0.0

      do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

        ave = ave+field_in(ii,jj)

      enddo ; enddo

      if (ave > 0.0) field_out(i,j)=1.0

    enddo ; enddo

  else

    write (mesg,*) " unknown sampling method: ",method

    call mom_error(fatal, "downsample_field_2d: "//trim(mesg)//" "//trim(diag%debug_str))

  endif

end subroutine downsample_field_2d

!> Allocate and compute the 2d down sampled mask

!! The masks are down sampled based on a minority rule, i.e., a coarse cell is open (1)

!! if at least one of the sub-cells are open, otherwise it's closed (0)

subroutine downsample_mask_2d(field_in, field_out, dl, isc_o, jsc_o, isd_o, jsd_o, &

                              isc_d, iec_d, jsc_d, jec_d, isd_d, ied_d, jsd_d, jed_d)

  integer, intent(in) :: isd_o !< Original data domain i-start index

  integer, intent(in) :: jsd_o !< Original data domain j-start index

  real, dimension(isd_o:,jsd_o:), intent(in) :: field_in !< Original field to be down sampled in arbitrary units [A]

  real, dimension(:,:), pointer :: field_out   !< Down sampled field mask [nondim]

  integer, intent(in) :: dl    !< Level of down sampling

  integer, intent(in) :: isc_o !< Original i-start index

  integer, intent(in) :: jsc_o !< Original j-start index

  integer, intent(in) :: isc_d !< Computational i-start index of down sampled data

  integer, intent(in) :: iec_d !< Computational i-end index of down sampled data

  integer, intent(in) :: jsc_d !< Computational j-start index of down sampled data

  integer, intent(in) :: jec_d !< Computational j-end index of down sampled data

  integer, intent(in) :: isd_d !< Data domain i-start index of down sampled data

  integer, intent(in) :: ied_d !< Data domain i-end index of down sampled data

  integer, intent(in) :: jsd_d !< Data domain j-start index of down sampled data

  integer, intent(in) :: jed_d !< Data domain j-end index of down sampled data

  ! Local variables

  integer :: i, j, ii, jj, i0, j0

  real    :: tot_non_zero  ! The sum of values in the down-scaled cell [A]

  ! down sampled mask = 0 unless the mask value of one of the down sampling cells is 1

  allocate(field_out(isd_d:ied_d,jsd_d:jed_d))

  field_out(:,:) = 0.0

  do j=jsc_d,jec_d ; do i=isc_d,iec_d

    i0 = isc_o+dl*(i-isc_d)

    j0 = jsc_o+dl*(j-jsc_d)

    tot_non_zero = 0.0

    do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

      tot_non_zero = tot_non_zero + field_in(ii,jj)

    enddo ; enddo

    if (tot_non_zero > 0.0) field_out(i,j)=1.0

  enddo ; enddo

end subroutine downsample_mask_2d

!> Allocate and compute the 3d down sampled mask

!! The masks are down sampled based on a minority rule, i.e., a coarse cell is open (1)

!! if at least one of the sub-cells are open, otherwise it's closed (0)

subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isd_o, jsd_o, &

                              isc_d, iec_d, jsc_d, jec_d, isd_d, ied_d, jsd_d, jed_d)

  integer, intent(in) :: isd_o !< Original data domain i-start index

  integer, intent(in) :: jsd_o !< Original data domain j-start index

  real, dimension(isd_o:,jsd_o:,:), intent(in) :: field_in !< Original field to be down sampled in arbitrary units [A]

  real, dimension(:,:,:), pointer :: field_out   !< down sampled field mask [nondim]

  integer, intent(in) :: dl    !< Level of down sampling

  integer, intent(in) :: isc_o !< Original i-start index

  integer, intent(in) :: jsc_o !< Original j-start index

  integer, intent(in) :: isc_d !< Computational i-start index of down sampled data

  integer, intent(in) :: iec_d !< Computational i-end index of down sampled data

  integer, intent(in) :: jsc_d !< Computational j-start index of down sampled data

  integer, intent(in) :: jec_d !< Computational j-end index of down sampled data

  integer, intent(in) :: isd_d !< Computational i-start index of down sampled data

  integer, intent(in) :: ied_d !< Computational i-end index of down sampled data

  integer, intent(in) :: jsd_d !< Computational j-start index of down sampled data

  integer, intent(in) :: jed_d !< Computational j-end index of down sampled data

  ! Local variables

  integer :: i, j, ii, jj, i0, j0, k, ks, ke

  real    :: tot_non_zero  ! The sum of values in the down-scaled cell [A]

  ! down sampled mask = 0 unless the mask value of one of the down sampling cells is 1

  ks = lbound(field_in,3) ; ke = ubound(field_in,3)

  allocate(field_out(isd_d:ied_d,jsd_d:jed_d,ks:ke))

  field_out(:,:,:) = 0.0

  do k=ks,ke ; do j=jsc_d,jec_d ; do i=isc_d,iec_d

    i0 = isc_o+dl*(i-isc_d)

    j0 = jsc_o+dl*(j-jsc_d)

    tot_non_zero = 0.0

    do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1

      tot_non_zero = tot_non_zero + field_in(ii,jj,k)

    enddo ; enddo

    if (tot_non_zero > 0.0) field_out(i,j,k)=1.0

  enddo ; enddo ; enddo

end subroutine downsample_mask_3d

!> Fakes a register of a diagnostic to find out if an obsolete

!! parameter appears in the diag_table.

logical function found_in_diagtable(diag, varName)

  type(diag_ctrl),            intent(in) :: diag       !< A structure used to control diagnostics.

  character(len=*),           intent(in) :: varname    !< The obsolete diagnostic name

  ! Local

  integer :: handle ! Integer handle returned from diag_manager

  ! We use register_static_field_fms() instead of register_static_field() so

  ! that the diagnostic does not appear in the available diagnostics list.

  handle = register_static_field_infra('ocean_model', varname, diag%axesT1%handles)

  found_in_diagtable = (handle>0)

end function found_in_diagtable

end module mom_diag_mediator