diag_axis.F90

Go to the documentation of this file.

!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************

MODULE diag_axis_mod
#include <fms_platform.h>
  ! <CONTACT EMAIL="seth.underwood@noaa.gov">
  !   Seth Underwood
  ! </CONTACT>

  ! <OVERVIEW> <TT>diag_axis_mod</TT> is an integral part
  !   of diag_manager_mod. It helps to create axis IDs
  !   that are used in register_diag_field.
  ! </OVERVIEW>

  ! <DESCRIPTION> Users first create axis ID by calling
  !   diag_axis_init, then use this axis ID in
  !   register_diag_field.
  ! </DESCRIPTION>

  USE mpp_domains_mod, ONLY: domainug, domain1d, domain2d, mpp_get_compute_domain,&
       & mpp_get_domain_components, null_domain1d, null_domain2d, null_domainug,&
       & OPERATOR(.NE.), mpp_get_global_domain, mpp_get_domain_name
  USE fms_mod, ONLY: error_mesg, write_version_number, lowercase, uppercase,&
       & fms_error_handler, fatal, note
  USE diag_data_mod, ONLY: diag_axis_type, max_subaxes, max_axes,&
       & max_num_axis_sets, max_axis_attributes, debug_diag_manager,&
       & first_send_data_call, diag_atttype
#ifdef use_netCDF
  USE netcdf, ONLY: nf90_int, nf90_float, nf90_char
#endif

  IMPLICIT NONE

  PRIVATE
  PUBLIC  diag_axis_init, get_diag_axis, get_domain1d, get_domain2d,&
       & get_axis_length, get_axis_global_length, diag_subaxes_init,&
       & get_diag_axis_cart, get_diag_axis_data, max_axes, get_axis_aux,&
       & get_tile_count, get_axes_shift, get_diag_axis_name,&
       & get_axis_num, get_diag_axis_domain_name, diag_axis_add_attribute,&
       & get_domainug, axis_compatible_check, axis_is_compressed, &
       & get_compressed_axes_ids, get_axis_reqfld

  ! Module variables
  ! Parameters
  ! Include variable "version" to be written to log file.
#include<file_version.h>

!----------
!ug support
  integer(INT_KIND),parameter,public :: diag_axis_nodomain = 0
  integer(INT_KIND),parameter,public :: diag_axis_2ddomain = 1
  integer(INT_KIND),parameter,public :: diag_axis_ugdomain = 2
!----------

  ! counter of number of axes defined
  INTEGER, DIMENSION(:), ALLOCATABLE :: num_subaxes
  INTEGER :: num_def_axes = 0

  ! storage for axis set names
  CHARACTER(len=128), DIMENSION(:), ALLOCATABLE, SAVE :: axis_sets
  INTEGER :: num_axis_sets = 0

  ! ---- global storage for all defined axes ----
  TYPE(diag_axis_type), ALLOCATABLE, SAVE :: axes(:)
  LOGICAL :: module_is_initialized = .false.

  ! <INTERFACE NAME="diag_axis_add_attribute">
  !   <OVERVIEW>
  !     Add a attribute to the diag axis
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE diag_axis_add_attribute(diag_axis_id, att_name, att_value)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Add an arbitrary attribute and value to the diagnostic axis.  Any number
  !     of attributes can be added to a given axis.  All attribute addition must
  !     be done before first <TT>send_data</TT> call.
  !
  !     If a real or integer attribute is already defined, a FATAL error will be called.
  !     If a character attribute is already defined, then it will be prepended to the
  !     existing attribute value.
  !   </DESCRIPTION>
  !   <IN NAME="diag_axis_id" TYPE="INTEGER" />
  !   <IN NAME="att_name" TYPE="CHARACTER(len=*)" />
  !   <IN NAME="att_value" TYPE="REAL|INTEGER|CHARACTER(len=*)" />
  INTERFACE diag_axis_add_attribute
     MODULE PROCEDURE diag_axis_add_attribute_scalar_r
     MODULE PROCEDURE diag_axis_add_attribute_scalar_i
     MODULE PROCEDURE diag_axis_add_attribute_scalar_c
     MODULE PROCEDURE diag_axis_add_attribute_r1d
     MODULE PROCEDURE diag_axis_add_attribute_i1d
  END INTERFACE diag_axis_add_attribute
  ! </INTERFACE>

CONTAINS

  ! <FUNCTION NAME="diag_axis_init">
  !   <OVERVIEW>
  !     Initialize the axis, and return the axis ID.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     INTEGER FUNCTION diag_axis_init(name, data, units, cart_name, long_name,
  !           direction, set_name, edges, Domain, Domain2, aux, req, tile_count)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     <TT>diag_axis_init</TT> initializes an axis and returns the axis ID that
  !     is to be used with <TT>register_diag_field</TT>.  This function also
  !     increments the axis counter and fills in the axes
  !   </DESCRIPTION>
  !   <IN NAME="name" TYPE="CHARACTER(len=*)">Short name for axis</IN>
  !   <IN NAME="data" TYPE="REAL, DIMENSION(:)">Array of coordinate values</IN>
  !   <IN NAME="units" TYPE="CHARACTER(len=*)">Units for the axis</IN>
  !   <IN NAME="cart_name" TYPE="CHARACTER(len=*)">
  !     Cartesian axis ("X", "Y", "Z", "T")
  !   </IN>
  !   <IN NAME="direction" TYPE="INTEGER, OPTIONAL" DEFAULT="0">
  !     Indicates the direction of the axis:
  !     <UL>
  !       <LI>Up if +1</LI>
  !       <LI>Down if -1</LI>
  !       <LI>Neither up or down if 0</LI>
  !     </UL>
  !   </IN>
  !   <IN NAME="long_name" TYPE="CHARACTER(len=*), OPTIONAL" DEFAULT="name">
  !     Long name for the axis.
  !   </IN>
  !   <IN NAME="edges" TYPE="INTEGER, OPTIONAL">
  !     Axis ID for the previously defined "edges axis"
  !   </IN>
  !   <IN NAME="Domain" TYPE="TYPE(domain1d), OPTIONAL" />
  !   <IN NAME="Domain2" TYPE="TYPE(domain2d), OPTIONAL" />
  !   <IN NAME="aux" TYPE="CHARACTER(len=*), OPTIONAL">
  !     Auxiliary name, can only be <TT>geolon_t</TT> or <TT>geolat_t</TT>
  !   </IN>
  !   <IN NAME="req" TYPE="CHARACTER(len=*), OPTIONAL">
  !     Required field names.
  !   </IN>
  !   <IN NAME="tile_count" TYPE="INTEGER, OPTIONAL" />
  INTEGER FUNCTION diag_axis_init(name, DATA, units, cart_name, long_name, direction,&
       & set_name, edges, Domain, Domain2, DomainU, aux, req, tile_count)
    CHARACTER(len=*), INTENT(in) :: name
    REAL, DIMENSION(:), INTENT(in) :: data
    CHARACTER(len=*), INTENT(in) :: units
    CHARACTER(len=*), INTENT(in) :: cart_name
    CHARACTER(len=*), INTENT(in), OPTIONAL :: long_name, set_name
    INTEGER, INTENT(in), OPTIONAL :: direction, edges
    TYPE(domain1d), INTENT(in), OPTIONAL :: domain
    TYPE(domain2d), INTENT(in), OPTIONAL :: domain2
    TYPE(domainug), INTENT(in), OPTIONAL :: domainu
    CHARACTER(len=*), INTENT(in), OPTIONAL :: aux, req
    INTEGER, INTENT(in), OPTIONAL :: tile_count

    TYPE(domain1d) :: domain_x, domain_y
    INTEGER :: ierr, axlen
    INTEGER :: i, set, tile
    INTEGER :: isc, iec, isg, ieg
    CHARACTER(len=128) :: emsg

    IF ( .NOT.module_is_initialized ) THEN
       CALL write_version_number("DIAG_AXIS_MOD", version)
    ENDIF

    IF ( PRESENT(tile_count)) THEN
       tile = tile_count
    ELSE
       tile = 1
    END IF

    ! Allocate the axes
    IF (.NOT. ALLOCATED(axis_sets)) ALLOCATE(axis_sets(max_num_axis_sets))
    IF (.NOT. ALLOCATED(axes)) ALLOCATE(axes(max_axes))
    IF (.NOT. ALLOCATED(num_subaxes)) THEN
       ALLOCATE(num_subaxes(max_axes))
       num_subaxes = 0
    END IF

    !---- is there an axis set? ----
    IF ( PRESENT(set_name) ) THEN
       set = get_axis_set_num(set_name)
       !---- add new set name ----
       IF (set == 0) THEN
          num_axis_sets = num_axis_sets + 1
          IF ( num_axis_sets > max_num_axis_sets ) THEN
             WRITE (emsg, fmt='("num_axis_sets (",I2,") exceeds max_num_axis_sets (",I2,"). ")')&
                  & num_axis_sets, max_num_axis_sets
             ! <ERROR STATUS="FATAL">
             !   num_axis_sets (<num_axis_sets>) exceeds max_num_axis_sets(<num_axis_sets>).
             !   Increase max_num_axis_sets via diag_manager_nml.
             ! </ERROR>
             CALL error_mesg('diag_axis_mod::diag_axis_init',&
                  & trim(emsg)//'  Increase max_num_axis_sets via diag_manager_nml.', fatal)
          END IF
          set = num_axis_sets
          axis_sets(set) = set_name
       END IF
    ELSE
       set = 0
    END IF

    !---- see if axis already exists --
    ! if this is time axis, return the ID of a previously defined
    ! if this is spatial axis, FATAL error
    DO i = 1, num_def_axes
       IF ( trim(name) == axes(i)%name ) THEN
          IF ( trim(name) == 'Stations' .OR. trim(name) == 'Levels') THEN
             diag_axis_init = i
             RETURN
          ELSE IF ( set == axes(i)%set ) THEN
             IF ( trim(lowercase(name)) == 'time' .OR.&
                  & trim(lowercase(cart_name)) == 't' .OR.&
                  & trim(lowercase(name)) == 'nv' .OR.&
                  & trim(lowercase(cart_name)) == 'n' ) THEN
                diag_axis_init = i
                RETURN
             ELSE IF ( (lowercase(cart_name) /= 'x' .AND. lowercase(cart_name) /= 'y')&
                  & .OR. tile /= axes(i)%tile_count) THEN
                ! <ERROR STATUS="FATAL">axis_name <NAME> and axis_set already exist.</ERROR>
                CALL error_mesg('diag_axis_mod::diag_axis_init',&
                     & 'axis_name '//trim(name)//' and axis_set already exist.', fatal)
             END IF
          END IF
       END IF
    END DO

    !---- register axis ----
    num_def_axes = num_def_axes + 1
    ! <ERROR STATUS="FATAL">max_axes exceeded, increase it via diag_manager_nml</ERROR>
    IF (num_def_axes > max_axes) CALL error_mesg ('diag_axis_mod::diag_axis_init',&
         & 'max_axes exceeded, increase via diag_manager_nml', fatal)
    diag_axis_init = num_def_axes

    !---- check and then save cart_name name ----
    IF ( trim(uppercase(cart_name)) == 'X' .OR.&
         & trim(uppercase(cart_name)) == 'Y' .OR.&
         & trim(uppercase(cart_name)) == 'Z' .OR.&
         & trim(uppercase(cart_name)) == 'T' .OR.&
         & trim(uppercase(cart_name)) == 'U' .OR.&
         & trim(uppercase(cart_name)) == 'N' ) THEN
       axes(diag_axis_init)%cart_name = trim(uppercase(cart_name))
    ELSE
       ! <ERROR STATUS="FATAL">Invalid cart_name name.</ERROR>
       CALL error_mesg('diag_axis_mod::diag_axis_init', 'Invalid cart_name name.', fatal)
    END IF

    !---- allocate storage for coordinate values of axis ----
    IF ( axes(diag_axis_init)%cart_name == 'T' ) THEN
       axlen = 0
    ELSE
       axlen = SIZE(DATA(:))
    END IF
    ALLOCATE ( axes(diag_axis_init)%data(1:axlen) )

    ! Initialize Axes(diag_axis_init)
    axes(diag_axis_init)%name   = trim(name)
    axes(diag_axis_init)%data   = DATA(1:axlen)
    axes(diag_axis_init)%units  = units
    axes(diag_axis_init)%length = axlen
    axes(diag_axis_init)%set    = set
    ! start and end are used in subaxes information only
    axes(diag_axis_init)%start = -1
    axes(diag_axis_init)%end = -1
    axes(diag_axis_init)%subaxis_name = ""
    axes(diag_axis_init)%shift = 0
    axes(diag_axis_init)%num_attributes = 0

    IF ( PRESENT(long_name) ) THEN
       axes(diag_axis_init)%long_name = long_name
    ELSE
       axes(diag_axis_init)%long_name = name
    END IF

    IF ( PRESENT(aux) ) THEN
       axes(diag_axis_init)%aux = trim(aux)
    ELSE
       axes(diag_axis_init)%aux = 'none'
    END IF

    IF ( PRESENT(req) ) THEN
       axes(diag_axis_init)%req = trim(req)
    ELSE
       axes(diag_axis_init)%req = 'none'
    END IF


    !---- axis direction (-1, 0, or +1) ----
    IF ( PRESENT(direction) )THEN
       IF ( abs(direction) /= 1 .AND. direction /= 0 )&
            ! <ERROR STATUS="FATAL">direction must be 0, +1, or -1</ERROR>
            & CALL error_mesg('diag_axis_mod::diag_axis_init', 'direction must be 0, +1 or -1', fatal)
       axes(diag_axis_init)%direction = direction
    ELSE
       axes(diag_axis_init)%direction = 0
    END IF

    !---- Handle the DomainU check
    IF (present(domainu) .AND. (PRESENT(domain2) .OR. PRESENT(domain)) ) THEN
       ! <ERROR STATUS="FATAL">Presence of DomainU and another Domain at the same time is prohibited</ERROR>
       CALL error_mesg('diag_axis_mod::diag_axis_init',&
            & 'Presence of DomainU and another Domain at the same time is prohibited', fatal)
    !---- domain2d type ----
    ELSE IF ( PRESENT(domain2) .AND. PRESENT(domain)) THEN
       ! <ERROR STATUS="FATAL">Presence of both Domain and Domain2 at the same time is prohibited</ERROR>
       CALL error_mesg('diag_axis_mod::diag_axis_init',&
            & 'Presence of both Domain and Domain2 at the same time is prohibited', fatal)
    ELSE IF ( PRESENT(domain2) .OR. PRESENT(domain)) THEN
       IF ( axes(diag_axis_init)%cart_name /= 'X' .AND. axes(diag_axis_init)%cart_name /= 'Y') THEN
          ! <ERROR STATUS="FATAL">Domain must not be present for an axis which is not in the X or Y direction.</ERROR>
          CALL error_mesg('diag_axis_mod::diag_axis_init',&
               & 'A Structured Domain must not be present for an axis which is not in the X or Y direction', fatal)
       END IF
    ELSE IF (present(domainu) .AND. axes(diag_axis_init)%cart_name /= 'U') THEN
          CALL error_mesg('diag_axis_mod::diag_axis_init',&
               & 'In the unstructured domain, the axis cart_name must be U', fatal)
    END IF

    axes(diag_axis_init)%tile_count = tile

    IF ( PRESENT(domain2) ) THEN
       axes(diag_axis_init)%Domain2 = domain2
       CALL mpp_get_domain_components(domain2, domain_x, domain_y, tile_count=tile_count)
       IF ( axes(diag_axis_init)%cart_name == 'X' ) axes(diag_axis_init)%Domain = domain_x
       IF ( axes(diag_axis_init)%cart_name == 'Y' ) axes(diag_axis_init)%Domain = domain_y
       axes(diag_axis_init)%DomainUG = null_domainug
    ELSE IF ( PRESENT(domain)) THEN
       !---- domain1d type ----
       axes(diag_axis_init)%Domain2 = null_domain2d ! needed since not 2-D domain
       axes(diag_axis_init)%Domain = domain
       axes(diag_axis_init)%DomainUG = null_domainug
    ELSE IF (present(domainu)) THEN
       axes(diag_axis_init)%Domain2 = null_domain2d
       axes(diag_axis_init)%Domain = null_domain1d
       axes(diag_axis_init)%DomainUG = domainu
    ELSE
       axes(diag_axis_init)%Domain2 = null_domain2d
       axes(diag_axis_init)%Domain = null_domain1d
       axes(diag_axis_init)%DomainUG = null_domainug
    END IF

    !--- set up the shift value for x-y axis
    IF ( axes(diag_axis_init)%Domain .NE. null_domain1d ) THEN
       CALL mpp_get_compute_domain(axes(diag_axis_init)%Domain, isc, iec)
       CALL mpp_get_global_domain(axes(diag_axis_init)%Domain, isg, ieg)
       IF ( axes(diag_axis_init)%length == ieg - isg + 2 ) THEN
          axes(diag_axis_init)%shift = 1
       END IF
    END IF

    !---- have axis edges been defined ? ----
    axes(diag_axis_init)%edges = 0
    IF (PRESENT(edges) ) THEN
       IF ( edges > 0 .AND. edges < num_def_axes ) THEN
          ierr=0
          IF ( axes(edges)%cart_name /= axes(diag_axis_init)%cart_name) ierr=1
          IF ( axes(edges)%length    /= axes(diag_axis_init)%length+1 ) ierr=ierr+2
          IF ( axes(edges)%set       /= axes(diag_axis_init)%set      ) ierr=ierr+4
          IF ( ierr > 0 )   THEN
             ! <ERROR STATUS="FATAL">Edges axis does not match axis (code <CODE>).</ERROR>
             WRITE (emsg,'("Edges axis does not match axis (code ",I1,").")') ierr
             CALL error_mesg('diag_axis_mod::diag_axis_init', emsg, fatal)
          END IF
          axes(diag_axis_init)%edges = edges
       ELSE
          ! <ERROR STATUS="FATAL">Edges axis is not defined.</ERROR>
          CALL error_mesg('diag_axis_mod::diag_axis_init', 'Edges axis is not defined', fatal)
       END IF
    END IF

    ! Module is now initialized
    module_is_initialized = .true.

  END FUNCTION diag_axis_init
  ! </FUNCTION>

  ! <FUNCTION NAME="diag_subaxes_init">
  !   <OVERVIEW>
  !     Create a subaxis on a parent axis.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     INTEGER FUNCTION diag_subaxes_init(axis, subdata, start_indx, end_indx,
  !           domain_1d, domain_2d)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Given the ID of a parent axis, create a subaxis and fill it with data,
  !     and return the ID of the corresponding subaxis.
  !
  !     The subaxis is defined on the parent axis from <TT>start_indx</TT>
  !     to <TT>end_indx</TT>.
  !   </DESCRIPTION>
  !   <IN NAME="axis" TYPE="INTEGER">ID of the parent axis</IN>
  !   <IN NAME="subdata" TYPE="REAL, DIMENSION(:)">Data of the subaxis</IN>
  !   <IN NAME="start_indx" TYPE="INTEGER">Start index of the subaxis</IN>
  !   <IN NAME="end_indx" TYPE="INTEGER">End index of the subaxis</IN>
  !   <IN NAME="domain_1d" TYPE="TYPE(domain1d), OPTIONAL" />
  !   <IN NAME="domain_2d" TYPE="TYPE(domain2d), OPTIONAL" />
  INTEGER FUNCTION diag_subaxes_init(axis, subdata, start_indx, end_indx, domain_2d)
    INTEGER, INTENT(in) :: axis
    REAL, DIMENSION(:), INTENT(in) :: subdata
    INTEGER, INTENT(in) :: start_indx
    INTEGER, INTENT(in) :: end_indx
    TYPE(domain2d), INTENT(in), OPTIONAL  :: domain_2d

    INTEGER :: i, nsub_axis, direction
    INTEGER :: xbegin, xend, ybegin, yend
    INTEGER :: ad_xbegin, ad_xend, ad_ybegin, ad_yend
    CHARACTER(len=128) :: name, nsub_name
    CHARACTER(len=128) :: units
    CHARACTER(len=128) :: cart_name
    CHARACTER(len=128) :: long_name
    CHARACTER(len=128) :: emsg
    LOGICAL :: subaxis_set, hasdomain

    ! there may be more than 1 subaxis on a parent axis, check for redundancy
    nsub_axis = 0
    subaxis_set = .false.

    IF ( PRESENT(domain_2d) ) THEN
       hasdomain = .true.
       CALL mpp_get_compute_domain(domain_2d, xbegin, xend, ybegin, yend)
    ELSE
       hasdomain = .false.
    END IF
    sa_search: DO i = 1, num_subaxes(axis)
       IF ( start_indx == axes(axis)%start(i) .AND. end_indx == axes(axis)%end(i) ) THEN
          IF ( hasdomain ) THEN
             CALL mpp_get_compute_domain(axes(axis)%subaxis_domain2(i), ad_xbegin, ad_xend, ad_ybegin, ad_yend)
             IF ( .NOT.((xbegin == ad_xbegin .AND. xend == ad_xend) .AND.&
                  & (ybegin == ad_ybegin .AND. yend == ad_yend)) ) THEN
                cycle sa_search
             END IF
          END IF
          nsub_axis = i
          subaxis_set = .true.    !subaxis already exists
          name = trim(axes(axis)%subaxis_name(nsub_axis))
          EXIT sa_search
       END IF
    END DO sa_search

    IF ( nsub_axis == 0 ) THEN  ! create new subaxis
       num_subaxes(axis) = num_subaxes(axis) + 1
       IF (num_subaxes(axis) > max_subaxes) THEN
          ! <ERROR STATUS="FATAL">max_subaxes (value <VALUE>) is too small.  Consider increasing max_subaxes.</ERROR>
          WRITE (emsg,'("max_subaxes (value ",I4,") is too small.  Consider increasing max_subaxes.")') max_subaxes
          CALL error_mesg('diag_axis_mod::diag_subaxes_init', emsg, fatal)
       END IF
       nsub_axis = num_subaxes(axis)
       axes(axis)%start(nsub_axis) = start_indx
       axes(axis)%end(nsub_axis)   = end_indx
       if ( hasdomain ) axes(axis)%subaxis_domain2(nsub_axis) = domain_2d
    END IF

    ! Create new name for the subaxis from name of parent axis
    ! If subaxis already exists, get the index and return
    IF(subaxis_set) THEN
       IF ( axes(axis)%set > 0 ) THEN
          diag_subaxes_init = get_axis_num(name, set_name=trim(axis_sets(axes(axis)%set)))
       ELSE
          diag_subaxes_init = get_axis_num(name)
       END IF
    ELSE
       ! get a new index for subaxis
       !::sdu:: Need a check to allow larger numbers in the index number.
       WRITE (nsub_name,'(I2.2)') nsub_axis
       name = trim(axes(axis)%name)//'_sub'//trim(nsub_name)
       axes(axis)%subaxis_name(nsub_axis) = name
       long_name = trim(axes(axis)%long_name)
       units = trim(axes(axis)%units)
       cart_name = trim(axes(axis)%cart_name)
       direction = axes(axis)%direction
       IF (axes(axis)%set > 0) THEN
          diag_subaxes_init =  diag_axis_init(trim(name), subdata, trim(units), trim(cart_name), trim(long_name),&
               & set_name=trim(axis_sets(axes(axis)%set)), direction=direction, domain2=domain_2d)
       ELSE
          diag_subaxes_init =  diag_axis_init(trim(name), subdata, trim(units), trim(cart_name), trim(long_name),&
               & direction=direction, domain2=domain_2d)
       END IF
    END IF
  END FUNCTION diag_subaxes_init
  ! </FUNCTION>

  ! <SUBROUTINE NAME="get_diag_axis">
  !   <OVERVIEW>
  !     Return information about the axis with index ID
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE get_diag_axis(id, name, units, long_name, cart_name,
  !          direction, edges, Domain, data, num_attributes, attributes)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Return information about the axis with index ID
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  !   <OUT NAME="name" TYPE="CHARACTER(len=*)">Short name for axis</OUT>
  !   <OUT NAME="units" TYPE="CHARACTER(len=*)">Units for axis</OUT>
  !   <OUT NAME="long_name" TYPE="CHARACTER(len=*)">Long name for axis</OUT>
  !   <OUT NAME="cart_name" TYPE="CHARACTER(len=*)">
  !     Cartesian axis ("x", "y", "z", "t", "u").
  !   </OUT>
  !   <OUT NAME="direction" TYPE="INTEGER">
  !     Direction of data. (See <TT>diag_axis_init</TT> for a description of
  !     allowed values)
  !   </OUT>
  !   <OUT NAME="edges" TYPE="INTEGER">
  !     Axis ID for the previously defined "edges axis".
  !   </OUT>
  !   <OUT NAME="Domain" TYPE="TYPE(domain1d)" />
  !   <OUT NAME="DomainU" TYPE="TYPE(domainUG)" />
  !   <OUT NAME="data" TYPE="REAL, DIMENSION(:)">
  !     Array of coordinate values for this axis.
  !   </OUT>
  SUBROUTINE get_diag_axis(id, name, units, long_name, cart_name,&
       & direction, edges, Domain, DomainU, DATA, num_attributes, attributes)
    CHARACTER(len=*), INTENT(out) :: name, units, long_name, cart_name
    INTEGER, INTENT(in) :: id
    TYPE(domain1d), INTENT(out) :: domain
    TYPE(domainug), INTENT(out) :: domainu
    INTEGER, INTENT(out) :: direction, edges
    REAL, DIMENSION(:), INTENT(out) :: data
    INTEGER, INTENT(out), OPTIONAL :: num_attributes
    TYPE(diag_atttype), ALLOCATABLE, DIMENSION(:), INTENT(out), OPTIONAL :: attributes

    INTEGER :: i, j, istat

    CALL valid_id_check(id, 'get_diag_axis')
    name      = axes(id)%name
    units     = axes(id)%units
    long_name = axes(id)%long_name
    cart_name = axes(id)%cart_name
    direction = axes(id)%direction
    edges     = axes(id)%edges
    domain    = axes(id)%Domain
    domainu   = axes(id)%DomainUG
    IF ( axes(id)%length > SIZE(DATA(:)) ) THEN
       ! <ERROR STATUS="FATAL">array data is too small.</ERROR>
       CALL error_mesg('diag_axis_mod::get_diag_axis', 'array data is too small', fatal)
    ELSE
       DATA(1:axes(id)%length) = axes(id)%data(1:axes(id)%length)
    END IF
    IF ( PRESENT(num_attributes) ) THEN
       num_attributes = axes(id)%num_attributes
    END IF
    IF ( PRESENT(attributes) ) THEN
       IF ( _allocated(axes(id)%attributes) ) THEN
          IF ( ALLOCATED(attributes) ) THEN
             ! If allocate, make sure attributes is large enough to hold Axis(id)%attributes
             IF ( axes(id)%num_attributes .GT. SIZE(attributes(:)) ) THEN
                CALL error_mesg('diag_axis_mod::get_diag_axis', 'array attribute is too small', fatal)
             END IF
          ELSE
             ! Allocate attributes
             ALLOCATE(attributes(axes(id)%num_attributes), stat=istat)
             IF ( istat .NE. 0 ) THEN
                CALL error_mesg('diag_axis_mod::get_diag_axis', 'Unable to allocate memory for attribute', fatal)
             END IF
          END IF
          DO i=1, axes(id)%num_attributes
             ! Unallocate all att arrays in preparation for new data
             IF ( _allocated(attributes(i)%fatt) ) THEN
                DEALLOCATE(attributes(i)%fatt)
             END IF
             IF ( _allocated(attributes(i)%iatt) ) THEN
                DEALLOCATE(attributes(i)%iatt)
             END IF

             ! Copy in attribute data
             attributes(i)%type = axes(id)%attributes(i)%type
             attributes(i)%len = axes(id)%attributes(i)%len
             attributes(i)%name = axes(id)%attributes(i)%name
             attributes(i)%catt = axes(id)%attributes(i)%catt
             ! Allocate fatt arrays (if needed), and copy in data
             IF ( _allocated(axes(id)%attributes(i)%fatt) ) THEN
                ALLOCATE(attributes(i)%fatt(SIZE(axes(id)%attributes(i)%fatt(:))), stat=istat)
                IF ( istat .NE. 0 ) THEN
                   CALL error_mesg('diag_axis_mod::get_diag_axis', 'Unable to allocate memory for attribute%fatt', fatal)
                END IF
                DO j=1, SIZE(attributes(i)%fatt(:))
                   attributes(i)%fatt(j) = axes(id)%attributes(i)%fatt(j)
                END DO
             END IF
             ! Allocate iatt arrays (if needed), and copy in data
             IF ( _allocated(axes(id)%attributes(i)%iatt) ) THEN
                ALLOCATE(attributes(i)%iatt(SIZE(axes(id)%attributes(i)%iatt(:))), stat=istat)
                IF ( istat .NE. 0 ) THEN
                   CALL error_mesg('diag_axis_mod::get_diag_axis', 'Unable to allocate memory for attribute%iatt', fatal)
                END IF
                DO j=1, SIZE(attributes(i)%iatt(:))
                   attributes(i)%iatt(j) = axes(id)%attributes(i)%iatt(j)
                END DO
             END IF
          END DO
       END IF
    END IF
  END SUBROUTINE get_diag_axis
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="get_diag_axis_cart">
  !   <OVERVIEW>
  !     Return the axis cartesian.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE get_diag_axis_cart(id, cart_name)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Return the axis cartesian ('X', 'Y', 'Z' or 'T') for the axis ID given.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  !   <OUT NAME="cart_name" TYPE="CHARACTER(len=*)">Cartesian axis</OUT>
  SUBROUTINE get_diag_axis_cart(id, cart_name)
    INTEGER, INTENT(in)           :: id
    CHARACTER(len=*), INTENT(out) :: cart_name

    CALL valid_id_check(id, 'get_diag_axis_cart')
    cart_name = axes(id)%cart_name
  END SUBROUTINE get_diag_axis_cart
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="get_diag_axis_data">
  !   <OVERVIEW>
  !     Return the axis data.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE get_diag_axis_data(id, data)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Return the axis data for the axis ID given.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  !   <OUT NAME="data" TYPE="REAL, DIMENSION(:)">Axis data</OUT>
  SUBROUTINE get_diag_axis_data(id, DATA)
    INTEGER, INTENT(in) :: id
    REAL, DIMENSION(:), INTENT(out) :: data

    CALL valid_id_check(id, 'get_diag_axis_data')
    IF (axes(id)%length > SIZE(DATA(:))) THEN
       ! <ERROR STATUS="FATAL">array data is too small</ERROR>
       CALL error_mesg('diag_axis_mod::get_diag_axis_data', 'array data is too small', fatal)
    ELSE
       DATA(1:axes(id)%length) = axes(id)%data
    END IF
  END SUBROUTINE get_diag_axis_data
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="get_diag_axis_name">
  !   <OVERVIEW>
  !     Return the short name of the axis.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE get_diag_axis_name (id, name)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Return the short name for the axis ID given.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  !   <OUT NAME="name" TYPE="CHARACTER(len=*)">Axis short name</OUT>
  SUBROUTINE get_diag_axis_name(id, name)
    INTEGER         , INTENT(in)  :: id
    CHARACTER(len=*), INTENT(out) :: name

    CALL valid_id_check(id, 'get_diag_axis_name')
    name = axes(id)%name
  END SUBROUTINE get_diag_axis_name
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="get_diag_axis_domain_name">
  !   <OVERVIEW>
  !     Return the name of the axis' domain
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE get_diag_axis_domain_name(id, name)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Retruns the name of the axis' domain.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  !   <OUT NAME="name" TYPE="CHARACTER(len=*)">Axis' domain name</OUT>
  SUBROUTINE get_diag_axis_domain_name(id, name)
    INTEGER, INTENT(in) :: id
    CHARACTER(len=*), INTENT(out) :: name

    CALL valid_id_check(id, 'get_diag_axis_domain_name')
    name = mpp_get_domain_name(axes(id)%domain2)
  END SUBROUTINE get_diag_axis_domain_name
  ! </SUBROUTINE>

  ! <FUNCTION NAME="get_axis_length">
  !   <OVERVIEW>
  !     Return the length of the axis.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     INTEGER FUNCTION get_axis_length(id)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Return the length of the axis ID given.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  INTEGER FUNCTION get_axis_length(id)
    INTEGER, INTENT(in) :: id

    INTEGER :: length

    CALL valid_id_check(id, 'get_axis_length')
    IF ( axes(id)%Domain .NE. null_domain1d ) THEN
       CALL mpp_get_compute_domain(axes(id)%Domain,size=length)
       !---one extra point is needed for some case. ( like symmetry domain )
       get_axis_length = length + axes(id)%shift
    ELSE
       get_axis_length = axes(id)%length
    END IF
  END FUNCTION get_axis_length
  ! </FUNCTION>

  ! <FUNCTION NAME="get_axis_aux">
  !   <OVERVIEW>
  !     Return the auxiliary name for the axis.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     CHARACTER(len=128) FUNCTION get_axis_aux(id)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Returns the auxiliary name for the axis.  The only possible values for
  !     the auxiliary names is <TT>geolon_t</TT> or <TT>geolat_t</TT>.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  CHARACTER(len=128) FUNCTION get_axis_aux(id)
    INTEGER, INTENT(in) :: id

    CALL valid_id_check(id, 'get_axis_aux')
    get_axis_aux =  axes(id)%aux
  END FUNCTION get_axis_aux
  ! </FUNCTION>

  ! <FUNCTION NAME="get_axis_reqfld">
  !   <OVERVIEW>
  !     Return the required field names for the axis.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     CHARACTER(len=128) FUNCTION get_axis_reqfld(id)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Returns the required field names for the axis.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  CHARACTER(len=128) FUNCTION get_axis_reqfld(id)
    INTEGER, INTENT(in) :: id

    CALL valid_id_check(id, 'get_axis_reqfld')
    get_axis_reqfld =  axes(id)%req
  END FUNCTION get_axis_reqfld
  ! </FUNCTION>

  ! <FUNCTION NAME="get_axis_global_length">
  !   <OVERVIEW>
  !     Return the global length of the axis.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     INTEGER FUNCTION get_axis_global_length (id)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Returns the global length of the axis ID given.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  INTEGER FUNCTION get_axis_global_length(id)
    INTEGER, INTENT(in) :: id

    CALL valid_id_check(id, 'get_axis_global_length')
    get_axis_global_length = axes(id)%length
  END FUNCTION get_axis_global_length
  ! </FUNCTION>

  ! <FUNCTION NAME="get_tile_count">
  !   <OVERVIEW>
  !     Return the tile count for the axis.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     INTEGER FUNCTION get_tile_count (ids)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Return the tile count for the axis IDs given.
  !   </DESCRIPTION>
  !   <IN NAME="ids" TYPE="INTEGER, DIMENSION(:)">
  !     Axis IDs.  Possible dimensions: 1 <= <TT>size(ids(:))</TT> <= 4.
  !   </IN>
  INTEGER FUNCTION get_tile_count(ids)
    INTEGER, DIMENSION(:), INTENT(in) :: ids

    INTEGER :: i, id, flag

    IF ( SIZE(ids(:)) < 1 ) THEN
       ! <ERROR STATUS="FATAL">input argument has incorrect size.</ERROR>
       CALL error_mesg('diag_axis_mod::get_tile_count', 'input argument has incorrect size', fatal)
    END IF
    get_tile_count = 1
    flag = 0
    DO i = 1, SIZE(ids(:))
       id = ids(i)
       CALL valid_id_check(id, 'get_tile_count')
       IF ( axes(id)%cart_name == 'X' .OR.  &
            axes(id)%cart_name == 'Y' ) flag = flag + 1
       !     --- both x/y axes found ---
       IF ( flag == 2 ) THEN
          get_tile_count = axes(id)%tile_count
          EXIT
       END IF
    END DO
  END FUNCTION get_tile_count
  ! </FUNCTION>

  ! <FUNCTION NAME="get_domain1d">
  !   <OVERVIEW>
  !     Return the 1D domain.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     TYPE(domain1d) FUNCTION get_domain1d(id)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Retrun the 1D domain for the axis ID given.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  TYPE(domain1d) function get_domain1d(id)
    INTEGER, INTENT(in) :: id

    CALL valid_id_check(id, 'get_domain1d')
    IF (axes(id)%Domain .NE. null_domain1d) THEN
       get_domain1d = axes(id)%Domain
    ELSE
       get_domain1d = null_domain1d
    ENDIF
  END FUNCTION get_domain1d
  ! </FUNCTION>

  ! <FUNCTION NAME="get_domain2d">
  !   <OVERVIEW>
  !     Return the 2D domain.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     TYPE(domain2d) FUNCTION get_domain2d(ids)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Return the 2D domain for the axis IDs given.
  !   </DESCRIPTION>
  !   <IN NAME="ids" TYPE="INTEGER, DIMENSION(:)">
  !     Axis IDs.  Possible dimensions: 1 <= <TT>size(ids(:))</TT> <= 4.
  !   </IN>
  TYPE(domain2d) function get_domain2d(ids)
    INTEGER, DIMENSION(:), INTENT(in) :: ids

    INTEGER :: i, id, flag

    IF ( SIZE(ids(:)) < 1 ) THEN
       ! <ERROR STATUS="FATAL">input argument has incorrect size.</ERROR>
       CALL error_mesg('diag_axis_mod::get_domain2d', 'input argument has incorrect size', fatal)
    END IF
    get_domain2d = null_domain2d
    flag = 0
    DO i = 1, SIZE(ids(:))
       id = ids(i)
       CALL valid_id_check(id, 'get_domain2d')
       IF ( axes(id)%cart_name == 'X' .OR. axes(id)%cart_name == 'Y' ) flag = flag + 1
       !     --- both x/y axes found ---
       IF ( flag == 2 ) THEN
          IF (axes(id)%Domain2 .NE. null_domain2d) get_domain2d = axes(id)%Domain2
          EXIT
       END IF
    END DO
  END FUNCTION get_domain2d

  ! <FUNCTION NAME="get_domainUG">
  !   <OVERVIEW>
  !     Return the UG domain.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     TYPE(domainUG) FUNCTION get_domainUG(id)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Retrun the 1D domain for the axis ID given.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
  TYPE(domainug) function get_domainug(id)
    INTEGER, INTENT(in) :: id

    CALL valid_id_check(id, 'get_domainUG')
    IF (axes(id)%DomainUG .NE. null_domainug) THEN
       get_domainug = axes(id)%DomainUG
    ELSE
       get_domainug = null_domainug
    ENDIF
  END FUNCTION get_domainug

  ! <SUBROUTINE NAME="axis_compatible_check">
  !   <OVERVIEW>
  !     Checks if the axes are compatible
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     INTEGER FUNCTION axis_compatible_check(id)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Checks if the axes are compatible
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis ID</IN>
!----------
!ug support
  function axis_compatible_check(id,varname) result(domain_type)

   !Inputs/Outputs
    integer,dimension(:),intent(in)  :: id          !<The array of axis IDs
    character(*),intent(in),optional :: varname     !<The name of the variable
    integer(INT_KIND)                :: domain_type !<DIAG_AXIS_NODOMAIN = no domain.
                                                    !<DIAG_AXIS_2DDOMAIN = structured domain.
                                                    !<DIAG_AXIS_UGDOMAIN = unstructured domain.

   !Local variables
    logical :: xory          !<XorY set to true if X or Y is found as a cart_name.
    logical :: ug            !<UG set to true if U is found as a cart_name.
    integer :: n             !<Looping index.
    logical :: uses_domain2d !<True if an axis is associated with a 2D domain.
    logical :: uses_domainug !<True if an axis is associated with an unstructured domain.

   !Initialize flags.
    xory = .false.
    ug = .false.
    uses_domain2d = .false.
    uses_domainug = .false.

   !Make sure that valid set of axes was passed, and determine the domain-type
   !associated with the axes.
    do n = 1,size(id)
        call valid_id_check(id(n), &
                            "axis_compatible_check")
        if (axes(id(n))%cart_name .eq. "X" .or. &
            axes(id(n))%cart_name .eq. "Y") then
            xory = .true.
        elseif (axes(id(n))%cart_name .eq. "U") then
            ug = .true.
        endif
        if (axes(id(n))%Domain2 .ne. null_domain2d) then
            uses_domain2d = .true.
        elseif (axes(id(n))%DomainUG .ne. null_domainug) then
            uses_domainug = .true.
        endif
    enddo
    if (ug .and. xory) then
        if (present(varname)) then
            call error_mesg("axis_compatible_check", &
                            "Can not use an unstructured grid with a "// &
                            "horizontal cartesian coordinate for the field " &
                            //trim(varname), &
                            fatal)
        else
            call error_mesg("axis_compatible_check", &
                            "Can not use an unstructured grid with a horizontal "// &
                            "cartesian coordinate", &
                            fatal)
        endif
    endif
    if (uses_domain2d .and. uses_domainug) then
        if (present(varname)) then
            call error_mesg("axis_compatible_check", &
                            "Can not use an unstructured grid with a"// &
                            "structured grid for the field "//trim(varname), &
                            fatal)
        else
            call error_mesg("axis_compatible_check", &
                            "Can not use an unstructured grid with a"// &
                            "structured grid.", &
                            fatal)
        endif
    endif
    if (uses_domain2d) then
        domain_type = diag_axis_2ddomain
    elseif (uses_domainug) then
        domain_type = diag_axis_ugdomain
    else
        domain_type = diag_axis_nodomain
    endif

    return
  end function axis_compatible_check
!----------

  ! </FUNCTION>

  ! <SUBROUTINE NAME="get_axes_shift">
  !   <OVERVIEW>
  !     Return the value of the shift.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE get_axes_shift(ids, ishift, jshift)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Return the value of the shift for the axis IDs given.
  !   </DESCRIPTION>
  !   <IN NAME="ids" TYPE="INTEGER, DIMENSION(:)">
  !     Axis IDs.  Possible dimensions: 1 <= <TT>size(ids(:))</TT> <= 4
  !   </IN>
  !   <OUT NAME="ishift" TYPE="INTEGER">X shift value.</OUT>
  !   <OUT NAME="jshift" TYPE="INTEGER">Y shift value.</OUT>
  SUBROUTINE get_axes_shift(ids, ishift, jshift)
    INTEGER, DIMENSION(:), INTENT(in) :: ids
    INTEGER, INTENT(out) :: ishift, jshift

    INTEGER :: i, id

    !-- get the value of the shift.
    ishift = 0
    jshift = 0
    DO i = 1, SIZE(ids(:))
       id = ids(i)
       CALL valid_id_check(id, 'get_axes_shift')
       SELECT CASE (axes(id)%cart_name)
       CASE ( 'X' )
          ishift = axes(id)%shift
       CASE ( 'Y' )
          jshift = axes(id)%shift
       END SELECT
    END DO
  END SUBROUTINE get_axes_shift
  ! </SUBROUTINE>

  ! <PRIVATE>
  ! <FUNCTION NAME="get_axis_num">
  !   <OVERVIEW>
  !     Returns index into axis table corresponding to a given axis name.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     INTEGER FUNCTION get_axis_num(axis_name, set_name)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Returns index into axis table corresponding to a given axis name.
  !   </DESCRIPTION>
  !   <IN NAME="axis_name" TYPE="CHARACTER(len=*)">Axis name.</IN>
  !   <IN NAME="set_name" TYPE="CHARACTER(len=*), OPTIONAL">Set name.</IN>
  INTEGER FUNCTION get_axis_num(axis_name, set_name)
    CHARACTER(len=*), INTENT(in) :: axis_name
    CHARACTER(len=*), INTENT(in), OPTIONAL :: set_name

    INTEGER :: set, n

    IF ( PRESENT(set_name) ) THEN
       set = get_axis_set_num(trim(set_name))
    ELSE
       set = 0
    END IF
    get_axis_num = 0
    DO n = 1, num_def_axes
       IF ( trim(axis_name) == trim(axes(n)%name) .AND. axes(n)%set == set ) THEN
          get_axis_num = n
          RETURN
       END IF
    END DO
  END FUNCTION get_axis_num
  ! </FUNCTION>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <FUNCTION NAME="get_axis_set_num">
  !   <OVERVIEW>
  !     Returns index in axis set table corresponding to a given axis set name.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     INTEGER FUNCTION get_axis_set_num(set_name)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Returns index in axis set table corresponding to a given axis set name.
  !   </DESCRIPTION>
  !   <IN NAME="set_name" TYPE="CHARACTER(len=*)">Set name.</IN>
  INTEGER FUNCTION get_axis_set_num (set_name)
    CHARACTER(len=*), INTENT(in) :: set_name

    INTEGER :: iset

    get_axis_set_num = 0
    DO iset = 1, num_axis_sets
       IF ( set_name == axis_sets(iset) ) THEN
          get_axis_set_num = iset
          RETURN
       END IF
    END DO
  END FUNCTION get_axis_set_num
  ! </FUNCTION>
  ! </PRIVATE>

  ! <PRIVATE>
  ! <SUBROUTINE NAME="valid_id_check">
  !   <OVERVIEW>
  !     Check to see if the axis id is a vaild id.
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE valid_id_check(id, routine_name)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Check to see if the given axis id is a valid id.  If the axis id is invalid,
  !     call a FATAL error.  If the ID is valid, just return.
  !   </DESCRIPTION>
  !   <IN NAME="id" TYPE="INTEGER">Axis id to check for validity</IN>
  !   <IN NAME="routine_name" TYPE="CHARACTER(len=*)">Name of the subroutine checking for a valid axis id.</IN>
  SUBROUTINE valid_id_check(id, routine_name)
    INTEGER, INTENT(in) :: id
    CHARACTER(len=*), INTENT(in) :: routine_name

    CHARACTER(len=5) :: emsg

    IF ( id < 1 .OR. id > num_def_axes) THEN
       ! <ERROR STATUS="FATAL">
       !   Illegal value for axis used (value <VALUE>).
       ! </ERROR>
       WRITE (emsg, '(I2)') id
       CALL error_mesg('diag_axis_mod::'//trim(routine_name),&
            & 'Illegal value for axis_id used (value '//trim(emsg)//').', fatal)
    END IF
  END SUBROUTINE valid_id_check
  ! </SUBROUTINE>
  ! </PRIVATE>

  SUBROUTINE diag_axis_attribute_init(diag_axis_id, name, type, cval, ival, rval)
    INTEGER, INTENT(in) :: diag_axis_id !< input field ID, obtained from diag_axis_mod::diag_axis_init.
    CHARACTER(len=*) :: name !< Name of the attribute
    INTEGER, INTENT(in) :: type !< NetCDF type (NF_FLOAT, NF_INT, NF_CHAR)
    CHARACTER(len=*), INTENT(in), OPTIONAL :: cval !< Character string attribute value
    INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: ival !< Integer attribute value(s)
    REAL, DIMENSION(:), INTENT(in), OPTIONAL :: rval !< Real attribute value(s)

    INTEGER :: istat, length, i, j, this_attribute, out_field
    CHARACTER(len=1024) :: err_msg

    IF ( .NOT.first_send_data_call ) THEN
       ! Call error due to unable to add attribute after send_data called
       ! <ERROR STATUS="FATAL">
       !   Attempting to add attribute <name> to axis <axis_name>
       !   after first send_data call.  Too late.
       ! </ERROR>
       CALL error_mesg('diag_manager_mod::diag_axis_add_attribute', 'Attempting to add attribute "'&
            &//trim(name)//'" to axis after first send_data call.  Too late.', fatal)
    END IF

    ! Simply return if diag_axis_id <= 0 --- not registered
    IF ( diag_axis_id .LE. 0 ) THEN
       RETURN
    ELSE IF ( diag_axis_id .GT. num_def_axes ) THEN
       ! Call error axis_id greater than num_def_axes.  Axis is unknown
       ! <ERROR STATUS="FATAL">
       !   Attempting to add attribute <name> to axis ID <axis_ID>, however ID unknown.
       ! </ERROR>
       WRITE(err_msg, '(I5)') diag_axis_id
       CALL error_mesg('diag_manager_mod::diag_axis_add_attribute', 'Attempting to add attribute "'&
            &//trim(name)//'" to axis ID "'//trim(err_msg)//'", however ID unknown.', fatal)

    ELSE
       ! Allocate memory for the attributes
       CALL attribute_init_axis(axes(diag_axis_id))

       ! Check if attribute already exists
       this_attribute = 0
       DO i=1, axes(diag_axis_id)%num_attributes
          IF ( trim(axes(diag_axis_id)%attributes(i)%name) .EQ. trim(name) ) THEN
             this_attribute = i
             EXIT
          END IF
       END DO

       IF ( this_attribute.NE.0 .AND. (type.EQ.nf90_int .OR. type.EQ.nf90_float) ) THEN
          ! <ERROR STATUS="FATAL">
          !   Attribute <name> already defined for axis <axis_name>.
          !   Contact the developers
          ! </ERROR>
          CALL error_mesg('diag_manager_mod::diag_axis_add_attribute',&
               & 'Attribute "'//trim(name)//'" already defined for axis "'&
               &//trim(axes(diag_axis_id)%name)//'".  Contact the developers.', fatal)
       ELSE IF ( this_attribute.NE.0 .AND. type.EQ.nf90_char .AND. debug_diag_manager ) THEN
          ! <ERROR STATUS="NOTE">
          !   Attribute <name> already defined for axis <axis_name>.
          !   Prepending.
          ! </ERROR>
          CALL error_mesg('diag_manager_mod::diag_axis_add_attribute',&
               & 'Attribute "'//trim(name)//'" already defined for axis"'&
               &//trim(axes(diag_axis_id)%name)//'".  Prepending.', note)
       ELSE
          ! Defining a new attribute
          ! Increase the number of field attributes
          this_attribute = axes(diag_axis_id)%num_attributes + 1
          ! Checking to see if num_attributes == max_axis_attributes, and return error message
          IF ( this_attribute .GT. max_axis_attributes ) THEN
             ! <ERROR STATUS="FATAL">
             !   Number of attributes exceeds max_axis_attributes for attribute <name> for axis <axis_name>.
             !   Increase diag_manager_nml:max_axis_attributes.
             ! </ERROR>
             CALL error_mesg('diag_manager_mod::diag_axis_add_attribute',&
                  & 'Number of attributes exceeds max_axis_attributes for attribute "'&
                  & //trim(name)//'" for axis "'//trim(axes(diag_axis_id)%name)&
                  & //'".  Increase diag_manager_nml:max_axis_attributes.',&
                  & fatal)
          ELSE
             axes(diag_axis_id)%num_attributes = this_attribute
             ! Set name and type
             axes(diag_axis_id)%attributes(this_attribute)%name = name
             axes(diag_axis_id)%attributes(this_attribute)%type = type
             ! Initialize catt to a blank string, as len_trim doesn't always work on an uninitialized string
             axes(diag_axis_id)%attributes(this_attribute)%catt = ''
          END IF
       END IF

       SELECT CASE (type)
       CASE (nf90_int)
          IF ( .NOT.PRESENT(ival) ) THEN
             ! <ERROR STATUS="FATAL">
             !   Number type claims INTEGER, but ival not present for attribute <name> for axis <axis_name>.
             !   Contact the developers.
             ! </ERROR>
             CALL error_mesg('diag_manager_mod::diag_axis_add_attribute',&
                  & 'Attribute type claims INTEGER, but ival not present for attribute "'&
                  & //trim(name)//'" for axis "'//trim(axes(diag_axis_id)%name)&
                  & //'". Contact the developers.', fatal)
          END IF
          length = SIZE(ival)
          ! Allocate iatt(:) to size of ival
          ALLOCATE(axes(diag_axis_id)%attributes(this_attribute)%iatt(length), stat=istat)
          IF ( istat.NE.0 ) THEN
             ! <ERROR STATUS="FATAL">
             !   Unable to allocate iatt for attribute <name> for axis <axis_name>
             ! </ERROR>
             CALL error_mesg('diag_manager_mod::diag_axis_add_attribute', 'Unable to allocate iatt for attribute "'&
                  & //trim(name)//'" for axis "'//trim(axes(diag_axis_id)%name)//'"', fatal)
          END IF
          ! Set remaining fields
          axes(diag_axis_id)%attributes(this_attribute)%len = length
          axes(diag_axis_id)%attributes(this_attribute)%iatt = ival
       CASE (nf90_float)
          IF ( .NOT.PRESENT(rval) ) THEN
             ! <ERROR STATUS="FATAL">
             !   Attribute type claims REAL, but rval not present for attribute <name> for axis <axis_name>.
             !   Contact the developers.
             ! </ERROR>
             CALL error_mesg('diag_manager_mod::diag_axis_add_attribute',&
                  & 'Attribute type claims REAL, but rval not present for attribute "'&
                  & //trim(name)//'" for axis "'//trim(axes(diag_axis_id)%name)&
                  & //'". Contact the developers.', fatal)
          END IF
          length = SIZE(rval)
          ! Allocate iatt(:) to size of rval
          ALLOCATE(axes(diag_axis_id)%attributes(this_attribute)%fatt(length), stat=istat)
          IF ( istat.NE.0 ) THEN
             ! <ERROR STATUS="FATAL">
             !   Unable to allocate fatt for attribute <name> for axis <axis_name>
             ! </ERROR>
             CALL error_mesg('diag_manager_mod::diag_axis_add_attribute', 'Unable to allocate fatt for attribute "'&
                  & //trim(name)//'" for axis "'//trim(axes(diag_axis_id)%name)&
                  & //'"', fatal)
          END IF
          ! Set remaining fields
          axes(diag_axis_id)%attributes(this_attribute)%len = length
          axes(diag_axis_id)%attributes(this_attribute)%fatt = rval
       CASE (nf90_char)
          IF ( .NOT.PRESENT(cval) ) THEN
             ! <ERROR STATUS="FATAL">
             !   Attribute type claims CHARACTER, but cval not present for attribute <name> for axis <axis_name>.
             !   Contact the developers.
             ! </ERROR>
             CALL error_mesg('diag_manager_mod::diag_axis_add_attribute',&
                  & 'Attribute type claims CHARACTER, but cval not present for attribute "'&
                  & //trim(name)//'" for axis "'//trim(axes(diag_axis_id)%name)&
                  & //'". Contact the developers.', fatal)
          END IF
          CALL prepend_attribute_axis(axes(diag_axis_id), trim(name), trim(cval))
       CASE default
          ! <ERROR STATUS="FATAL">
          !   Unknown attribute type for attribute <name> for axis <axis_name>.
          !   Contact the developers.
          ! </ERROR>
          CALL error_mesg('diag_manager_mod::diag_axis_add_attribute', 'Unknown attribute type for attribute "'&
               & //trim(name)//'" for axis "'//trim(axes(diag_axis_id)%name)&
               & //'". Contact the developers.', fatal)
       END SELECT
    END IF
  END SUBROUTINE diag_axis_attribute_init

  ! <SUBROUTINE NAME="diag_axis_add_attribute_scalar_r" INTERFACE="diag_axis_add_attribute">
  !   <IN NAME="diag_axis_id" TYPE="INTEGER" />
  !   <IN NAME="att_name" TYPE="CHARACTER(len=*)" />
  !   <IN NAME="att_value" TYPE="REAL" />
  SUBROUTINE diag_axis_add_attribute_scalar_r(diag_axis_id, att_name, att_value)
    INTEGER, INTENT(in) :: diag_axis_id
    CHARACTER(len=*), INTENT(in) :: att_name
    REAL, INTENT(in) :: att_value

    CALL diag_axis_add_attribute_r1d(diag_axis_id, att_name, (/ att_value /))
  END SUBROUTINE diag_axis_add_attribute_scalar_r
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="diag_axis_add_attribute_scalar_i" INTERFACE="diag_axis_add_attribute">
  !   <IN NAME="diag_axis_id" TYPE="INTEGER" />
  !   <IN NAME="att_name" TYPE="CHARACTER(len=*)" />
  !   <IN NAME="att_value" TYPE="INTEGER" />
  SUBROUTINE diag_axis_add_attribute_scalar_i(diag_axis_id, att_name, att_value)
    INTEGER, INTENT(in) :: diag_axis_id
    CHARACTER(len=*), INTENT(in) :: att_name
    INTEGER, INTENT(in) :: att_value

    CALL diag_axis_add_attribute_i1d(diag_axis_id, att_name, (/ att_value /))
  END SUBROUTINE diag_axis_add_attribute_scalar_i
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="diag_axis_add_attribute_scalar_c" INTERFACE="diag_axis_add_attribute">
  !   <IN NAME="diag_axis_id" TYPE="INTEGER" />
  !   <IN NAME="att_name" TYPE="CHARACTER(len=*)" />
  !   <IN NAME="att_value" TYPE="CHARACTER(len=*)" />
  SUBROUTINE diag_axis_add_attribute_scalar_c(diag_axis_id, att_name, att_value)
    INTEGER, INTENT(in) :: diag_axis_id
    CHARACTER(len=*), INTENT(in) :: att_name
    CHARACTER(len=*), INTENT(in) :: att_value

    CALL diag_axis_attribute_init(diag_axis_id, att_name, nf90_char, cval=att_value)
  END SUBROUTINE diag_axis_add_attribute_scalar_c
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="diag_axis_add_attribute_r1d" INTERFACE="diag_axis_add_attribute">
  !   <IN NAME="diag_axis_id" TYPE="INTEGER" />
  !   <IN NAME="att_name" TYPE="CHARACTER(len=*)" />
  !   <IN NAME="att_value" TYPE="REAL, DIMENSION(:)" />
  SUBROUTINE diag_axis_add_attribute_r1d(diag_axis_id, att_name, att_value)
    INTEGER, INTENT(in) :: diag_axis_id
    CHARACTER(len=*), INTENT(in) :: att_name
    REAL, DIMENSION(:), INTENT(in) :: att_value

    INTEGER :: num_attributes, len
    CHARACTER(len=512) :: err_msg

    CALL diag_axis_attribute_init(diag_axis_id, att_name, nf90_float, rval=att_value)
  END SUBROUTINE diag_axis_add_attribute_r1d
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="diag_axis_add_attribute_i1d" INTERFACE="diag_axis_add_attribute">
  !   <IN NAME="diag_axis_id" TYPE="INTEGER" />
  !   <IN NAME="att_name" TYPE="CHARACTER(len=*)" />
  !   <IN NAME="att_value" TYPE="INTEGER, DIMENSION(:)" />
  SUBROUTINE diag_axis_add_attribute_i1d(diag_axis_id, att_name, att_value)
    INTEGER, INTENT(in) :: diag_axis_id
    CHARACTER(len=*), INTENT(in) :: att_name
    INTEGER, DIMENSION(:), INTENT(in) :: att_value

    CALL diag_axis_attribute_init(diag_axis_id, att_name, nf90_int, ival=att_value)
  END SUBROUTINE diag_axis_add_attribute_i1d
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="attribute_init_axis" INTERFACE="attribute_init">
  !   <OVERVIEW>
  !     Allocates the atttype in out_axis
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE attribute_init(out_axis, err_msg)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Allocates memory in out_file for the attributes.  Will <TT>FATAL</TT> if err_msg is not included
  !     in the subroutine call.
  !   </DESCRIPTION>
  !   <INOUT NAME="out_file" TYPE="TYPE(file_type)">output file to allocate memory for attribute</INOUT>
  !   <OUT NAME="err_msg" TYPE="CHARACTER(len=*), OPTIONAL">Error message, passed back to calling function</OUT>
  SUBROUTINE attribute_init_axis(out_axis, err_msg)
    TYPE(diag_axis_type), INTENT(inout) :: out_axis
    CHARACTER(LEN=*), INTENT(out), OPTIONAL :: err_msg

    INTEGER :: istat

    ! Initialize err_msg
    IF ( PRESENT(err_msg) ) err_msg = ''

    ! Allocate memory for the attributes
    IF ( .NOT._allocated(out_axis%attributes) ) THEN
       ALLOCATE(out_axis%attributes(max_axis_attributes), stat=istat)
       IF ( istat.NE.0 ) THEN
          ! <ERROR STATUS="FATAL">
          !   Unable to allocate memory for diag axis attributes
          ! </ERROR>
          IF ( fms_error_handler('diag_util_mod::attribute_init_axis', 'Unable to allocate memory for diag axis attributes', err_msg) ) THEN
             RETURN
          END IF
       ELSE
          ! Set equal to 0.  It will be increased when attributes added
          out_axis%num_attributes = 0
       END IF
    END IF
  END SUBROUTINE attribute_init_axis
  ! </SUBROUTINE>

  ! <SUBROUTINE NAME="prepend_attribute_axis">
  !   <OVERVIEW>
  !     Prepends the attribute value to an already existing attribute.  If the
  !     attribute isn't yet defined, then creates a new attribute
  !   </OVERVIEW>
  !   <TEMPLATE>
  !     SUBROUTINE prepend_attribute(out_file, attribute, prepend_value)
  !   </TEMPLATE>
  !   <DESCRIPTION>
  !     Checks if the attribute already exists for <TT>out_axis</TT>.  If it exists,
  !     then prepend the <TT>prepend_value</TT>, otherwise, create the attribute
  !     with the <TT>prepend_value</TT>. <TT>err_msg</TT> indicates no duplicates found.
  !   </DESCRIPTION>
  !   <INOUT NAME="out_axis" TYPE="TYPE(diag_axis_type)">diagnostic axis that will get the attribute</INOUT>
  !   <IN NAME="att_name" TYPE="CHARACTER(len=*)">Name of the attribute</IN>
  !   <IN NAME="prepend_value" TYPE="CHARACTER(len=*)">Value to prepend</IN>
  !   <OUT NAME="err_msg" TYPE="CHARACTER(len=*), OPTIONAL">Error message, passed back to calling routine</OUT>
  SUBROUTINE prepend_attribute_axis(out_axis, att_name, prepend_value, err_msg)
    TYPE(diag_axis_type), INTENT(inout) :: out_axis
    CHARACTER(len=*), INTENT(in) :: att_name, prepend_value
    CHARACTER(len=*), INTENT(out) , OPTIONAL :: err_msg

    INTEGER :: length, i, this_attribute
    CHARACTER(len=512) :: err_msg_local

    ! Initialize string variables
    err_msg_local = ''
    IF ( PRESENT(err_msg) ) err_msg = ''

    ! Make sure the attributes for this out file have been initialized
    CALL attribute_init_axis(out_axis, err_msg_local)
    IF ( trim(err_msg_local) .NE. '' ) THEN
       IF ( fms_error_handler('diag_util_mod::prepend_attribute_axis', trim(err_msg_local), err_msg) ) THEN
          RETURN
       END IF
    END IF

    ! Find if attribute exists
    this_attribute = 0
    DO i=1, out_axis%num_attributes
       IF ( trim(out_axis%attributes(i)%name) .EQ. trim(att_name) ) THEN
          this_attribute = i
          EXIT
       END IF
    END DO

    IF ( this_attribute > 0 ) THEN
       IF ( out_axis%attributes(this_attribute)%type .NE. nf90_char ) THEN
          ! <ERROR STATUS="FATAL">
          !   Attribute <name> is not a character attribute.
          ! </ERROR>
          IF ( fms_error_handler('diag_util_mod::prepend_attribute_axis',&
               & 'Attribute "'//trim(att_name)//'" is not a character attribute.',&
               & err_msg) ) THEN
             RETURN
          END IF
       END IF
    ELSE
       ! Defining a new attribute
       ! Increase the number of file attributes
       this_attribute = out_axis%num_attributes + 1
       IF ( this_attribute .GT. max_axis_attributes ) THEN
          ! <ERROR STATUS="FATAL">
          !   Number of attributes exceeds max_axis_attributes for attribute <name>.
          !   Increase diag_manager_nml:max_axis_attributes.
          ! </ERROR>
          IF ( fms_error_handler('diag_util_mod::prepend_attribute_axis',&
               & 'Number of attributes exceeds max_axis_attributes for attribute "'&
               &//trim(att_name)//'".  Increase diag_manager_nml:max_axis_attributes.',&
               & err_msg) ) THEN
             RETURN
          END IF
       ELSE
          out_axis%num_attributes = this_attribute
          ! Set name and type
          out_axis%attributes(this_attribute)%name = att_name
          out_axis%attributes(this_attribute)%type = nf90_char
          ! Initialize catt to a blank string, as len_trim doesn't always work on an uninitialized string
          out_axis%attributes(this_attribute)%catt = ''
       END IF
    END IF

    ! Only add string only if not already defined
    IF ( index(trim(out_axis%attributes(this_attribute)%catt), trim(prepend_value)).EQ.0 ) THEN
       ! Check if new string length goes beyond the length of catt
       length = len_trim(trim(prepend_value)//" "//trim(out_axis%attributes(this_attribute)%catt))
       IF ( length.GT.len(out_axis%attributes(this_attribute)%catt) ) THEN
          ! <ERROR STATUS="FATAL">
          !   Prepend length for attribute <name> is longer than allowed.
          ! </ERROR>
          IF ( fms_error_handler('diag_util_mod::prepend_attribute_file',&
               & 'Prepend length for attribute "'//trim(att_name)//'" is longer than allowed.',&
               & err_msg) ) THEN
             RETURN
          END IF
       END IF
       ! Set files
       out_axis%attributes(this_attribute)%catt =&
            & trim(prepend_value)//' '//trim(out_axis%attributes(this_attribute)%catt)
       out_axis%attributes(this_attribute)%len = length
    END IF
  END SUBROUTINE prepend_attribute_axis
  ! </SUBROUTINE>

  ! given an axis, returns TRUE if the axis uses compression-by-gathering: that is, if
  ! this is an axis for fields on unstructured grid
  logical function axis_is_compressed(id)
    integer, intent(in) :: id

    integer :: i

    CALL valid_id_check(id, 'axis_is_compressed')

    axis_is_compressed = .false.
    if (.not._allocated(axes(id)%attributes)) return
    do i = 1, axes(id)%num_attributes
       if (trim(axes(id)%attributes(i)%name)=='compress') then
          axis_is_compressed = .true.
          return
       endif
    enddo
  end function axis_is_compressed


  ! given an index of compressed-by-gathering axis, return an array of axes used in
  ! compression. It is a fatal error to call it on axis that is not compressed
  subroutine get_compressed_axes_ids(id, r)
    integer, intent(in)  :: id
    integer, intent(out), allocatable :: r(:)

    integer iatt, k, k1, k2, n
    logical :: space

    character(*), parameter :: tag = 'get_compressed_axes_ids'

    CALL valid_id_check(id, tag)

    associate(axis=>axes(id))
    if (.not._allocated(axis%attributes)) call error_mesg(tag, &
       'attempt to get compression dimensions from axis "'//trim(axis%name)//'" which is not compressed (does not have any attributes)', fatal)

    iatt = 0
    do k = 1,axis%num_attributes
       if (trim(axis%attributes(k)%name)=='compress') then
          iatt = k; exit ! from loop
       endif
    enddo

    if (iatt == 0) call error_mesg(tag, &
       'attempt to get compression dimensions from axis "'//trim(axis%name)//&
       '" which is not compressed (does not have "compress" attributes).', fatal)
    if (axis%attributes(iatt)%type/=nf90_char) call error_mesg(tag, &
       'attempt to get compression dimensions from axis "'//trim(axis%name)//&
       '" but the axis attribute "compress" has incorrect type.', fatal)

    ! parse the "compress" attribute
    ! calculate the number of compression axes
    space = .true.; n=0
    do k = 1, len(axis%attributes(iatt)%catt)
       if (space.and.(axis%attributes(iatt)%catt(k:k)/=' ')) then
          n = n+1
       endif
       space = (axis%attributes(iatt)%catt(k:k)==' ')
    enddo

    allocate(r(n))
    ! make array of compression axes indices. Go from the last to the first to get the
    ! array in FORTRAN order: they are listed in "compress" attribute  C order (fastest
    ! dimension last)
    k2 = 0
    do k = n, 1, -1
       do k1 = k2+1, len(axis%attributes(iatt)%catt)
          if (axis%attributes(iatt)%catt(k1:k1)/=' ') exit
       enddo
       do k2 = k1+1, len(axis%attributes(iatt)%catt)
          if (axis%attributes(iatt)%catt(k2:k2)==' ') exit
       enddo
       r(k) = get_axis_num(axis%attributes(iatt)%catt(k1:k2),axis_sets(axis%set))
       if (r(k)<=0) call error_mesg(tag, &
           'compression dimension "'//trim(axis%attributes(iatt)%catt(k1:k2))//&
           '" not found among the axes of set "'//trim(axis_sets(axis%set))//'".', fatal)
    enddo
    end associate ! axis
  end subroutine get_compressed_axes_ids
END MODULE diag_axis_mod