diag_util.F90

!***********************************************************************
!*                   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/>.
!***********************************************************************
!> @defgroup diag_util_mod diag_util_mod
!> @ingroup diag_manager
!! @brief Functions and subroutines necessary for the <TT>diag_manager_mod</TT>.
!! @author Seth Underwood
 
MODULE diag_util_mod
 
use platform_mod
use,intrinsic :: iso_fortran_env, only: real128
use,intrinsic :: iso_c_binding, only: c_double,c_float,c_int64_t, &
                                      c_int32_t,c_int16_t,c_intptr_t
 
  !   <FUTURE>
  !     Make an interface <TT>check_bounds_are_exact</TT> for the subroutines <TT>check_bounds_are_exact_static</TT>
  !     and
  !     <TT>check_bounds_are_exact_dynamic</TT>.
  !     <PRE>
  !       INTERFACE check_bounds_are_exact
  !         MODULE PROCEDURE check_bounds_are_exact_static
  !         MODULE PROCEDURE check_bounds_are_exact_dynamic
  !       END INTERFACE check_bounds_are_exact
  !     </PRE>
  !   </FUTURE>
 
  USE diag_data_mod, ONLY: output_fields, input_fields, files, do_diag_field_log, diag_log_unit,&
       & very_large_axis_length, time_zero, very_large_file_freq, end_of_run, every_time,&
       & diag_seconds, diag_minutes, diag_hours, diag_days, diag_months, diag_years, get_base_time,&
       & time_unit_list, max_files, get_base_year, get_base_month, get_base_day, get_base_hour, get_base_minute,&
       & get_base_second, num_files, max_files, max_fields_per_file, max_out_per_in_field,&
       & max_input_fields,num_input_fields, max_output_fields, num_output_fields, coord_type,&
       & mix_snapshot_average_fields, global_descriptor, cmor_missing_value, use_cmor, pack_size,&
       & debug_diag_manager, flush_nc_files, output_field_type, max_field_attributes, max_file_attributes,&
       & file_type, prepend_date, region_out_use_alt_value, glo_reg_val, glo_reg_val_alt,&
       & diag_field_not_found, diag_init_time
  USE diag_data_mod, ONLY: fileobju, fileobj, fnum_for_domain, fileobjnd
  USE diag_axis_mod, ONLY: get_diag_axis_data, get_axis_global_length, get_diag_axis_cart,&
       & get_domain1d, get_domain2d, diag_subaxes_init, diag_axis_init, get_diag_axis, get_axis_aux,&
       & get_axes_shift, get_diag_axis_name, get_diag_axis_domain_name, get_domainug, &
       & get_axis_reqfld, axis_is_compressed, get_compressed_axes_ids
  USE diag_output_mod, ONLY: diag_output_init, write_axis_meta_data,&
       & write_field_meta_data, done_meta_data, diag_flush
  USE diag_output_mod, ONLY: diag_field_write, diag_write_time !<fms2_io
  USE diag_grid_mod, ONLY: get_local_indexes
  USE fms_diag_time_utils_mod, ONLY: diag_time_inc, get_time_string, get_date_dif
  USE fms_mod, ONLY: error_mesg, fatal, warning, note, mpp_pe, mpp_root_pe, lowercase, fms_error_handler,&
       & string, write_version_number
  USE mpp_domains_mod,ONLY: domain1d, domain2d, mpp_get_compute_domain, null_domain1d, null_domain2d,&
       & OPERATOR(.NE.), OPERATOR(.EQ.), mpp_modify_domain, mpp_get_domain_components,&
       & mpp_get_ntile_count, mpp_get_current_ntile, mpp_get_tile_id, mpp_mosaic_defined, mpp_get_tile_npes,&
       & domainug, null_domainug
  USE time_manager_mod,ONLY: time_type, OPERATOR(==), OPERATOR(>), no_calendar, increment_date,&
       & increment_time, get_calendar_type, get_date, get_time, leap_year, OPERATOR(-),&
       & OPERATOR(<), OPERATOR(>=), OPERATOR(<=), OPERATOR(==)
  USE mpp_mod, ONLY: mpp_npes
  USE constants_mod, ONLY: seconds_per_day, seconds_per_hour, seconds_per_minute
  USE fms2_io_mod
  USE fms_diag_bbox_mod, ONLY: fmsdiagibounds_type
  USE netcdf, ONLY: nf90_char
 
  IMPLICIT NONE
  PRIVATE
  PUBLIC get_subfield_size, log_diag_field_info, update_bounds, check_out_of_bounds,&
       & check_bounds_are_exact_dynamic, check_bounds_are_exact_static, init_file, diag_time_inc,&
       & find_input_field, init_input_field, init_output_field, diag_data_out, write_static,&
       & check_duplicate_output_fields, get_date_dif, get_subfield_vert_size, sync_file_times,&
       & prepend_attribute, attribute_init, diag_util_init,&
       & fms_diag_check_out_of_bounds, &
       & fms_diag_check_bounds_are_exact_dynamic, fms_diag_check_bounds_are_exact_static, get_file_start_time
 
 
  !> @brief Prepend a value to a string attribute in the output field or output file.
  !> @ingroup diag_util_mod
  INTERFACE prepend_attribute
     MODULE PROCEDURE prepend_attribute_field
     MODULE PROCEDURE prepend_attribute_file
  END INTERFACE prepend_attribute
 
  !> @brief Allocates the atttype in out_file.
  !> @ingroup diag_util_mod
  INTERFACE attribute_init
     MODULE PROCEDURE attribute_init_field
     MODULE PROCEDURE attribute_init_file
  END INTERFACE attribute_init
 
  INTERFACE fms_diag_check_out_of_bounds
    module procedure fms_diag_check_out_of_bounds_r4
    module procedure fms_diag_check_out_of_bounds_r8
  END INTERFACE fms_diag_check_out_of_bounds
 
 
!> @addtogroup diag_util_mod
!> @{
 
  ! Include variable "version" to be written to log file.
#include <file_version.h>
 
  LOGICAL :: module_initialized = .false.
 
  character(len=1), public :: field_log_separator = '|' !< separator used for csv-style log of registered fields
                                                !! set by nml in diag_manager init
 
 
CONTAINS
 
  !> @brief Write the version number of this file to the log file.
  SUBROUTINE diag_util_init()
    IF (module_initialized) THEN
       RETURN
    END IF
 
    ! Write version number out to log file
    call write_version_number("DIAG_UTIL_MOD", version)
  END SUBROUTINE diag_util_init
 
  !> @brief Get the size, start, and end indices for output fields.
  SUBROUTINE get_subfield_size(axes, outnum)
    INTEGER, INTENT(in) :: axes(:) !< axes of the input_field
    INTEGER, INTENT(in) :: outnum  !< position in array output_fields
 
    REAL, ALLOCATABLE   :: global_lat(:), global_lon(:), global_depth(:)
    INTEGER :: global_axis_size, global_axis_sizey
    INTEGER :: i,xbegin,xend,ybegin,yend,xbegin_l,xend_l,ybegin_l,yend_l
    CHARACTER(len=1) :: cart
    TYPE(domain2d) :: domain2, domain2_new
    TYPE(domain1d) :: domain1, domain1x, domain1y
    REAL :: start(3) !< start coordinates in 3 axes
    REAL :: end(3) !< end coordinates in 3 axes
    INTEGER :: gstart_indx(3)!< global start indices of output domain in 3 axes
    INTEGER :: gend_indx(3) !< global end indices of output domain in 3 axes
    REAL, ALLOCATABLE :: subaxis_x(:) !< containing local coordinates in x,y,z axes
    REAL, ALLOCATABLE :: subaxis_y(:) !< containing local coordinates in x,y,z axes
    REAL, ALLOCATABLE :: subaxis_z(:) !< containing local coordinates in x,y,z axes
    CHARACTER(len=128) :: msg
    INTEGER :: ishift, jshift
    INTEGER :: grv !< Value used to determine if the region defined in the diag_table is for the whole
                   !! axis, or a sub-axis
    CHARACTER(len=128), DIMENSION(2) :: axis_domain_name
 
    !initilization for local output
    ! initially out of (lat/lon/depth) range
    start = -1.e10
    end = -1.e10
    gstart_indx = -1
    gend_indx = -1
 
    ! get the value to compare to determine if writing full axis data
    IF ( region_out_use_alt_value ) THEN
       grv = glo_reg_val_alt
    ELSE
       grv = glo_reg_val
    END IF
 
    ! get axis data (lat, lon, depth) and indices
    start = output_fields(outnum)%output_grid%start
    end = output_fields(outnum)%output_grid%end
 
    CALL get_diag_axis_domain_name(axes(1), axis_domain_name(1))
    CALL get_diag_axis_domain_name(axes(2), axis_domain_name(2))
 
    IF (   index(lowercase(axis_domain_name(1)), 'cubed') == 0 .AND. &
         & index(lowercase(axis_domain_name(2)), 'cubed') == 0 ) THEN
       DO i = 1, SIZE(axes(:))
          global_axis_size = get_axis_global_length(axes(i))
          output_fields(outnum)%output_grid%subaxes(i) = -1
          CALL get_diag_axis_cart(axes(i), cart)
          SELECT CASE(cart)
          CASE ('X')
             ! <ERROR STATUS="FATAL">wrong order of axes.  X should come first.</ERROR>
             IF( i.NE.1 ) CALL error_mesg('diag_util_mod::get_subfield_size',&
                  & 'wrong order of axes, X should come first',fatal)
             ALLOCATE(global_lon(global_axis_size))
             CALL get_diag_axis_data(axes(i),global_lon)
             IF( int(start(i)) == grv .AND. int(end(i)) == grv ) THEN
                gstart_indx(i) = 1
                gend_indx(i) = global_axis_size
                output_fields(outnum)%output_grid%subaxes(i) = axes(i)
             ELSE
                gstart_indx(i) = get_index(start(i),global_lon)
                gend_indx(i) = get_index(end(i),global_lon)
             END IF
             ALLOCATE(subaxis_x(gstart_indx(i):gend_indx(i)))
             subaxis_x=global_lon(gstart_indx(i):gend_indx(i))
          CASE ('Y')
             ! <ERROR STATUS="FATAL">wrong order of axes, Y should come second.</ERROR>
             IF( i.NE.2 ) CALL error_mesg('diag_util_mod::get_subfield_size',&
                  & 'wrong order of axes, Y should come second',fatal)
             ALLOCATE(global_lat(global_axis_size))
             CALL get_diag_axis_data(axes(i),global_lat)
             IF( int(start(i)) == grv .AND. int(end(i)) == grv ) THEN
                gstart_indx(i) = 1
                gend_indx(i) = global_axis_size
                output_fields(outnum)%output_grid%subaxes(i) = axes(i)
             ELSE
                gstart_indx(i) = get_index(start(i),global_lat)
                gend_indx(i) = get_index(end(i),global_lat)
             END IF
             ALLOCATE(subaxis_y(gstart_indx(i):gend_indx(i)))
             subaxis_y=global_lat(gstart_indx(i):gend_indx(i))
          CASE ('Z')
             ! <ERROR STATUS="FATAL">wrong values in vertical axis of region</ERROR>
             IF ( start(i)*end(i)<0. ) CALL error_mesg('diag_util_mod::get_subfield_size',&
                  & 'wrong values in vertical axis of region',fatal)
             IF ( start(i)>=0. .AND. end(i)>0. ) THEN
                ALLOCATE(global_depth(global_axis_size))
                CALL get_diag_axis_data(axes(i),global_depth)
                gstart_indx(i) = get_index(start(i),global_depth)
                gend_indx(i) = get_index(end(i),global_depth)
                ALLOCATE(subaxis_z(gstart_indx(i):gend_indx(i)))
                subaxis_z=global_depth(gstart_indx(i):gend_indx(i))
                output_fields(outnum)%output_grid%subaxes(i) =&
                     & diag_subaxes_init(axes(i),subaxis_z, gstart_indx(i),gend_indx(i))
                DEALLOCATE(subaxis_z,global_depth)
             ELSE ! regional vertical axis is the same as global vertical axis
                gstart_indx(i) = 1
                gend_indx(i) = global_axis_size
                output_fields(outnum)%output_grid%subaxes(i) = axes(i)
                ! <ERROR STATUS="FATAL">i should equal 3 for z axis</ERROR>
                IF( i /= 3 ) CALL error_mesg('diag_util_mod::get_subfield_size',&
                     & 'i should equal 3 for z axis', fatal)
             END IF
          CASE default
             ! <ERROR STATUS="FATAL">Wrong axis_cart</ERROR>
             CALL error_mesg('diag_util_mod::get_subfield_size', 'Wrong axis_cart', fatal)
          END SELECT
       END DO
 
       DO i = 1, SIZE(axes(:))
          IF ( gstart_indx(i) == -1 .OR. gend_indx(i) == -1 ) THEN
             ! <ERROR STATUS="FATAL">
             !   can not find gstart_indx/gend_indx for <output_fields(outnum)%output_name>,
             !   check region bounds for axis <i>.
             ! </ERROR>
             WRITE(msg,'(A,I2)') ' check region bounds for axis ', i
             CALL error_mesg('diag_util_mod::get_subfield_size', 'can not find gstart_indx/gend_indx for '&
                  & //trim(output_fields(outnum)%output_name)//','//trim(msg), fatal)
          END IF
       END DO
    ELSE ! cubed sphere
       ! get the i and j start and end indexes
       CALL get_local_indexes(lonstart=start(1), lonend=end(1), &
            &                 latstart=start(2), latend=end(2), &
            &                 istart=gstart_indx(1), iend=gend_indx(1), &
            &                 jstart=gstart_indx(2), jend=gend_indx(2))
       global_axis_size =  get_axis_global_length(axes(1))
       ALLOCATE(global_lon(global_axis_size))
       global_axis_sizey = get_axis_global_length(axes(2))
       ALLOCATE(global_lat(global_axis_sizey))
       CALL get_diag_axis_data(axes(1),global_lon)
       CALL get_diag_axis_data(axes(2),global_lat)
 
       !Potential fix for out-of-bounds error for global_lon and global_lat.
       IF ((gstart_indx(1) .GT. 0 .AND. gstart_indx(2) .GT. 0) .AND. &
           (gstart_indx(1) .LE. global_axis_size .AND. gstart_indx(2) .LE. global_axis_sizey) .AND. &
           (gend_indx(1) .GT. 0 .AND. gend_indx(2) .GT. 0) .AND. &
           (gend_indx(1) .LE. global_axis_size .AND. gend_indx(2) .LE. global_axis_sizey)) THEN
          ALLOCATE(subaxis_x(gstart_indx(1):gend_indx(1)))
          ALLOCATE(subaxis_y(gstart_indx(2):gend_indx(2)))
          subaxis_x=global_lon(gstart_indx(1):gend_indx(1))
          subaxis_y=global_lat(gstart_indx(2):gend_indx(2))
       END IF
 
       ! Now deal with the Z component
       IF ( SIZE(axes(:)) > 2 ) THEN
          global_axis_size = get_axis_global_length(axes(3))
          output_fields(outnum)%output_grid%subaxes(3) = -1
          CALL get_diag_axis_cart(axes(3), cart)
          ! <ERROR STATUS="FATAL">
          !   axis(3) should be Z-axis
          ! </ERROR>
          IF ( lowercase(cart) /= 'z' ) CALL error_mesg('diag_util_mod::get_subfield_size', &
               &'axis(3) should be Z-axis', fatal)
          ! <ERROR STATUS="FATAL">
          !   wrong values in vertical axis of region
          ! </ERROR>
          IF ( start(3)*end(3)<0. ) CALL error_mesg('diag_util_mod::get_subfield_size',&
               & 'wrong values in vertical axis of region',fatal)
          IF ( start(3)>=0. .AND. end(3)>0. ) THEN
             ALLOCATE(global_depth(global_axis_size))
             CALL get_diag_axis_data(axes(3),global_depth)
             gstart_indx(3) = get_index(start(3),global_depth)
             IF( start(3) == 0.0 )  gstart_indx(3) = 1
             gend_indx(3) = get_index(end(3),global_depth)
             IF( start(3) >= maxval(global_depth) ) gstart_indx(3)= global_axis_size
             IF( end(3)   >= maxval(global_depth) ) gend_indx(3)  = global_axis_size
 
             ALLOCATE(subaxis_z(gstart_indx(3):gend_indx(3)))
             subaxis_z=global_depth(gstart_indx(3):gend_indx(3))
             output_fields(outnum)%output_grid%subaxes(3) =&
                  & diag_subaxes_init(axes(3),subaxis_z, gstart_indx(3),gend_indx(3))
             DEALLOCATE(subaxis_z,global_depth)
          ELSE ! regional vertical axis is the same as global vertical axis
             gstart_indx(3) = 1
             gend_indx(3) = global_axis_size
             output_fields(outnum)%output_grid%subaxes(3) = axes(3)
          END IF
       END IF
    END IF
 
    ! get domain and compute_domain(xbegin,xend,ybegin,yend)
    xbegin = -1
    xend = -1
    ybegin = -1
    yend = -1
 
    domain2 = get_domain2d(axes)
    IF ( domain2 .NE. null_domain2d ) THEN
       CALL mpp_get_compute_domain(domain2, xbegin, xend, ybegin, yend)
       CALL mpp_get_domain_components(domain2, domain1x, domain1y)
    ELSE
       DO i = 1, min(SIZE(axes(:)),2)
          domain1 = get_domain1d(axes(i))
          IF ( domain1 .NE. null_domain1d ) THEN
             CALL get_diag_axis_cart(axes(i),cart)
             SELECT CASE(cart)
             CASE ('X')
                domain1x = get_domain1d(axes(i))
                CALL mpp_get_compute_domain(domain1x, xbegin, xend)
             CASE ('Y')
                domain1y = get_domain1d(axes(i))
                CALL mpp_get_compute_domain(domain1y, ybegin, yend)
             CASE default ! do nothing here
             END SELECT
          ELSE
             ! <ERROR STATUS="FATAL">No domain available</ERROR>
             CALL error_mesg('diag_util_mod::get_subfield_size', 'NO domain available', fatal)
          END IF
       END DO
    END IF
 
    CALL get_axes_shift(axes, ishift, jshift)
    xend = xend+ishift
    yend = yend+jshift
 
    IF ( xbegin == -1 .OR. xend == -1 .OR. ybegin == -1 .OR. yend == -1 ) THEN
       ! <ERROR STATUS="FATAL">wrong compute domain indices</ERROR>
       CALL error_mesg('diag_util_mod::get_subfield_size', 'wrong compute domain indices',fatal)
    END IF
 
    ! get the area containing BOTH compute domain AND local output area
    IF( gstart_indx(1) > xend .OR. xbegin > gend_indx(1) ) THEN
       output_fields(outnum)%output_grid%l_start_indx(1) = -1
       output_fields(outnum)%output_grid%l_end_indx(1) = -1
       output_fields(outnum)%need_compute = .false. ! not involved
    ELSEIF ( gstart_indx(2) > yend .OR. ybegin > gend_indx(2) ) THEN
       output_fields(outnum)%output_grid%l_start_indx(2) = -1
       output_fields(outnum)%output_grid%l_end_indx(2) = -1
       output_fields(outnum)%need_compute = .false. ! not involved
    ELSE
       output_fields(outnum)%output_grid%l_start_indx(1) = max(xbegin, gstart_indx(1))
       output_fields(outnum)%output_grid%l_start_indx(2) = max(ybegin, gstart_indx(2))
       output_fields(outnum)%output_grid%l_end_indx(1) = min(xend, gend_indx(1))
       output_fields(outnum)%output_grid%l_end_indx(2) = min(yend, gend_indx(2))
       output_fields(outnum)%need_compute = .true.  ! involved in local output
    END IF
 
    IF ( output_fields(outnum)%need_compute ) THEN
       ! need to modify domain1d and domain2d for subaxes
       xbegin_l = output_fields(outnum)%output_grid%l_start_indx(1)
       xend_l = output_fields(outnum)%output_grid%l_end_indx(1)
       ybegin_l = output_fields(outnum)%output_grid%l_start_indx(2)
       yend_l = output_fields(outnum)%output_grid%l_end_indx(2)
       CALL mpp_modify_domain(domain2, domain2_new, xbegin_l,xend_l, ybegin_l,yend_l,&
            & gstart_indx(1),gend_indx(1), gstart_indx(2),gend_indx(2))
 
       output_fields(outnum)%output_grid%subaxes(1) =&
            & diag_subaxes_init(axes(1),subaxis_x, gstart_indx(1),gend_indx(1),domain2_new)
       output_fields(outnum)%output_grid%subaxes(2) =&
            & diag_subaxes_init(axes(2),subaxis_y, gstart_indx(2),gend_indx(2),domain2_new)
       DO i = 1, SIZE(axes(:))
          IF ( output_fields(outnum)%output_grid%subaxes(i) == -1 ) THEN
             ! <ERROR STATUS="FATAL">
             !   <output_fields(outnum)%output_name> error at i = <i>
             ! </ERROR>
             WRITE(msg,'(a,"/",I4)') 'at i = ',i
             CALL error_mesg('diag_util_mod::get_subfield_size '//trim(output_fields(outnum)%output_name),&
                  'error '//trim(msg), fatal)
          END IF
       END DO
 
       ! local start index should start from 1
       output_fields(outnum)%output_grid%l_start_indx(1) = max(xbegin, gstart_indx(1)) - xbegin + 1
       output_fields(outnum)%output_grid%l_start_indx(2) = max(ybegin, gstart_indx(2)) - ybegin + 1
       output_fields(outnum)%output_grid%l_end_indx(1) = min(xend, gend_indx(1)) - xbegin + 1
       output_fields(outnum)%output_grid%l_end_indx(2) = min(yend, gend_indx(2)) - ybegin + 1
       IF ( SIZE(axes(:))>2 ) THEN
          output_fields(outnum)%output_grid%l_start_indx(3) = gstart_indx(3)
          output_fields(outnum)%output_grid%l_end_indx(3) = gend_indx(3)
       ELSE
          output_fields(outnum)%output_grid%l_start_indx(3) = 1
          output_fields(outnum)%output_grid%l_end_indx(3) = 1
       END IF
    END IF
    IF ( ALLOCATED(subaxis_x) ) DEALLOCATE(subaxis_x, global_lon)
    IF ( ALLOCATED(subaxis_y) ) DEALLOCATE(subaxis_y, global_lat)
  END SUBROUTINE get_subfield_size
 
  !> @brief Get size, start and end indices for output fields.
  SUBROUTINE get_subfield_vert_size(axes, outnum)
    INTEGER, DIMENSION(:), INTENT(in) :: axes !< axes of the input_field
    INTEGER, INTENT(in) :: outnum  !< position in array output_fields
 
    REAL, DIMENSION(3) :: start !< start and end coordinates in 3 axes
    REAL, DIMENSION(3) :: end !< start and end coordinates in 3 axes
    REAL, ALLOCATABLE, DIMENSION(:) :: global_depth
    REAL, ALLOCATABLE, DIMENSION(:) :: subaxis_z !< containing local coordinates in x,y,z axes
    INTEGER :: i, global_axis_size
    INTEGER, DIMENSION(3) :: gstart_indx !< global start and end indices of output domain in 3 axes
    INTEGER, DIMENSION(3) :: gend_indx !< global start and end indices of output domain in 3 axes
    CHARACTER(len=1) :: cart
    CHARACTER(len=128) :: msg
!----------
!ug support
    integer :: vert_dim_num
!----------
 
    !initilization for local output
    start = -1.e10
    end = -1.e10 ! initially out of (lat/lon/depth) range
    gstart_indx = -1
    gend_indx=-1
 
    ! get axis data (lat, lon, depth) and indices
    start= output_fields(outnum)%output_grid%start
    end = output_fields(outnum)%output_grid%end
 
!----------
!ug support
    vert_dim_num = 3
!----------
    DO i = 1, SIZE(axes(:))
       global_axis_size = get_axis_global_length(axes(i))
       output_fields(outnum)%output_grid%subaxes(i) = -1
       CALL get_diag_axis_cart(axes(i), cart)
       SELECT CASE(cart)
       CASE ('X')
          ! <ERROR STATUS="FATAL">wrong order of axes, X should come first</ERROR>
          IF ( i.NE.1 ) CALL error_mesg('diag_util_mod::get_subfield_vert_size',&
               & 'wrong order of axes, X should come first',fatal)
          gstart_indx(i) = 1
          gend_indx(i) = global_axis_size
          output_fields(outnum)%output_grid%subaxes(i) = axes(i)
       CASE ('Y')
          ! <ERROR STATUS="FATAL">wrong order of axes, Y should come second</ERROR>
          IF( i.NE.2 ) CALL error_mesg('diag_util_mod::get_subfield_vert_size',&
               & 'wrong order of axes, Y should come second',fatal)
          gstart_indx(i) = 1
          gend_indx(i) = global_axis_size
          output_fields(outnum)%output_grid%subaxes(i) = axes(i)
!----------
!ug support
       case ("U")
           if (i .ne. 1) then
               call error_mesg("diag_util_mod::get_subfield_vert_size", &
                               "the unstructured axis must be the first dimension.", &
                               fatal)
           endif
           gstart_indx(i) = 1
           gend_indx(i) = global_axis_size
           output_fields(outnum)%output_grid%subaxes(i) = axes(i)
           vert_dim_num = 2
           start(vert_dim_num) = start(3)
           end(vert_dim_num) = end(3)
!----------
       CASE ('Z')
!----------
!ug support
          if (i .ne. vert_dim_num) then
              call error_mesg("diag_util_mod::get_subfield_vert_size",&
                              "i should equal vert_dim_num for z axis", &
                              fatal)
          endif
!----------
          ! <ERROR STATUS="FATAL">wrong values in vertical axis of region</ERROR>
          IF( start(i)*end(i) < 0. ) CALL error_mesg('diag_util_mod::get_subfield_vert_size',&
               & 'wrong values in vertical axis of region',fatal)
          IF( start(i) >= 0. .AND. end(i) > 0. ) THEN
             ALLOCATE(global_depth(global_axis_size))
             CALL get_diag_axis_data(axes(i),global_depth)
             gstart_indx(i) = get_index(start(i),global_depth)
             IF( start(i) == 0.0 )  gstart_indx(i) = 1
 
             gend_indx(i) = get_index(end(i),global_depth)
             IF( start(i) >= maxval(global_depth) ) gstart_indx(i)= global_axis_size
             IF( end(i)   >= maxval(global_depth) ) gend_indx(i)  = global_axis_size
 
             ALLOCATE(subaxis_z(gstart_indx(i):gend_indx(i)))
             subaxis_z=global_depth(gstart_indx(i):gend_indx(i))
             output_fields(outnum)%output_grid%subaxes(i) = &
                    diag_subaxes_init(axes(i),subaxis_z, gstart_indx(i),gend_indx(i))
             DEALLOCATE(subaxis_z,global_depth)
          ELSE !   vertical axis is the same as global vertical axis
             gstart_indx(i) = 1
             gend_indx(i) = global_axis_size
             output_fields(outnum)%output_grid%subaxes(i) = axes(i)
          END IF
       CASE default
          ! <ERROR STATUS="FATAL">Wrong axis_cart</ERROR>
          CALL error_mesg('diag_util_mod::get_subfield_vert_size', 'Wrong axis_cart', fatal)
       END SELECT
    END DO
 
    DO i = 1,SIZE(axes(:))
       IF ( gstart_indx(i) == -1 .OR. gend_indx(i) == -1 ) THEN
          ! <ERROR STATUS="FATAL">
          !   can not find gstart_indx/gend_indx for <output_fields(outnum)%output_name>
          !   check region bounds for axis
          ! </ERROR>
          WRITE(msg,'(A,I2)') ' check region bounds for axis ', i
          CALL error_mesg('diag_util_mod::get_subfield_vert_size', 'can not find gstart_indx/gend_indx for '&
               & //trim(output_fields(outnum)%output_name)//','//trim(msg), fatal)
       END IF
    END DO
 
    DO i= 1, 2
       output_fields(outnum)%output_grid%l_start_indx(i) = gstart_indx(i)
       output_fields(outnum)%output_grid%l_end_indx(i)   = gend_indx(i)
    END DO
 
    IF( SIZE(axes(:)) > 2 ) THEN
       output_fields(outnum)%output_grid%l_start_indx(3) = gstart_indx(3)
       output_fields(outnum)%output_grid%l_end_indx(3)   = gend_indx(3)
    ELSE
       output_fields(outnum)%output_grid%l_start_indx(3) = 1
       output_fields(outnum)%output_grid%l_end_indx(3)   = 1
    END IF
  END SUBROUTINE get_subfield_vert_size
 
  !> @brief Find index <TT>i</TT> of array such that <TT>array(i)</TT> is closest to number.
  INTEGER FUNCTION get_index(number, array)
    REAL, INTENT(in) :: number !<
    REAL, INTENT(in), DIMENSION(:) :: array !<
 
    INTEGER :: i, n
    LOGICAL :: found
 
    n = SIZE(array(:))
    ! check if array is monotonous
    DO i = 2, n-1
       IF( (array(i-1)<array(i).AND.array(i)>array(i+1)) .OR. (array(i-1)>array(i).AND.array(i)<array(i+1))) THEN
          ! <ERROR STATUS="FATAL">array NOT monotonously ordered</ERROR>
          CALL error_mesg('diag_util_mod::get_index', 'array NOT monotonously ordered',fatal)
       END IF
    END DO
    get_index = -1
    found = .false.
    ! search in increasing array
    DO i = 1, n-1
       IF ( (array(i)<=number).AND.(array(i+1)>= number) ) THEN
          IF( number - array(i) <= array(i+1) - number ) THEN
             get_index = i
             found=.true.
          ELSE
             get_index = i+1
             found=.true.
          ENDIF
          EXIT
       END IF
    END DO
    ! if not found, search in decreasing array
    IF( .NOT.found ) THEN
       DO i = 1, n-1
          IF ( (array(i)>=number).AND.(array(i+1)<= number) ) THEN
             IF ( array(i)-number <= number-array(i+1) ) THEN
                get_index = i
                found = .true.
             ELSE
                get_index = i+1
                found = .true.
             END IF
             EXIT
          END IF
       END DO
    END IF
    ! if still not found, is it less than the first element
    ! or greater than last element? (Increasing Array)
    ! But it must be within 2x the axis spacing
    ! i.e. array(1)-(array(3)-array(1)).LT.number .AND. or 2*array(1)-array(3).LT.number
    IF ( .NOT. found ) THEN
       IF ( 2*array(1)-array(3).LT.number .AND. number.LT.array(1) ) THEN
          get_index = 1
          found = .true.
       ELSE IF ( array(n).LT.number .AND. number.LT.2*array(n)-array(n-2) ) THEN
          get_index = n
          found = .true.
       ELSE
          found = .false.
       END IF
    END IF
 
    ! if still not found, is it greater than the first element
    ! or less than the last element? (Decreasing Array)
    ! But it must be within 2x the axis spacing (see above)
    IF ( .NOT. found ) THEN
       IF ( 2*array(1)-array(3).GT.number .AND. number.GT.array(1) ) THEN
          get_index = 1
          found = .true.
       ELSE IF ( array(n).GT.number .AND. number.GT.2*array(n)-array(n-2) ) THEN
          get_index = n
          found = .true.
       ELSE
          found = .false.
       END IF
    END IF
  END FUNCTION get_index
 
  !> @brief Writes brief diagnostic field info to the log file.
  !! @details If the <TT>do_diag_field_log</TT> namelist parameter is .TRUE.,
  !!     then a line briefly describing diagnostic field is added to
  !!     the log file.  Normally users should not call this subroutine
  !!     directly, since it is called by register_static_field and
  !!     register_diag_field if do_not_log is not set to .TRUE..  It is
  !!     used, however, in LM3 to avoid excessive logs due to the
  !!     number of fields registered for each of the tile types.  LM3
  !!     code uses a do_not_log parameter in the registration calls,
  !!     and subsequently calls this subroutine to log field information
  !!     under a generic name.
  SUBROUTINE log_diag_field_info(module_name, field_name, axes, long_name, units,&
                              & missing_value, range, dynamic )
    CHARACTER(len=*), INTENT(in) :: module_name !< Module name
    CHARACTER(len=*), INTENT(in) :: field_name !< Field name
    INTEGER, DIMENSION(:), INTENT(in) :: axes !< Axis IDs
    CHARACTER(len=*), OPTIONAL, INTENT(in) :: long_name !< Long name for field.
    CHARACTER(len=*), OPTIONAL, INTENT(in) :: units !< Unit of field.
    CLASS(*), OPTIONAL, INTENT(in) :: missing_value !< Missing value value.
    CLASS(*), DIMENSION(:), OPTIONAL, INTENT(IN) :: range !< Valid range of values for field.
    LOGICAL, OPTIONAL, INTENT(in) :: dynamic !< <TT>.TRUE.</TT> if field is not static.
 
    ! ---- local vars
    CHARACTER(len=256) :: lmodule, lfield, lname, lunits
    CHARACTER(len=64)  :: lmissval, lmin, lmax
    CHARACTER(len=8)   :: numaxis, timeaxis
    CHARACTER(len=1)   :: sep = '|'
    CHARACTER(len=256) :: axis_name, axes_list
    INTEGER :: i
    REAL :: missing_value_use !< Local copy of missing_value
    REAL, DIMENSION(2) :: range_use !< Local copy of range
 
    IF ( .NOT.do_diag_field_log ) RETURN
    IF ( mpp_pe().NE.mpp_root_pe() ) RETURN
 
    ! Fatal error if range is present and its extent is not 2.
    IF ( PRESENT(range) ) THEN
      IF ( SIZE(range) .NE. 2 ) THEN
        CALL error_mesg('diag_util_mod::fms_log_field_info', 'extent of range should be 2', fatal)
      END IF
    END IF
 
    lmodule = trim(module_name)
    lfield = trim(field_name)
 
    IF ( PRESENT(long_name) ) THEN
      lname  = trim(long_name)
    ELSE
      lname  = ''
    END IF
 
    IF ( PRESENT(units) ) THEN
      lunits = trim(units)
    ELSE
      lunits = ''
    END IF
 
    WRITE (numaxis,'(i1)') SIZE(axes)
 
    IF (PRESENT(missing_value)) THEN
      IF ( use_cmor ) THEN
        WRITE (lmissval,*) cmor_missing_value
      ELSE
        SELECT TYPE (missing_value)
          TYPE IS (real(kind=r4_kind))
            missing_value_use = missing_value
          TYPE IS (real(kind=r8_kind))
            missing_value_use = real(missing_value)
          CLASS DEFAULT
            CALL error_mesg ('diag_util_mod::log_diag_field_info',&
              & 'The missing_value is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
        END SELECT
        WRITE (lmissval,*) missing_value_use
      END IF
    ELSE
      lmissval = ''
    ENDIF
 
    IF ( PRESENT(range) ) THEN
      SELECT TYPE (range)
      TYPE IS (real(kind=r4_kind))
          range_use = range
      TYPE IS (real(kind=r8_kind))
          range_use = real(range)
      CLASS DEFAULT
          CALL error_mesg ('diag_util_mod::log_diag_field_info',&
              & 'The range is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
      END SELECT
      WRITE (lmin,*) range_use(1)
      WRITE (lmax,*) range_use(2)
    ELSE
      lmin = ''
      lmax = ''
    END IF
 
    IF ( PRESENT(dynamic) ) THEN
      IF (dynamic) THEN
          timeaxis = 'T'
      ELSE
          timeaxis = 'F'
      END IF
    ELSE
      timeaxis = ''
    END IF
 
    axes_list=''
    DO i = 1, SIZE(axes)
       CALL get_diag_axis_name(axes(i),axis_name)
       IF ( trim(axes_list) /= '' ) axes_list = trim(axes_list)//','
       axes_list = trim(axes_list)//trim(axis_name)
    END DO
 
    WRITE (diag_log_unit,'(777a)') &
        & trim(lmodule),  field_log_separator, trim(lfield),  field_log_separator, trim(lname),    field_log_separator,&
        & trim(lunits),   field_log_separator, trim(numaxis), field_log_separator, trim(timeaxis), field_log_separator,&
        & trim(lmissval), field_log_separator, trim(lmin),    field_log_separator, trim(lmax),     field_log_separator,&
        & trim(axes_list)
  END SUBROUTINE log_diag_field_info
 
 
 
  !> @brief Update the <TT>output_fields</TT> x, y, and z min and max boundaries (array indices)
  !! with the six specified bounds values.
  SUBROUTINE update_bounds(out_num, lower_i, upper_i, lower_j, upper_j, lower_k, upper_k)
    INTEGER, INTENT(in) :: out_num !< output field ID
    INTEGER, INTENT(in) :: lower_i !< Lower i bound.
    INTEGER, INTENT(in) :: upper_i !< Upper i bound.
    INTEGER, INTENT(in) :: lower_j !< Lower j bound.
    INTEGER, INTENT(in) :: upper_j !< Upper j bound.
    INTEGER, INTENT(in) :: lower_k !< Lower k bound.
    INTEGER, INTENT(in) :: upper_k !< Upper k bound.
    CALL output_fields(out_num)%buff_bounds%update_bounds &
      & ( lower_i, upper_i, lower_j, upper_j, lower_k, upper_k )
  END SUBROUTINE update_bounds
 
 
 
  !> @brief Compares the bounding indices of an array specified in "current_bounds"
!! to the corresponding lower and upper bounds specified in "bounds"
!! Comparison is done by the two user specified input functions lowerb_comp and upperb_comp.
!! If any compariosn function returns true, then, after filling error_str, this routine also returns
!! true. The suplied comparison functions should return true for errors : for indices out of bounds,
!! or indices are not equal when expected to be equal.
LOGICAL FUNCTION compare_buffer_bounds_to_size(current_bounds, bounds, error_str, lowerb_comp, upperb_comp)
  TYPE (fmsdiagibounds_type), INTENT(in) :: current_bounds  !<A bounding box holding the current bounds
                                                                     !! of an array.
  TYPE (fmsdiagibounds_type), INTENT(in):: bounds !< The bounding box to check against.
  CHARACTER(*), INTENT(out) :: error_str !< The return status, which is set to non-empty message
                                             !! if the check fails.
 
  !> @brief Interface lowerb_comp should be used for comparison to lower bounds of buffer.
  INTERFACE
     LOGICAL FUNCTION lowerb_comp(a , b)
       INTEGER, INTENT(IN) :: a !< One of the two args that are to be compared to each other.
       INTEGER, INTENT(IN) :: b !< One of the two args that are to be compared to each other.
     END FUNCTION lowerb_comp
  END INTERFACE
 
  !> @brief Interface lowerb_comp should be used for comparison to upper bounds of buffer.
  INTERFACE
     LOGICAL FUNCTION upperb_comp(a, b)
       INTEGER, INTENT(IN) :: a !< One of the two args that are to be compared to each other.
       INTEGER, INTENT(IN) :: b !< One of the two args that are to be compared to each other.
     END FUNCTION upperb_comp
  END INTERFACE
 
  compare_buffer_bounds_to_size = .false.
 
  IF (lowerb_comp( bounds%get_imin() , current_bounds%get_imin()) .OR. &
       upperb_comp( bounds%get_imax() , current_bounds%get_imax()).OR.&
       lowerb_comp( bounds%get_jmin() , current_bounds%get_jmin()) .OR.&
       upperb_comp( bounds%get_jmax() , current_bounds%get_jmax()) .OR.&
       lowerb_comp( bounds%get_kmin() , current_bounds%get_kmin()) .OR.&
       upperb_comp( bounds%get_kmax() , current_bounds%get_kmax())) THEN
    compare_buffer_bounds_to_size = .true.
     error_str ='Buffer bounds=   :   ,   :   ,   :     Actual bounds=   :   ,   :   ,   :   '
     WRITE(error_str(15:17),'(i3)') current_bounds%get_imin()
     WRITE(error_str(19:21),'(i3)') current_bounds%get_imax()
     WRITE(error_str(23:25),'(i3)') current_bounds%get_jmin()
     WRITE(error_str(27:29),'(i3)') current_bounds%get_jmax()
     WRITE(error_str(31:33),'(i3)') current_bounds%get_kmin()
     WRITE(error_str(35:37),'(i3)') current_bounds%get_kmax()
     WRITE(error_str(54:56),'(i3)') bounds%get_imin()
     WRITE(error_str(58:60),'(i3)') bounds%get_imax()
     WRITE(error_str(62:64),'(i3)') bounds%get_jmin()
     WRITE(error_str(66:68),'(i3)') bounds%get_jmax()
     WRITE(error_str(70:72),'(i3)') bounds%get_kmin()
     WRITE(error_str(74:76),'(i3)') bounds%get_kmax()
  ELSE
    compare_buffer_bounds_to_size = .false.
    error_str = ''
  END IF
END FUNCTION compare_buffer_bounds_to_size
 
!> @brief return true iff a<b.
LOGICAL FUNCTION a_lessthan_b(a , b)
 INTEGER, INTENT(IN) :: a !< The first of the two integer args that are to be compared to each other.
 INTEGER, INTENT(IN) :: b !< The first of the two integer args that are to be compared to each other.
 a_lessthan_b = a < b
END FUNCTION a_lessthan_b
 
!> @brief return true iff a>b.
LOGICAL FUNCTION a_greaterthan_b(a, b)
 INTEGER, INTENT(IN) :: a !< The first of the two integer args that are to be compared to each other.
 INTEGER, INTENT(IN) :: b !< The first of the two integer args that are to be compared to each other.
 a_greaterthan_b = a > b
END FUNCTION a_greaterthan_b
 
!> @brief return true iff a /= b
LOGICAL FUNCTION a_noteq_b(a, b)
INTEGER, INTENT(IN) :: a !< The first of the two integer args that are to be compared to each other.
INTEGER, INTENT(IN) :: b !< The first of the two integer args that are to be compared to each other.
a_noteq_b = a /= b
END FUNCTION a_noteq_b
 
  !> @brief Checks if the array indices for <TT>output_fields(out_num)</TT> are outside the
  !! <TT>output_fields(out_num)%buffer</TT> upper and lower bounds.
  !! If there is an error then error message will be filled.
SUBROUTINE check_out_of_bounds(out_num, diag_field_id, err_msg)
  INTEGER, INTENT(in) :: out_num !< Output field ID number.
  INTEGER, INTENT(in) :: diag_field_id !< Input field ID number.
  CHARACTER(len=*), INTENT(out) :: err_msg !< Return status of <TT>check_out_of_bounds</TT>.  An empty
                                           !! error string indicates the x, y, and z indices are not outside the
 
  CHARACTER(len=128) :: error_string1, error_string2
  LOGICAL :: out_of_bounds = .true.
  TYPE (fmsdiagibounds_type) :: array_bounds
  associate(buff_bounds => output_fields(out_num)%buff_bounds)
 
    CALL array_bounds%reset_bounds_from_array_4D(output_fields(out_num)%buffer)
 
    out_of_bounds = compare_buffer_bounds_to_size(array_bounds, buff_bounds, &
     & error_string2, a_lessthan_b, a_greaterthan_b)
 
    IF (out_of_bounds .EQV. .true.) THEN
      WRITE(error_string1,'(a,"/",a)')  trim(input_fields(diag_field_id)%module_name), &
        & trim(output_fields(out_num)%output_name)
      err_msg = 'module/output_field='//trim(error_string1)//&
          & '  Bounds of buffer exceeded.  '//trim(error_string2)
      !   imax, imin, etc need to be reset in case the program is not terminated.
      call buff_bounds%reset(very_large_axis_length, 0)
    ELSE
      err_msg = ''
    END IF
  end associate
END SUBROUTINE check_out_of_bounds
 
 !> @brief Checks if the array indices for <TT>output_fields(out_num)</TT> are outside the
  !! <TT>output_fields(out_num)%buffer</TT> upper and lower bounds.
  !! If there is an error then error message will be filled.
SUBROUTINE fms_diag_check_out_of_bounds_r4(ofb, bounds, output_name, module_name, err_msg)
  REAL(kind=r4_kind), INTENT (in), DIMENSION(:,:,:,:,:) :: ofb !< The output field buffer to check
  TYPE (fmsDiagIbounds_type), INTENT(inout) :: bounds !< The bounding box to check against
  CHARACTER(:), ALLOCATABLE, INTENT(in) :: output_name !< output name for placing in error message
  CHARACTER(:), ALLOCATABLE, INTENT(in) :: module_name !< module name for placing in error message
  CHARACTER(len=*), INTENT(inout) :: err_msg !< Return status of <TT>check_out_of_bounds</TT>.  An empty
                                           !! error string indicates the x, y, and z indices are not outside the
 
  CHARACTER(len=128) :: error_string1, error_string2
  LOGICAL :: out_of_bounds = .true.
  TYPE (fmsDiagIbounds_type) :: array_bounds
 
  CALL array_bounds%reset_bounds_from_array_5D(ofb)
 
  out_of_bounds = compare_buffer_bounds_to_size(array_bounds, bounds, &
     & error_string2, a_lessthan_b, a_greaterthan_b)
 
  IF (out_of_bounds .EQV. .true.) THEN
     WRITE(error_string1,'(a,"/",a)')  trim(module_name), trim(output_name)
     err_msg = 'module/output_field='//trim(error_string1)//&
          & '  Bounds of buffer exceeded.  '//trim(error_string2)
     !   imax, imin, etc need to be reset in case the program is not terminated.
     call bounds%reset(very_large_axis_length,0)
  ELSE
     err_msg = ''
  END IF
END SUBROUTINE fms_diag_check_out_of_bounds_r4
 
 !> @brief Checks if the array indices for output_field buffer (ofb) are outside the
  !! are outside the bounding box (bounds).
  !! If there is an error then error message will be filled.
 
SUBROUTINE fms_diag_check_out_of_bounds_r8(ofb, bounds, output_name, module_name, err_msg)
  REAL(kind=r8_kind), INTENT (in), DIMENSION(:,:,:,:,:) :: ofb !< The output field buffer to check
  TYPE (fmsDiagIbounds_type), INTENT(inout) :: bounds !< The bounding box to check against
  CHARACTER(:), ALLOCATABLE, INTENT(in) :: output_name  !< output name for placing in error message
  CHARACTER(:), ALLOCATABLE, INTENT(in) :: module_name   !< module name for placing in error message
  CHARACTER(len=*), INTENT(out) :: err_msg !< Return status of <TT>check_out_of_bounds</TT>.  An empty
                                           !! error string indicates the x, y, and z indices are not outside the
 
  CHARACTER(len=128) :: error_string1, error_string2
  LOGICAL :: out_of_bounds = .true.
  TYPE (fmsDiagIbounds_type) :: array_bounds  !<A bounding box holdstore the current bounds ofb
 
  CALL array_bounds%reset_bounds_from_array_5D(ofb)
 
  out_of_bounds = compare_buffer_bounds_to_size(array_bounds, bounds, &
     &  error_string2, a_lessthan_b, a_greaterthan_b)
 
  IF (out_of_bounds .EQV. .true.) THEN
     WRITE(error_string1,'(a,"/",a)')  trim(module_name), trim(output_name)
     err_msg = 'module/output_field='//trim(error_string1)//&
          & '  Bounds of buffer exceeded.  '//trim(error_string2)
     !   imax, imin, etc need to be reset in case the program is not terminated.
     call bounds%reset(very_large_axis_length,0)
  ELSE
     err_msg = ''
  END IF
END SUBROUTINE fms_diag_check_out_of_bounds_r8
 
 
!> @brief Checks that array indices specified in the bounding box "current_bounds"
!! are identical to those in the bounding box "bounds" match exactly. The check
!! occurs only when the time changed.
!! If there is an error then error message will be filled.
SUBROUTINE fms_diag_check_bounds_are_exact_dynamic(current_bounds, bounds, output_name, module_name, &
     &  Time, field_prev_Time, err_msg)
  TYPE (fmsdiagibounds_type), INTENT(in) :: current_bounds !<A bounding box holding the current bounds
                                                           !! of an array.
  TYPE (fmsdiagibounds_type), INTENT(inout) :: bounds !<The bounding box to check against
  CHARACTER(:), ALLOCATABLE, INTENT(in) :: output_name  !< output name for placing in error message
  CHARACTER(:), ALLOCATABLE, INTENT(in) :: module_name  !< module name for placing in error message
  TYPE(time_type), INTENT(in)  ::  time !< Time to use in check.  The check is only performed if
  !! <TT>output_fields(out_num)%Time_of_prev_field_data</TT> is not
  !! equal to <TT>Time</TT> or <TT>Time_zero</TT>.
  TYPE(time_type), INTENT(inout) :: field_prev_time  !< <TT>output_fields(out_num)%Time_of_prev_field_data</TT>
  CHARACTER(len=*), INTENT(out) :: err_msg !< Return status of <TT>check_bounds_are_exact_dynamic</TT>.
  !! An empty error string indicates the x, y, and z indices are
  !!     equal to the buffer array boundaries.
 
  CHARACTER(len=128) :: error_string1, error_string2
  LOGICAL :: do_check
  LOGICAL :: lims_not_exact
 
  err_msg = ''
 
  ! Check bounds only when the value of Time changes. When windows are used,
  ! a change in Time indicates that a new loop through the windows has begun,
  !  so a check of the previous loop can be done.
  IF ( time == field_prev_time ) THEN
     do_check = .false.
  ELSE
     IF ( field_prev_time == time_zero ) THEN
        ! It may or may not be OK to check, I don't know how to tell.
        ! Check will be done on subsequent calls anyway.
        do_check = .false.
     ELSE
        do_check = .true.
     END IF
     field_prev_time = time
  END IF
 
  IF ( do_check ) THEN
     lims_not_exact = compare_buffer_bounds_to_size(current_bounds, bounds, &
          & error_string2, a_noteq_b, a_noteq_b)
     IF( lims_not_exact .eqv. .true.) THEN
        WRITE(error_string1,'(a,"/",a)') trim(module_name), trim(output_name)
        err_msg = trim(error_string1)//' Bounds of data do not match those of buffer. '//trim(error_string2)
     END IF
     call bounds%reset(very_large_axis_length, 0)
  END IF
END SUBROUTINE fms_diag_check_bounds_are_exact_dynamic
 
 
!> @brief This is an adaptor to the check_bounds_are_exact_dynamic_modern function to
!! maintain an interface servicing the legacy diag_manager.
SUBROUTINE check_bounds_are_exact_dynamic(out_num, diag_field_id, Time, err_msg)
  INTEGER, INTENT(in) :: out_num !< Output field ID number.
  INTEGER, INTENT(in) :: diag_field_id !< Input field ID number.
  TYPE(time_type), INTENT(in) :: time !< Time to use in check.  The check is only performed if
                                      !! <TT>output_fields(out_num)%Time_of_prev_field_data</TT> is not
                                      !! equal to <TT>Time</TT> or <TT>Time_zero</TT>.
  CHARACTER(len=*), INTENT(out) :: err_msg !< Return status of <TT>check_bounds_are_exact_dynamic</TT>.
                                           !! An empty error string indicates the x, y, and z indices are
                                           !! equal to the buffer array boundaries.
  CHARACTER(:), ALLOCATABLE :: output_name  !< output name for placing in error message
  CHARACTER(:), ALLOCATABLE :: module_name  !< module name for placing in error message
  TYPE (fmsdiagibounds_type) :: current_bounds  !< a bounding box to store the current bounds of the array.
 
  output_name = output_fields(out_num)%output_name
  module_name = input_fields(diag_field_id)%module_name
 
  CALL current_bounds%reset_bounds_from_array_4D(output_fields(out_num)%buffer)
 
  CALL fms_diag_check_bounds_are_exact_dynamic(current_bounds, output_fields(out_num)%buff_bounds, &
       &  output_name, module_name, &
       & time, output_fields(out_num)%Time_of_prev_field_data, err_msg)
 
END SUBROUTINE check_bounds_are_exact_dynamic
 
 
  !> @brief Check if the array indices for <TT>output_fields(out_num)</TT> are equal to the
  !! <TT>output_fields(out_num)%buffer</TT> upper and lower bounds.
  SUBROUTINE check_bounds_are_exact_static(out_num, diag_field_id, err_msg)
    INTEGER, INTENT(in) :: out_num !< Output field ID
    INTEGER, INTENT(in) :: diag_field_id !< Input field ID.
    CHARACTER(len=*), INTENT(out) :: err_msg !< The return status, which is set to non-empty message
                                                  !! if the check fails.
    CHARACTER(:), ALLOCATABLE :: output_name !< output name for placing in error message
    CHARACTER(:), ALLOCATABLE :: module_name !< output name for placing in error message
    TYPE (fmsdiagibounds_type) :: current_bounds !< a bounding box to store the current bounds of the array.
 
    output_name = output_fields(out_num)%output_name
    module_name = input_fields(diag_field_id)%module_name
 
    CALL current_bounds%reset_bounds_from_array_4D(output_fields(out_num)%buffer)
 
    CALL fms_diag_check_bounds_are_exact_static(current_bounds, output_fields(out_num)%buff_bounds, &
       &  output_name, module_name, err_msg)
  END SUBROUTINE check_bounds_are_exact_static
 
 
  !> @brief Check if the array indices specified in the bounding box "current_bounds" are equal to those
  !! specified in the bounding box "bounds" output_fields are equal to the buffer upper and lower bounds.
    !! If there is an error then error message will be filled.
  SUBROUTINE fms_diag_check_bounds_are_exact_static(current_bounds, bounds, output_name, module_name, err_msg)
    TYPE (fmsdiagibounds_type), INTENT(in) :: current_bounds  !<A bounding box holding the current bounds
                                                             !! of the array.
    TYPE (fmsdiagibounds_type), INTENT(inout) :: bounds  !<The bounding box to check against
    CHARACTER(:), ALLOCATABLE, INTENT(in) :: output_name !< output name for placing in error message
    CHARACTER(:), ALLOCATABLE, INTENT(in) :: module_name  !< module name for placing in error message
    CHARACTER(len=*), INTENT(out) :: err_msg !< The return status, which is set to non-empty message
                                               !! if the check fails.
 
    CHARACTER(len=128)  :: error_string1, error_string2
    LOGICAL :: lims_not_exact
 
    err_msg = ''
    lims_not_exact = compare_buffer_bounds_to_size(current_bounds, bounds,  &
         & error_string2, a_noteq_b, a_noteq_b)
    IF( lims_not_exact .eqv. .true.) THEN
       WRITE(error_string1,'(a,"/",a)') trim(module_name), trim(output_name)
       err_msg = trim(error_string1)//' Bounds of data do not match those of buffer. '//trim(error_string2)
    END IF
  call bounds%reset(very_large_axis_length, 0)
  END SUBROUTINE fms_diag_check_bounds_are_exact_static
 
 
  !> @brief Initialize the output file.
  SUBROUTINE init_file(name, output_freq, output_units, format, time_units, long_name, tile_count,&
       & new_file_freq, new_file_freq_units, start_time, file_duration, file_duration_units, filename_time_bounds)
    CHARACTER(len=*), INTENT(in) :: name !< File name.
    CHARACTER(len=*), INTENT(in) :: long_name !< Long name for time axis.
    INTEGER, INTENT(in) :: output_freq !< How often data is to be written to the file.
    INTEGER, INTENT(in) :: output_units !< The output frequency unit.  (MIN, HOURS, DAYS, etc.)
    INTEGER, INTENT(in) :: format !< Number type/kind the data is to be written out to the file.
    INTEGER, INTENT(in) :: time_units !< Time axis units.
    INTEGER, INTENT(in) :: tile_count !< Tile number.
    INTEGER, INTENT(in), OPTIONAL :: new_file_freq !< How often a new file is to be created.
    INTEGER, INTENT(in), OPTIONAL :: new_file_freq_units !< The new file frequency unit.  (MIN, HOURS, DAYS, etc.)</IN>
    INTEGER, INTENT(in), OPTIONAL :: file_duration !< How long file is to be used.
    INTEGER, INTENT(in), OPTIONAL :: file_duration_units !< File duration unit.  (MIN, HOURS, DAYS, etc.)
    TYPE(time_type), INTENT(in), OPTIONAL :: start_time !< Time when the file is to start
    CHARACTER(len=*), INTENT(in), OPTIONAL :: filename_time_bounds !< Set time bound file name to
                                                                   ! use "begin" "middle" or "end"
    INTEGER :: new_file_freq1, new_file_freq_units1
    INTEGER :: file_duration1, file_duration_units1
    INTEGER :: n
    LOGICAL :: same_file_err !< .FALSE. indicates that if the file name had
                             !! previously been registered, this new file
                             !! contained differences from the previous.
    REAL, DIMENSION(1) :: tdata
    CHARACTER(len=128) :: time_units_str
 
    ! Check if this file has already been defined
    same_file_err=.false. ! To indicate that if this file was previously defined
                          ! no differences in this registration was detected.
    DO n=1,num_files
      IF ( trim(files(n)%name) == trim(name) ) THEN
        ! File is defined, check if all inputs are the same
        ! Start with the required parameters
        IF ( files(n)%output_freq.NE.output_freq .OR.&
           & files(n)%output_units.NE.output_units .OR.&
           & files(n)%format.NE.format .OR.&
           & files(n)%time_units.NE.time_units .OR.&
           & trim(files(n)%long_name).NE.trim(long_name) .OR.&
           & files(n)%tile_count.NE.tile_count ) THEN
           same_file_err=.true.
        END IF
 
        ! Now check the OPTIONAL parameters
        IF ( PRESENT(new_file_freq) ) THEN
           IF ( files(n)%new_file_freq.NE.new_file_freq ) THEN
             same_file_err=.true.
           END IF
        END IF
 
        IF ( PRESENT(new_file_freq_units) ) THEN
          IF ( files(n)%new_file_freq_units.NE.new_file_freq_units ) THEN
            same_file_err=.true.
          END IF
        END IF
 
        IF ( PRESENT(start_time) ) THEN
          IF ( files(n)%start_time==start_time ) THEN
            same_file_err=.true.
          END IF
        END IF
 
        IF ( PRESENT(file_duration) ) THEN
          IF ( files(n)%duration.NE.file_duration) THEN
            same_file_err=.true.
          END IF
        END IF
 
        IF ( PRESENT(file_duration_units) ) THEN
          IF ( files(n)%duration_units.NE.file_duration_units ) THEN
            same_file_err=.true.
          END IF
        END IF
 
        ! If the same file was defined twice, simply return, else FATAL
        IF ( same_file_err ) THEN
          ! Something in this file is not identical to the previously defined
          ! file of the same name.  FATAL
          CALL error_mesg('diag_util_mod::init_file',&
                  & 'The file "'//trim(name)//'" is defined multiple times in&
                  & the diag_table.', fatal)
        ELSE
          ! Issue a note that the same file is defined multiple times
          CALL error_mesg('diag_util_mod::init_file',&
                  & 'The file "'//trim(name)//'" is defined multiple times in&
                  & the diag_table.', note)
          ! Return to the calling function
          RETURN
        END IF
      END IF
    END DO
 
    ! Get a number for this file
    num_files = num_files + 1
    IF ( num_files >= max_files ) THEN
       ! <ERROR STATUS="FATAL">
       !   max_files exceeded, increase max_files via the max_files variable
       !   in the namelist diag_manager_nml.
       ! </ERROR>
       CALL error_mesg('diag_util_mod::init_file',&
            & ' max_files exceeded, increase max_files via the max_files variable&
            & in the namelist diag_manager_nml.', fatal)
    END IF
 
    IF ( PRESENT(new_file_freq) ) THEN
       new_file_freq1 = new_file_freq
    ELSE
       new_file_freq1 = very_large_file_freq
    END IF
 
    IF ( PRESENT(new_file_freq_units) ) THEN
       new_file_freq_units1 = new_file_freq_units
    ELSE IF ( get_calendar_type() == no_calendar ) THEN
       new_file_freq_units1 = diag_days
    ELSE
       new_file_freq_units1 = diag_years
    END IF
 
    IF ( PRESENT(file_duration) ) THEN
       file_duration1 = file_duration
    ELSE
       file_duration1 = new_file_freq1
    END IF
 
    IF ( PRESENT(file_duration_units) ) THEN
       file_duration_units1 = file_duration_units
    ELSE
       file_duration_units1 = new_file_freq_units1
    END IF
 
    files(num_files)%tile_count = tile_count
    files(num_files)%name = trim(name)
    files(num_files)%output_freq = output_freq
    files(num_files)%output_units = output_units
    files(num_files)%format = FORMAT
    files(num_files)%time_units = time_units
    files(num_files)%long_name = trim(long_name)
    files(num_files)%num_fields = 0
    files(num_files)%local = .false.
    files(num_files)%last_flush = get_base_time()
    files(num_files)%file_unit = -1
    files(num_files)%new_file_freq = new_file_freq1
    files(num_files)%new_file_freq_units = new_file_freq_units1
    files(num_files)%duration = file_duration1
    files(num_files)%duration_units = file_duration_units1
!> Initialize the times to 0
    files(num_files)%rtime_current = -1.0
    files(num_files)%time_index = 0
    files(num_files)%filename_time_bounds = filename_time_bounds
 
    IF ( PRESENT(start_time) ) THEN
       files(num_files)%start_time = start_time
    ELSE
       files(num_files)%start_time = get_base_time()
    END IF
    files(num_files)%next_open=diag_time_inc(files(num_files)%start_time,new_file_freq1,new_file_freq_units1)
    files(num_files)%close_time = diag_time_inc(files(num_files)%start_time,file_duration1, file_duration_units1)
    IF ( files(num_files)%close_time>files(num_files)%next_open ) THEN
       ! <ERROR STATUS="FATAL">
       !   close time GREATER than next_open time, check file duration,
       !   file frequency in <files(num_files)%name>
       ! </ERROR>
       CALL error_mesg('diag_util_mod::init_file', 'close time GREATER than next_open time, check file duration,&
            & file frequency in '//files(num_files)%name, fatal)
    END IF
 
    ! add time_axis_id and time_bounds_id here
    WRITE(time_units_str, 11) trim(time_unit_list(files(num_files)%time_units)), get_base_year(),&
         & get_base_month(), get_base_day(), get_base_hour(), get_base_minute(), get_base_second()
11  FORMAT(a, ' since ', i4.4, '-', i2.2, '-', i2.2, ' ', i2.2, ':', i2.2, ':', i2.2)
    files(num_files)%time_axis_id = diag_axis_init(trim(long_name), tdata, time_units_str, 'T',&
         & trim(long_name) , set_name=trim(name) )
    !---- register axis for storing time boundaries
    files(num_files)%time_bounds_id = diag_axis_init( 'nv',(/1.,2./),'none','N','vertex number',&
         & set_name='nv')
  END SUBROUTINE init_file
 
  !> @brief Synchronize the file's start and close times with the model start and end times.
  !! @details <TT>sync_file_times</TT> checks to see if the file start time is less than the
  !!     model's init time (passed in as the only argument).  If it is less, then the
  !!     both the file start time and end time are synchronized using the passed in initial time
  !!     and the duration as calculated by the <TT>diag_time_inc</TT> function.  <TT>sync_file_times</TT>
  !!     will also increase the <TT>next_open</TT> until it is greater than the init_time.
  SUBROUTINE sync_file_times(file_id, init_time, err_msg)
    INTEGER, INTENT(in) :: file_id !< The file ID
    TYPE(time_type), INTENT(in) :: init_time !< Initial time use for the synchronization.
    CHARACTER(len=*), OPTIONAL, INTENT(out) :: err_msg !< Return error message
 
    CHARACTER(len=128) :: msg
 
    IF ( PRESENT(err_msg) ) err_msg = ''
 
    IF ( files(file_id)%start_time < init_time ) THEN
       ! Sync the start_time of the file with the initial time of the model
       files(file_id)%start_time = init_time
       ! Sync the file's close time also
       files(file_id)%close_time = diag_time_inc(files(file_id)%start_time,&
            & files(file_id)%duration, files(file_id)%duration_units)
    END IF
 
    ! Need to increase next_open until it is greate than init_time
    DO WHILE ( files(file_id)%next_open <= init_time )
       files(file_id)%next_open = diag_time_inc(files(file_id)%next_open,&
            & files(file_id)%new_file_freq, files(file_id)%new_file_freq_units, err_msg=msg)
       IF ( msg /= '' ) THEN
          IF ( fms_error_handler('diag_util_mod::sync_file_times',&
               & ' file='//trim(files(file_id)%name)//': '//trim(msg), err_msg) ) RETURN
       END IF
    END DO
  END SUBROUTINE sync_file_times
 
  !> @brief Return the file number for file name and tile.
  !! @return Integer find_file
  INTEGER FUNCTION find_file(name, tile_count)
    INTEGER, INTENT(in) :: tile_count !< Tile number.
    CHARACTER(len=*), INTENT(in) :: name !< File name.
 
    INTEGER :: i
 
    find_file = -1
    DO i = 1, num_files
       IF( trim(files(i)%name) == trim(name) .AND. tile_count == files(i)%tile_count ) THEN
          find_file = i
          RETURN
       END IF
    END DO
  END FUNCTION find_file
 
  !> @brief Return the field number for the given module name, field name, and tile number.
  !! @return Integer find_input_field
  INTEGER FUNCTION find_input_field(module_name, field_name, tile_count)
    CHARACTER(len=*), INTENT(in) :: module_name !< Module name.
    CHARACTER(len=*), INTENT(in) :: field_name !< field name.
    INTEGER, INTENT(in) :: tile_count !< Tile number.
 
    INTEGER :: i
 
    find_input_field = diag_field_not_found ! Default return value if not found.
    DO i = 1, num_input_fields
       IF(tile_count == input_fields(i)%tile_count .AND.&
            & trim(input_fields(i)%module_name) == trim(module_name) .AND.&
            & lowercase(trim(input_fields(i)%field_name)) == lowercase(trim(field_name))) THEN
          find_input_field = i
          RETURN
       END IF
    END DO
  END FUNCTION find_input_field
 
  !> @brief Initialize the input field.
  SUBROUTINE init_input_field(module_name, field_name, tile_count)
    CHARACTER(len=*),  INTENT(in) :: module_name !< Module name.
    CHARACTER(len=*),  INTENT(in) :: field_name !< Input field name.
    INTEGER, INTENT(in) :: tile_count !< Tile number.
 
    ! Get a number for this input_field if not already set up
    IF ( find_input_field(module_name, field_name, tile_count) < 0 ) THEN
       num_input_fields = num_input_fields + 1
       IF ( num_input_fields > max_input_fields ) THEN
          ! <ERROR STATUS="FATAL">max_input_fields exceeded, increase it via diag_manager_nml</ERROR>
          CALL error_mesg('diag_util_mod::init_input_field',&
               & 'max_input_fields exceeded, increase it via diag_manager_nml', fatal)
       END IF
    ELSE
       ! If this is already initialized do not need to do anything
       RETURN
    END IF
 
    input_fields(num_input_fields)%module_name = trim(module_name)
    input_fields(num_input_fields)%field_name = trim(field_name)
    input_fields(num_input_fields)%num_output_fields = 0
    ! Set flag that this field has not been registered
    input_fields(num_input_fields)%register = .false.
    input_fields(num_input_fields)%local = .false.
    input_fields(num_input_fields)%standard_name = 'none'
    input_fields(num_input_fields)%tile_count = tile_count
    input_fields(num_input_fields)%numthreads = 1
    input_fields(num_input_fields)%active_omp_level = 0
    input_fields(num_input_fields)%time = time_zero
  END SUBROUTINE init_input_field
 
  !> @brief Initialize the output field.
  SUBROUTINE init_output_field(module_name, field_name, output_name, output_file,&
       & time_method, pack, tile_count, local_coord)
    CHARACTER(len=*), INTENT(in) :: module_name !< Module name.
    CHARACTER(len=*), INTENT(in) :: field_name !< Output field name.
    CHARACTER(len=*), INTENT(in) :: output_name !< Output name written to file.
    CHARACTER(len=*), INTENT(in) :: output_file !< File where field should be written.
    CHARACTER(len=*), INTENT(in) :: time_method !< Data reduction method.
                                                !! See <LINK SRC="diag_manager.html">diag_manager_mod</LINK>
                                                !! for valid methods.</IN>
    INTEGER, INTENT(in) :: pack !< Packing method.
    INTEGER, INTENT(in) :: tile_count !< Tile number.
    TYPE(coord_type), INTENT(in), OPTIONAL :: local_coord !< Region to be written.
                                                          !! If missing, then all data to be written.</IN>
    INTEGER :: out_num, in_num, file_num, file_num_tile1
    INTEGER :: num_fields, i, method_selected, l1
    INTEGER :: ioerror
    REAL :: pow_value
    INTEGER :: grv !< Value used to determine if the region defined in the diag_table is for the whole
                   !! axis, or a sub-axis
    CHARACTER(len=128) :: error_msg
    CHARACTER(len=50) :: t_method
    character(len=256) :: tmp_name
 
    ! Value to use to determine if a region is to be output on the full axis, or sub-axis
    ! get the value to compare to determine if writing full axis data
    IF ( region_out_use_alt_value ) THEN
       grv = glo_reg_val_alt
    ELSE
       grv = glo_reg_val
    END IF
 
 
    ! Get a number for this output field
    num_output_fields = num_output_fields + 1
    IF ( num_output_fields > max_output_fields ) THEN
       ! <ERROR STATUS="FATAL">max_output_fields = <max_output_fields> exceeded.  Increase via diag_manager_nml</ERROR>
       WRITE (unit=error_msg,fmt=*) max_output_fields
       CALL error_mesg('diag_util_mod::init_output_field', 'max_output_fields = '//trim(error_msg)//' exceeded.&
            &  Increase via diag_manager_nml', fatal)
    END IF
    out_num = num_output_fields
 
    ! First, find the index to the associated input field
    in_num = find_input_field(module_name, field_name, tile_count)
    IF ( in_num < 0 ) THEN
       IF ( tile_count > 1 ) THEN
          WRITE (error_msg,'(A,"/",A,"/",A)') trim(module_name),trim(field_name),&
               & "tile_count="//trim(string(tile_count))
       ELSE
          WRITE (error_msg,'(A,"/",A)') trim(module_name),trim(field_name)
       END IF
       ! <ERROR STATUS="FATAL">module_name/field_name <module_name>/<field_name>[/tile_count=<tile_count>]
       ! NOT registered</ERROR>
       CALL error_mesg('diag_util_mod::init_output_field',&
            & 'module_name/field_name '//trim(error_msg)//' NOT registered', fatal)
    END IF
 
    ! Add this output field into the list for this input field
    input_fields(in_num)%num_output_fields =&
         & input_fields(in_num)%num_output_fields + 1
    IF ( input_fields(in_num)%num_output_fields > max_out_per_in_field ) THEN
       ! <ERROR STATUS="FATAL">
       !   MAX_OUT_PER_IN_FIELD = <MAX_OUT_PER_IN_FIELD> exceeded for <module_name>/<field_name>,
       !   increase MAX_OUT_PER_IN_FIELD
       !   in the diag_manager_nml namelist.
       ! </ERROR>
       WRITE (unit=error_msg,fmt=*) max_out_per_in_field
       CALL error_mesg('diag_util_mod::init_output_field',&
        & 'MAX_OUT_PER_IN_FIELD exceeded for '//trim(module_name)//"/"//trim(field_name)//&
        &', increase MAX_OUT_PER_IN_FIELD in the diag_manager_nml namelist', fatal)
    END IF
    input_fields(in_num)%output_fields(input_fields(in_num)%num_output_fields) = out_num
 
    ! Also put pointer to input field in this output field
    output_fields(out_num)%input_field = in_num
 
    ! Next, find the number for the corresponding file
    IF ( trim(output_file).EQ.'null' ) THEN
       file_num = max_files
    ELSE
       file_num = find_file(output_file, 1)
       IF ( file_num < 0 ) THEN
          ! <ERROR STATUS="FATAL">
          !   file <file_name> is NOT found in the diag_table.
          ! </ERROR>
          CALL error_mesg('diag_util_mod::init_output_field', 'file '&
               & //trim(output_file)//' is NOT found in the diag_table', fatal)
       END IF
       IF ( tile_count > 1 ) THEN
          file_num_tile1 = file_num
          file_num = find_file(output_file, tile_count)
          IF(file_num < 0) THEN
             CALL init_file(files(file_num_tile1)%name, files(file_num_tile1)%output_freq,&
                  & files(file_num_tile1)%output_units, files(file_num_tile1)%format,&
                  & files(file_num_tile1)%time_units, files(file_num_tile1)%long_name,&
                  & tile_count, files(file_num_tile1)%new_file_freq,&
                  & files(file_num_tile1)%new_file_freq_units, files(file_num_tile1)%start_time,&
                  & files(file_num_tile1)%duration, files(file_num_tile1)%duration_units  )
             file_num = find_file(output_file, tile_count)
             IF ( file_num < 0 ) THEN
                ! <ERROR STATUS="FATAL">
                !   file <output_file> is not initialized for tile_count = <tile_count>
                ! </ERROR>
                CALL error_mesg('diag_util_mod::init_output_field', 'file '//trim(output_file)//&
                     & ' is not initialized for tile_count = '//trim(string(tile_count)), fatal)
             END IF
          END IF
       END IF
    END IF
 
    ! Insert this field into list for this file
    files(file_num)%num_fields = files(file_num)%num_fields + 1
    IF ( files(file_num)%num_fields > max_fields_per_file ) THEN
       WRITE (unit=error_msg, fmt=*) max_fields_per_file
       ! <ERROR STATUS="FATAL">
       !   MAX_FIELDS_PER_FILE = <MAX_FIELDS_PER_FILE> exceeded.  Increase MAX_FIELDS_PER_FILE in diag_data.F90.
       ! </ERROR>
       CALL error_mesg('diag_util_mod::init_output_field',&
            & 'MAX_FIELDS_PER_FILE = '//trim(error_msg)// &
             & ' exceeded.  Increase MAX_FIELDS_PER_FILE in diag_data.F90.', fatal)
    END IF
    num_fields = files(file_num)%num_fields
    files(file_num)%fields(num_fields) = out_num
 
    ! Set the file for this output field
    output_fields(out_num)%output_file = file_num
 
    ! Enter the other data for this output field
    output_fields(out_num)%output_name = trim(output_name)
    output_fields(out_num)%pack = pack
    output_fields(out_num)%pow_value = 1
    output_fields(out_num)%num_axes = 0
    output_fields(out_num)%total_elements = 0
    output_fields(out_num)%region_elements = 0
 
    call output_fields(out_num)%buff_bounds%reset(very_large_axis_length, 0)
 
    ! initialize the size of the diurnal axis to 1
    output_fields(out_num)%n_diurnal_samples = 1
 
    ! Initialize all time method to false
    method_selected = 0
    output_fields(out_num)%time_average = .false.
    output_fields(out_num)%time_rms = .false.
    output_fields(out_num)%time_min = .false.
    output_fields(out_num)%time_max = .false.
    output_fields(out_num)%time_sum = .false.
    output_fields(out_num)%time_ops = .false.
    output_fields(out_num)%written_once = .false.
 
    t_method = lowercase(time_method)
    ! cannot time average fields output every time
    IF ( files(file_num)%output_freq == every_time ) THEN
       output_fields(out_num)%time_average = .false.
       method_selected = method_selected+1
       t_method = 'point'
    ELSEIF ( index(t_method,'diurnal') == 1 ) THEN
       ! get the integer number from the t_method
       READ (unit=t_method(8:len_trim(t_method)), fmt=*, iostat=ioerror) output_fields(out_num)%n_diurnal_samples
       IF ( ioerror /= 0 ) THEN
          ! <ERROR STATUS="FATAL">
          !   could not find integer number of diurnal samples in string "<t_method>"
          ! </ERROR>
          CALL error_mesg('diag_util_mod::init_output_field',&
               & 'could not find integer number of diurnal samples in string "' //trim(t_method)//'"', fatal)
       ELSE IF ( output_fields(out_num)%n_diurnal_samples <= 0 ) THEN
          ! <ERROR STATUS="FATAL">
          !   The integer value of diurnal samples must be greater than zero.
          ! </ERROR>
          CALL error_mesg('diag_util_mod::init_output_field',&
               & 'The integer value of diurnal samples must be greater than zero.', fatal)
       END IF
       output_fields(out_num)%time_average = .true.
       method_selected = method_selected+1
       t_method='mean'
    ELSEIF ( index(t_method,'pow') == 1 ) THEN
       ! Get the integer number from the t_method
       READ (unit=t_method(4:len_trim(t_method)), fmt=*, iostat=ioerror) pow_value
       IF ( ioerror /= 0 .OR. output_fields(out_num)%pow_value < 1  .OR. floor(pow_value) /= ceiling(pow_value) ) THEN
          ! <ERROR STATUS="FATAL">
          !   Invalid power number in time operation "<t_method>".  Must be a positive integer.
          ! </ERROR>
          CALL error_mesg('diag_util_mod::init_output_field',&
               & 'Invalid power number in time operation "'//trim(t_method)//'".  Must be a positive integer', fatal)
       END IF
       output_fields(out_num)%pow_value = int(pow_value)
       output_fields(out_num)%time_average = .true.
       method_selected = method_selected+1
       t_method = 'mean_pow('//t_method(4:len_trim(t_method))//')'
    ELSE
       SELECT CASE(trim(t_method))
       CASE ( '.true.', 'mean', 'average', 'avg' )
          output_fields(out_num)%time_average = .true.
          method_selected = method_selected+1
          t_method = 'mean'
       CASE ( 'rms' )
          output_fields(out_num)%time_average = .true.
          output_fields(out_num)%time_rms = .true.
          output_fields(out_num)%pow_value = 2.0
          method_selected = method_selected+1
          t_method = 'root_mean_square'
       CASE ( '.false.', 'none', 'point' )
          output_fields(out_num)%time_average = .false.
          method_selected = method_selected+1
          t_method = 'point'
       CASE ( 'maximum', 'max' )
          output_fields(out_num)%time_max = .true.
          l1 = len_trim(output_fields(out_num)%output_name)
          if (l1 .ge. 3) then
              tmp_name = trim(adjustl(output_fields(out_num)%output_name(l1-2:l1)))
              IF (lowercase(trim(tmp_name)) /= 'max' ) then
                  output_fields(out_num)%output_name = trim(output_name)//'_max'
               endif
          endif
          method_selected = method_selected+1
          t_method = 'max'
       CASE ( 'minimum', 'min' )
          output_fields(out_num)%time_min = .true.
          l1 = len_trim(output_fields(out_num)%output_name)
          if (l1 .ge. 3) then
              tmp_name = trim(adjustl(output_fields(out_num)%output_name(l1-2:l1)))
              IF (lowercase(trim(tmp_name)) /= 'min' ) then
                  output_fields(out_num)%output_name = trim(output_name)//'_min'
              endif
          endif
          method_selected = method_selected+1
          t_method = 'min'
       CASE ( 'sum', 'cumsum' )
          output_fields(out_num)%time_sum = .true.
          l1 = len_trim(output_fields(out_num)%output_name)
          IF ( output_fields(out_num)%output_name(l1-2:l1) /= 'sum' )&
               & output_fields(out_num)%output_name = trim(output_name)//'_sum'
          method_selected = method_selected+1
          t_method = 'sum'
       END SELECT
    END IF
 
    ! reconcile logical flags
    output_fields(out_num)%time_ops = output_fields(out_num)%time_min.OR.output_fields(out_num)%time_max&
         & .OR.output_fields(out_num)%time_average .OR. output_fields(out_num)%time_sum
 
    output_fields(out_num)%phys_window = .false.
    ! need to initialize grid_type = -1(start, end, l_start_indx,l_end_indx etc...)
    IF ( PRESENT(local_coord) ) THEN
       input_fields(in_num)%local = .true.
       input_fields(in_num)%local_coord = local_coord
       IF ( int(local_coord%xbegin) ==  grv .AND. int(local_coord%xend) == grv .AND.&
            & int(local_coord%ybegin) == grv .AND. int(local_coord%yend) == grv ) THEN
          output_fields(out_num)%local_output = .false.
          output_fields(out_num)%need_compute = .false.
          output_fields(out_num)%reduced_k_range = .true.
       ELSE
          output_fields(out_num)%local_output = .true.
          output_fields(out_num)%need_compute = .false.
          output_fields(out_num)%reduced_k_range = .false.
       END IF
 
       output_fields(out_num)%output_grid%start(1) = local_coord%xbegin
       output_fields(out_num)%output_grid%start(2) = local_coord%ybegin
       output_fields(out_num)%output_grid%start(3) = local_coord%zbegin
       output_fields(out_num)%output_grid%end(1) = local_coord%xend
       output_fields(out_num)%output_grid%end(2) = local_coord%yend
       output_fields(out_num)%output_grid%end(3) = local_coord%zend
       DO i = 1, 3
          output_fields(out_num)%output_grid%l_start_indx(i) = -1
          output_fields(out_num)%output_grid%l_end_indx(i) = -1
          output_fields(out_num)%output_grid%subaxes(i) = -1
       END DO
    ELSE
       output_fields(out_num)%local_output = .false.
       output_fields(out_num)%need_compute = .false.
       output_fields(out_num)%reduced_k_range = .false.
    END IF
 
    ! <ERROR STATUS="FATAL">
    !   improper time method in diag_table for output field <output_name>
    ! </ERROR>
    IF ( method_selected /= 1 ) CALL error_mesg('diag_util_mod::init_output_field',&
         &'improper time method in diag_table for output field:'//trim(output_name),fatal)
 
    output_fields(out_num)%time_method = trim(t_method)
 
    ! allocate counters: NOTE that for simplicity we always allocate them, even
    ! if they are superceeded by 4D "counter" array. This isn't most memory
    ! efficient, approach, but probably tolerable since they are so small anyway
    ALLOCATE(output_fields(out_num)%count_0d(output_fields(out_num)%n_diurnal_samples))
    ALLOCATE(output_fields(out_num)%num_elements(output_fields(out_num)%n_diurnal_samples))
    output_fields(out_num)%count_0d(:) = 0
    output_fields(out_num)%num_elements(:) = 0
    output_fields(out_num)%num_attributes = 0
  END SUBROUTINE init_output_field
 
  !> @brief Open file for output, and write the meta data.
  SUBROUTINE opening_file(file, time, filename_time)
    ! WARNING: Assumes that all data structures are fully initialized
    INTEGER, INTENT(in) :: file !< File ID.
    TYPE(time_type), INTENT(in) :: time !< Time for the file time stamp
    TYPE(time_type), INTENT(in), optional :: filename_time !< Time used in setting the filename when
                                                           !! writting periodic files
 
    TYPE(time_type) :: fname_time !< Time used in setting the filename when writting periodic files
    REAL, DIMENSION(2) :: open_file_data
    INTEGER :: j, field_num, input_field_num, num_axes, k
    INTEGER :: field_num1
    INTEGER :: position
    INTEGER :: dir, edges
    INTEGER :: year, month, day, hour, minute, second
    INTEGER, DIMENSION(1) :: time_axis_id, time_bounds_id
    ! size of this axes array must be at least max num. of
    ! axes per field + 2; the last two elements are for time
    ! and time bounds dimensions
    INTEGER, DIMENSION(6) :: axes
    INTEGER, ALLOCATABLE  :: axesc(:) ! indices if compressed axes associated with the field
    LOGICAL :: time_ops, aux_present, match_aux_name, req_present, match_req_fields
    CHARACTER(len=7) :: avg_name = 'average'
    CHARACTER(len=128) :: time_units, timeb_units, avg, error_string, filename, aux_name, req_fields, fieldname
    CHARACTER(len=128) :: suffix, base_name
    CHARACTER(len=32) :: time_name, timeb_name,time_longname, timeb_longname, cart_name
    CHARACTER(len=256) :: fname
    CHARACTER(len=24) :: start_date
    TYPE(domain1d) :: domain
    TYPE(domain2d) :: domain2
    TYPE(domainug) :: domainU
    INTEGER :: is, ie, last, ind
    class(fmsnetcdffile_t), pointer    :: fileob
    integer :: actual_num_axes !< The actual number of axes to write including time
 
    aux_present = .false.
    match_aux_name = .false.
    req_present = .false.
    match_req_fields = .false.
 
    ! Here is where time_units string must be set up; time since base date
    WRITE (time_units, 11) trim(time_unit_list(files(file)%time_units)), get_base_year(),&
         & get_base_month(), get_base_day(), get_base_hour(), get_base_minute(), get_base_second()
11  FORMAT(a, ' since ', i4.4, '-', i2.2, '-', i2.2, ' ', i2.2, ':', i2.2, ':', i2.2)
    base_name = files(file)%name
    IF ( files(file)%new_file_freq < very_large_file_freq ) THEN
       position = index(files(file)%name, '%')
       IF ( position > 0 )  THEN
          base_name = base_name(1:position-1)
       ELSE
          ! <ERROR STATUS="FATAL">
          !   filename <files(file)%name> does not contain % for time stamp string
          ! </ERROR>
          CALL error_mesg('diag_util_mod::opening_file',&
               & 'file name '//trim(files(file)%name)//' does not contain % for time stamp string', fatal)
       END IF
       if (present(filename_time)) then
          fname_time = filename_time
       else
          fname_time = time
       endif
       suffix = get_time_string(files(file)%name, fname_time)
    ELSE
       suffix = ' '
    END IF
 
    ! Add ensemble ID to filename
    fname=base_name
    call get_instance_filename(fname, base_name)
 
    ! Set the filename
    filename = trim(base_name)//trim(suffix)
 
    ! prepend the file start date if prepend_date == .TRUE.
    IF ( prepend_date ) THEN
       call get_date(diag_init_time, year, month, day, hour, minute, second)
       write (start_date, '(1I20.4, 2I2.2)') year, month, day
 
       filename = trim(adjustl(start_date))//'.'//trim(filename)
    END IF
 
    ! Loop through all fields with this file to output axes
    ! JWD: This is a klooge; need something more robust
    domain2 = null_domain2d
    domainu = null_domainug
    DO j = 1, files(file)%num_fields
       field_num = files(file)%fields(j)
       if (output_fields(field_num)%local_output .AND. .NOT. output_fields(field_num)%need_compute) cycle
       num_axes = output_fields(field_num)%num_axes
       IF ( num_axes > 1 ) THEN
          domain2 = get_domain2d( output_fields(field_num)%axes(1:num_axes) )
          domainu = get_domainug( output_fields(field_num)%axes(1) )
          IF ( domain2 .NE. null_domain2d ) EXIT
       ELSEIF (num_axes == 1) THEN
          if (domainu .EQ. null_domainug) then
               domainu = get_domainug( output_fields(field_num)%axes(num_axes) )
          endif
       END IF
    END DO
 
    IF (domainu .NE. null_domainug .AND. domain2 .NE. null_domain2d) THEN
        CALL error_mesg('diag_util_mod::opening_file', &
                        'Domain2 and DomainU are somehow both set.', &
                            fatal)
    ENDIF
 
    IF ( allocated(files(file)%attributes) ) THEN
                CALL diag_output_init(filename, global_descriptor,&
                & files(file)%file_unit, domain2, domainu,&
                & fileobj(file),fileobju(file), fileobjnd(file), fnum_for_domain(file),&
                & attributes=files(file)%attributes(1:files(file)%num_attributes))
    ELSE
                CALL diag_output_init(filename, global_descriptor,&
                & files(file)%file_unit, domain2,domainu, &
                & fileobj(file),fileobju(file),fileobjnd(file),fnum_for_domain(file))
    END IF
    !> update fnum_for_domain with the correct domain
    files(file)%bytes_written = 0
    ! Does this file contain time_average fields?
    time_ops = .false.
    DO j = 1, files(file)%num_fields
       field_num = files(file)%fields(j)
       IF ( output_fields(field_num)%time_ops ) THEN
          time_ops = .true.
          EXIT
       END IF
    END DO
    ! Loop through all fields with this file to output axes
    DO j = 1, files(file)%num_fields
       field_num = files(file)%fields(j)
       input_field_num = output_fields(field_num)%input_field
       IF (.NOT.input_fields(input_field_num)%register) THEN
          WRITE (error_string,'(A,"/",A)') trim(input_fields(input_field_num)%module_name),&
               & trim(input_fields(input_field_num)%field_name)
          IF(mpp_pe() .EQ. mpp_root_pe()) THEN
             CALL error_mesg('diag_util_mod::opening_file',&
                  & 'module/field_name ('//trim(error_string)//') NOT registered', warning)
          END IF
          cycle
       END IF
       if (output_fields(field_num)%local_output .AND. .NOT. output_fields(field_num)%need_compute) cycle
 
       ! Put the time axis in the axis field
       num_axes = output_fields(field_num)%num_axes
       axes(1:num_axes) = output_fields(field_num)%axes(1:num_axes)
       ! make sure that axis_id are not -1
       DO k = 1, num_axes
          IF ( axes(k) < 0 ) THEN
             WRITE(error_string,'(a)') output_fields(field_num)%output_name
             ! <ERROR STATUS="FATAL">
             !   ouptut_name <output_fields(field_num)%output_name> has axis_id = -1
             ! </ERROR>
             CALL error_mesg('diag_util_mod::opening_file','output_name '//trim(error_string)//&
                  & ' has axis_id = -1', fatal)
          END IF
       END DO
       ! check if aux is present in any axes
       IF ( .NOT.aux_present ) THEN
          DO k = 1, num_axes
             aux_name = get_axis_aux(axes(k))
             IF ( trim(aux_name) /= 'none' ) THEN
                aux_present = .true.
                EXIT
             END IF
          END DO
       END IF
       ! check if required fields are present in any axes
       IF ( .NOT.req_present ) THEN
          DO k = 1, num_axes
             req_fields = get_axis_reqfld(axes(k))
             IF ( trim(req_fields) /= 'none' ) THEN
                CALL error_mesg('diag_util_mod::opening_file','required fields found: '//&
                               &trim(req_fields)//' in file '//trim(files(file)%name),note)
                req_present = .true.
                EXIT
             END IF
          END DO
       END IF
 
       axes(num_axes + 1) = files(file)%time_axis_id
!> Allocate the is_time_axis_registered field and set it to false for the first trip
       if (.not. allocated(files(file)%is_time_axis_registered)) then
          allocate(files(file)%is_time_axis_registered)
          files(file)%is_time_axis_registered = .false.
       endif
       if (time_ops) then
            !< If the file contains time_average fields write the "time" and "nv" dimension
            actual_num_axes = num_axes + 2
            axes(num_axes + 2) = files(file)%time_bounds_id
       else
            !< If the file doesn't contain time_average fields write the "time" dimension
            actual_num_axes = num_axes + 1
       endif
 
       if (fnum_for_domain(file) == "2d") then
          CALL write_axis_meta_data(files(file)%file_unit, axes(1:actual_num_axes),fileobj(file), time_ops=time_ops, &
                                   time_axis_registered=files(file)%is_time_axis_registered)
       elseif (fnum_for_domain(file) == "nd") then
          CALL write_axis_meta_data(files(file)%file_unit, axes(1:actual_num_axes),fileobjnd(file), time_ops=time_ops,&
                                   time_axis_registered=files(file)%is_time_axis_registered)
       elseif (fnum_for_domain(file) == "ug") then
          CALL write_axis_meta_data(files(file)%file_unit, axes(1:actual_num_axes),fileobju(file), time_ops=time_ops, &
                                   time_axis_registered=files(file)%is_time_axis_registered)
       endif
 
       ! write metadata for axes used  in compression-by-gathering, e.g. for unstructured
       ! grid
       DO k = 1, num_axes
          IF (axis_is_compressed(axes(k))) THEN
             CALL get_compressed_axes_ids(axes(k), axesc) ! returns allocatable array
             if (fnum_for_domain(file) == "ug") then
                 CALL write_axis_meta_data(files(file)%file_unit, axesc,fileobju(file), &
                                   time_axis_registered=files(file)%is_time_axis_registered)
             else
                 CALL error_mesg('diag_util_mod::opening_file::'//trim(filename), "Compressed "//&
                     "dimensions are only allowed with axis in the unstructured dimension", fatal)
             endif
             DEALLOCATE(axesc)
          ENDIF
       ENDDO
    END DO
 
    ! Looking for the first NON-static field in a file
    field_num1 = files(file)%fields(1)
    DO j = 1, files(file)%num_fields
       field_num = files(file)%fields(j)
       IF ( output_fields(field_num)%time_ops ) THEN
          field_num1 = field_num
          EXIT
       END IF
    END DO
    nfields_loop: DO j = 1, files(file)%num_fields
       field_num = files(file)%fields(j)
       input_field_num = output_fields(field_num)%input_field
       IF (.NOT.input_fields(input_field_num)%register) cycle
       IF (output_fields(field_num)%local_output .AND. .NOT. output_fields(field_num)%need_compute) cycle
       ! Make sure that 1 file contains either time_average or instantaneous fields
       ! cannot have both time_average and instantaneous in 1 file
       IF ( .NOT.mix_snapshot_average_fields ) THEN
          IF ( (output_fields(field_num)%time_ops.NEQV.output_fields(field_num1)%time_ops) .AND.&
               & .NOT.output_fields(field_num1)%static .AND. .NOT.output_fields(field_num)%static) THEN
             IF ( mpp_pe() == mpp_root_pe() ) THEN
                ! <ERROR STATUS="FATAL">
                !   <files(file)%name> can NOT have BOTH time average AND instantaneous fields.
                !   Create a new file or set mix_snapshot_average_fields=.TRUE. in the namelist diag_manager_nml.
                ! </ERROR>
                CALL error_mesg('diag_util_mod::opening_file','file '//trim(files(file)%name)// &
                     &' can NOT have BOTH time average AND instantaneous fields.'//&
                     &' Create a new file or set mix_snapshot_average_fields=.TRUE. in the namelist diag_manager_nml.'&
                     &, fatal)
             END IF
          END IF
       END IF
       ! check if any field has the same name as aux_name
       IF ( aux_present .AND. .NOT.match_aux_name ) THEN
          fieldname = output_fields(field_num)%output_name
          IF ( index(aux_name, trim(fieldname)) > 0 ) match_aux_name = .true.
       END IF
       ! check if any field has the same name as req_fields
       IF ( req_present .AND. .NOT.match_req_fields ) THEN
          fieldname = output_fields(field_num)%output_name
          is = 1; last = len_trim(req_fields)
          DO
            ind = index(req_fields(is:last),' ')
            IF (ind .eq. 0) ind = last-is+2
            ie = is+(ind-2)
            if (req_fields(is:ie) .EQ. trim(fieldname)) then
              match_req_fields = .true.
             !CALL error_mesg('diag_util_mod::opening_file','matched required field: '//TRIM(fieldname),NOTE)
              EXIT
            END IF
            is = is+ind
            if (is .GT. last) EXIT
          END DO
       END IF
 
       ! Put the time axis in the axis field
       num_axes = output_fields(field_num)%num_axes
       axes(1:num_axes) = output_fields(field_num)%axes(1:num_axes)
       IF ( .NOT.output_fields(field_num)%static ) THEN
          num_axes=num_axes+1
          axes(num_axes) = files(file)%time_axis_id
       END IF
       IF(output_fields(field_num)%time_average) THEN
          avg = avg_name
       ELSE IF(output_fields(field_num)%time_max) THEN
          avg = avg_name
       ELSE IF(output_fields(field_num)%time_min) THEN
          avg = avg_name
       ELSE
          avg = " "
       END IF
!> Use the correct file object
       if (fnum_for_domain(file) == "2d") then
          fileob => fileobj(file)
       elseif (fnum_for_domain(file) == "nd") then
          fileob => fileobjnd(file)
       elseif (fnum_for_domain(file) == "ug") then
          fileob => fileobju(file)
       endif
       IF ( input_fields(input_field_num)%missing_value_present ) THEN
          IF ( len_trim(input_fields(input_field_num)%interp_method) > 0 ) THEN
             output_fields(field_num)%f_type = write_field_meta_data(files(file)%file_unit,&
                  & output_fields(field_num)%output_name, axes(1:num_axes),&
                  & input_fields(input_field_num)%units,&
                  & input_fields(input_field_num)%long_name,&
                  & input_fields(input_field_num)%range, output_fields(field_num)%pack,&
                  & input_fields(input_field_num)%missing_value, avg_name = avg,&
                  & time_method=output_fields(field_num)%time_method,&
                  & standard_name = input_fields(input_field_num)%standard_name,&
                  & interp_method = input_fields(input_field_num)%interp_method,&
                  & attributes=output_fields(field_num)%attributes,&
                  & num_attributes=output_fields(field_num)%num_attributes,&
                  & use_ugdomain=files(file)%use_domainUG , &
                  & fileob=fileob)
          ELSE
             output_fields(field_num)%f_type = write_field_meta_data(files(file)%file_unit,&
                  & output_fields(field_num)%output_name, axes(1:num_axes),&
                  & input_fields(input_field_num)%units,&
                  & input_fields(input_field_num)%long_name,&
                  & input_fields(input_field_num)%range, output_fields(field_num)%pack,&
                  & input_fields(input_field_num)%missing_value, avg_name = avg,&
                  & time_method=output_fields(field_num)%time_method,&
                  & standard_name = input_fields(input_field_num)%standard_name,&
                  & attributes=output_fields(field_num)%attributes,&
                  & num_attributes=output_fields(field_num)%num_attributes,&
                  & use_ugdomain=files(file)%use_domainUG , &
                  & fileob=fileob)
 
          END IF
          ! NEED TO TAKE CARE OF TIME AVERAGING INFO TOO BOTH CASES
       ELSE
          IF ( len_trim(input_fields(input_field_num)%interp_method) > 0 ) THEN
             output_fields(field_num)%f_type = write_field_meta_data(files(file)%file_unit,&
                  & output_fields(field_num)%output_name, axes(1:num_axes),&
                  & input_fields(input_field_num)%units,&
                  & input_fields(input_field_num)%long_name,&
                  & input_fields(input_field_num)%range, output_fields(field_num)%pack,&
                  & avg_name = avg,&
                  & time_method=output_fields(field_num)%time_method,&
                  & standard_name = input_fields(input_field_num)%standard_name,&
                  & interp_method = input_fields(input_field_num)%interp_method,&
                  & attributes=output_fields(field_num)%attributes,&
                  & num_attributes=output_fields(field_num)%num_attributes,&
                  & use_ugdomain=files(file)%use_domainUG , &
                  & fileob=fileob)
 
          ELSE
             output_fields(field_num)%f_type = write_field_meta_data(files(file)%file_unit,&
                  & output_fields(field_num)%output_name, axes(1:num_axes),&
                  & input_fields(input_field_num)%units,&
                  & input_fields(input_field_num)%long_name,&
                  & input_fields(input_field_num)%range, output_fields(field_num)%pack,&
                  & avg_name = avg,&
                  & time_method=output_fields(field_num)%time_method,&
                  & standard_name = input_fields(input_field_num)%standard_name,&
                  & attributes=output_fields(field_num)%attributes,&
                  & num_attributes=output_fields(field_num)%num_attributes,&
                  & use_ugdomain=files(file)%use_domainUG , &
                  & fileob=fileob)
 
          END IF
       END IF
    END DO nfields_loop
    ! If any of the fields in the file are time averaged, need to output the axes
    ! Use double precision since time axis is double precision
    IF ( time_ops ) THEN
       time_axis_id(1) = files(file)%time_axis_id
       files(file)%f_avg_start = write_field_meta_data(files(file)%file_unit,&
            & avg_name // '_T1', time_axis_id, time_units,&
            & "Start time for average period", pack=pack_size , &
            & fileob=fileob)
       files(file)%f_avg_end = write_field_meta_data(files(file)%file_unit,&
            & avg_name // '_T2', time_axis_id, time_units,&
            & "End time for average period", pack=pack_size , &
            & fileob=fileob)
       files(file)%f_avg_nitems = write_field_meta_data(files(file)%file_unit,&
            & avg_name // '_DT', time_axis_id,&
            & trim(time_unit_list(files(file)%time_units)),&
            & "Length of average period", pack=pack_size , &
            & fileob=fileob)
    END IF
 
    IF ( time_ops ) THEN
       time_axis_id(1) = files(file)%time_axis_id
       time_bounds_id(1) = files(file)%time_bounds_id
       CALL get_diag_axis( time_axis_id(1), time_name, time_units, time_longname,&
            & cart_name, dir, edges, domain, domainu, open_file_data)
       CALL get_diag_axis( time_bounds_id(1), timeb_name, timeb_units, timeb_longname,&
            & cart_name, dir, edges, domain, domainu, open_file_data)
       ! CF Compliance requires the unit on the _bnds axis is the same as 'time'
       files(file)%f_bounds =  write_field_meta_data(files(file)%file_unit,&
            & trim(time_name)//'_bnds', (/time_bounds_id,time_axis_id/),&
            & time_units, trim(time_name)//' axis boundaries', pack=pack_size , &
            & fileob=fileob)
    END IF
    ! Let lower levels know that all meta data has been sent
    call done_meta_data(files(file)%file_unit)
 
    IF( aux_present .AND. .NOT.match_aux_name ) THEN
       ! <ERROR STATUS="WARNING">
       !   one axis has auxiliary but the corresponding field is NOT
       !   found in file <file_name>
       ! </ERROR>
       IF ( mpp_pe() == mpp_root_pe() ) CALL error_mesg('diag_util_mod::opening_file',&
            &'one axis has auxiliary but the corresponding field is NOT found in file '// &
             & trim(files(file)%name), warning)
    END IF
    IF( req_present .AND. .NOT.match_req_fields ) THEN
       ! <ERROR STATUS="FATAL">
       !   one axis has required fields but the corresponding field is NOT
       !   found in file <file_name>
       ! </ERROR>
       IF ( mpp_pe() == mpp_root_pe() ) CALL error_mesg('diag_util_mod::opening_file',&
                  &'one axis has required fields ('//trim(req_fields)//') but the '// &
                  &'corresponding fields are NOT found in file '//trim(files(file)%name), fatal)
    END IF
    ! Clean up pointer
    if (associated(fileob)) nullify(fileob)
  END SUBROUTINE opening_file
 
  !> @brief Write data out to file, and if necessary flush the buffers.
  SUBROUTINE diag_data_out(file, field, dat, time, final_call_in, static_write_in, filename_time)
    INTEGER, INTENT(in) :: file !< File ID.
    INTEGER, INTENT(in) :: field !< Field ID.
    REAL, DIMENSION(:,:,:,:), INTENT(inout) :: dat !< Data to write out.
    TYPE(time_type), INTENT(in) :: time !< Current model time.
    LOGICAL, OPTIONAL, INTENT(in):: final_call_in !< <TT>.TRUE.</TT> if this is the last write for file.
    LOGICAL, OPTIONAL, INTENT(in):: static_write_in !< <TT>.TRUE.</TT> if static fields are to be written to file.
    type(time_type), intent(in), optional :: filename_time !< Time used in setting the filename when
                                                           !! writting periodic files
 
    LOGICAL :: final_call, do_write, static_write
    REAL :: dif, time_data(2, 1, 1, 1), dt_time(1, 1, 1, 1), start_dif, end_dif
    REAL :: time_in_file !< Time in file at the beginning of this call
 
    !< Save the current time in the file. If the time in the file is not the same as the
    !! current time, files(file)%rtime_current will be updated
    time_in_file = files(file)%rtime_current
 
    do_write = .true.
    final_call = .false.
    IF ( PRESENT(final_call_in) ) final_call = final_call_in
    static_write = .false.
    IF ( PRESENT(static_write_in) ) static_write = static_write_in
!> dif is the time as a real that is evaluated
    dif = get_date_dif(time, get_base_time(), files(file)%time_units)
 
    ! get file_unit, open new file and close curent file if necessary
    IF ( .NOT.static_write .OR. files(file)%file_unit < 0 ) &
       CALL check_and_open(file, time, do_write, filename_time=filename_time)
    IF ( .NOT.do_write ) RETURN  ! no need to write data
!> Set up the time index and write the correct time value to the time array
    if (dif > files(file)%rtime_current) then
     files(file)%time_index = files(file)%time_index + 1
     files(file)%rtime_current = dif
     if (fnum_for_domain(file) == "2d") then
          call diag_write_time (fileobj(file), files(file)%rtime_current, files(file)%time_index,   &
                                time_name=fileobj(file)%time_name)
     elseif (fnum_for_domain(file) == "ug") then
          call diag_write_time (fileobju(file), files(file)%rtime_current, files(file)%time_index,  &
                                time_name=fileobju(file)%time_name)
     elseif (fnum_for_domain(file) == "nd") then
          call diag_write_time (fileobjnd(file), files(file)%rtime_current, files(file)%time_index, &
                                time_name=fileobjnd(file)%time_name)
     else
          call error_mesg("diag_util_mod::diag_data_out","Error opening the file "//files(file)%name,fatal)
     endif
    elseif (dif < files(file)%rtime_current .and. .not.(static_write) ) then
     call error_mesg("diag_util_mod::diag_data_out","The time for the file "//trim(files(file)%name)//&
                    " has gone backwards. There may be missing values for some of the variables",note)
    endif
!> Write data
    call diag_field_write (output_fields(field)%output_name, dat, static_write, file, fileobju, &
                         fileobj, fileobjnd, fnum_for_domain(file), time_in=files(file)%time_index)
    ! record number of bytes written to this file
    files(file)%bytes_written = files(file)%bytes_written +&
         & (SIZE(dat,1)*SIZE(dat,2)*SIZE(dat,3))*(8/output_fields(field)%pack)
    IF ( .NOT.output_fields(field)%written_once ) output_fields(field)%written_once = .true.
    ! *** inserted this line because start_dif < 0 for static fields ***
    IF ( .NOT.output_fields(field)%static ) THEN
       start_dif = get_date_dif(output_fields(field)%last_output, get_base_time(),files(file)%time_units)
       IF ( .NOT.mix_snapshot_average_fields ) THEN
          end_dif = get_date_dif(output_fields(field)%next_output, get_base_time(), files(file)%time_units)
       ELSE
          end_dif = dif
       END IF
    END IF
 
    if (files(file)%rtime_current > time_in_file) then !< If time was written in this call
       if (output_fields(field)%time_ops) then !< If this is a time_average field
          ! Output the axes if this is first time-averaged field
          time_data(1, 1, 1, 1) = start_dif
          call diag_field_write (files(file)%f_avg_start%fieldname, time_data(1:1,:,:,:), static_write, file, &
                                 fileobju, fileobj, fileobjnd, &
                                 fnum_for_domain(file), time_in=files(file)%time_index)
          time_data(2, 1, 1, 1) = end_dif
          call diag_field_write (files(file)%f_avg_end%fieldname, time_data(2:2,:,:,:), static_write, file, &
                                 fileobju, fileobj, fileobjnd, &
                                 fnum_for_domain(file), time_in=files(file)%time_index)
          ! Compute the length of the average
          dt_time(1, 1, 1, 1) = end_dif - start_dif
          call diag_field_write (files(file)%f_avg_nitems%fieldname, dt_time(1:1,:,:,:), static_write, file, &
                                 fileobju, fileobj, fileobjnd, &
                                 fnum_for_domain(file), time_in=files(file)%time_index)
          ! Include boundary variable for CF compliance
          call diag_field_write (files(file)%f_bounds%fieldname, time_data(1:2,:,:,:), static_write, file, &
                                 fileobju, fileobj, fileobjnd, &
                                 fnum_for_domain(file), time_in=files(file)%time_index)
       END IF
    END IF
 
    ! If write time is greater (equal for the last call) than last_flush for this file, flush it
    IF ( final_call ) THEN
       IF ( time >= files(file)%last_flush ) THEN
          files(file)%last_flush = time
       END IF
    ELSE
       IF ( time > files(file)%last_flush .AND. (flush_nc_files.OR.debug_diag_manager) ) THEN
          call diag_flush(file, fileobju, fileobj, fileobjnd, fnum_for_domain(file))
          files(file)%last_flush = time
       END IF
    END IF
  END SUBROUTINE diag_data_out
 
  !> @brief Checks if it is time to open a new file.
  !! @details Checks if it is time to open a new file. If yes, it first closes the
  !!     current file, opens a new file and returns file_unit
  !!     previous diag_manager_end is replaced by closing_file and output_setup by opening_file.
  SUBROUTINE check_and_open(file, time, do_write, filename_time)
    INTEGER, INTENT(in) :: file !<File ID.
    TYPE(time_type), INTENT(in) :: time !< Current model time.
    LOGICAL, INTENT(out) :: do_write !< <TT>.TRUE.</TT> if file is expecting more data to write,
                                     !! <TT>.FALSE.</TT> otherwise.
    TYPE(time_type), INTENT(in), optional :: filename_time !< Time used in setting the filename when
                                                           !! writting periodic files
 
    IF ( time >= files(file)%start_time ) THEN
       IF ( files(file)%file_unit < 0 ) THEN ! need to open a new file
          CALL opening_file(file, time, filename_time=filename_time)
          do_write = .true.
       ELSE
          do_write = .true.
          IF ( time > files(file)%close_time .AND. time < files(file)%next_open ) THEN
             do_write = .false. ! file still open but receives NO MORE data
          ELSE IF ( time > files(file)%next_open ) THEN ! need to close current file and open a new one
             CALL write_static(file)  ! write all static fields and close this file
             CALL opening_file(file, time, filename_time=filename_time)
             files(file)%time_index = 0 !< Reset the number of times in the files back to 0
             files(file)%start_time = files(file)%next_open
             files(file)%close_time =&
                  & diag_time_inc(files(file)%start_time,files(file)%duration, files(file)%duration_units)
             files(file)%next_open =&
                  & diag_time_inc(files(file)%next_open, files(file)%new_file_freq,&
                  & files(file)%new_file_freq_units)
             IF ( files(file)%close_time > files(file)%next_open ) THEN
                ! <ERROR STATUS="FATAL">
                !   <file_name> has close time GREATER than next_open time,
                !   check file duration and frequency
                ! </ERROR>
                CALL error_mesg('diag_util_mod::check_and_open',&
                     & files(file)%name// &
                      & ' has close time GREATER than next_open time, check file duration and frequency',fatal)
             END IF
          END IF ! no need to open new file, simply return file_unit
       END IF
    ELSE
       do_write = .false.
    END IF
  END SUBROUTINE check_and_open
 
  !> @brief Output all static fields in this file
  SUBROUTINE write_static(file)
    INTEGER, INTENT(in) :: file !< File ID.
    INTEGER :: j, i, input_num
 
    DO j = 1, files(file)%num_fields
       i = files(file)%fields(j)
       input_num = output_fields(i)%input_field
       ! skip fields that were not registered
       IF ( .NOT.input_fields(input_num)%register ) cycle
       IF ( output_fields(i)%local_output .AND. .NOT. output_fields(i)%need_compute) cycle
       ! only output static fields here
       IF ( .NOT.output_fields(i)%static ) cycle
       CALL diag_data_out(file, i, output_fields(i)%buffer, files(file)%last_flush, .true., .true.)
    END DO
!! New FMS_IO close
      ! File is stil open.  This is to protect when the diag_table has no Fields
      ! going to this file, and it was never opened (b/c diag_data_out was not
      ! called)
      if (fnum_for_domain(file) == "2d" )then
          if (check_if_open(fileobj(file))) call close_file (fileobj(file) )
      elseif (fnum_for_domain(file) == "nd") then
          if (check_if_open(fileobjnd(file)) ) then
               call close_file (fileobjnd(file))
          endif
      elseif (fnum_for_domain(file) == "ug") then
          if (check_if_open(fileobju(file))) call close_file (fileobju(file))
      endif
      files(file)%file_unit = -1
  END SUBROUTINE write_static
 
  !> @brief Checks to see if <TT>output_name</TT> and <TT>output_file</TT> are unique in <TT>output_fields</TT>.
  SUBROUTINE check_duplicate_output_fields(err_msg)
    CHARACTER(len=*), INTENT(out), OPTIONAL :: err_msg !< Error message.  If empty, then no duplicates found.
 
    INTEGER :: i, j, tmp_file
    CHARACTER(len=128) :: tmp_name
    CHARACTER(len=256) :: err_msg_local
 
    IF ( PRESENT(err_msg) ) err_msg=''
    ! Do the checking when more than 1 output_fileds present
    IF ( num_output_fields <= 1 ) RETURN
    err_msg_local = ''
 
    i_loop: DO i = 1, num_output_fields-1
       tmp_name = trim(output_fields(i)%output_name)
       tmp_file =  output_fields(i)%output_file
       DO j = i+1, num_output_fields
          IF ( (tmp_name == trim(output_fields(j)%output_name)) .AND. &
               &(tmp_file == output_fields(j)%output_file)) THEN
             err_msg_local = ' output_field "'//trim(tmp_name)//&
                  &'" duplicated in file "'//trim(files(tmp_file)%name)//'"'
             EXIT i_loop
          END IF
       END DO
    END DO i_loop
    IF ( err_msg_local /= '' ) THEN
       IF ( fms_error_handler(' ERROR in diag_table',err_msg_local,err_msg) ) RETURN
    END IF
  END SUBROUTINE check_duplicate_output_fields
 
  !> @brief Allocates the atttype in out_field
  SUBROUTINE attribute_init_field(out_field, err_msg)
    TYPE(output_field_type), INTENT(inout) :: out_field !< output field to allocate memory for attribute
    CHARACTER(LEN=*), INTENT(out), OPTIONAL :: err_msg !< Error message, passed back to calling function
 
    INTEGER :: istat
 
    ! Need to initialize err_msg if present
    IF ( PRESENT(err_msg) ) err_msg = ''
 
    ! Allocate memory for the attributes
    IF ( .NOT.allocated(out_field%attributes) ) THEN
       ALLOCATE(out_field%attributes(max_field_attributes), stat=istat)
       IF ( istat.NE.0 ) THEN
          ! <ERROR STATUS="FATAL">
          !   Unable to allocate memory for attribute <name> to module/input_field <module_name>/<field_name>
          ! </ERROR>
          IF ( fms_error_handler('diag_util_mod::attribute_init_field',&
               & 'Unable to allocate memory for attributes', err_msg) ) THEN
             RETURN
          END IF
       ELSE
          ! Set equal to 0.  It will be increased when attributes added
          out_field%num_attributes = 0
       END IF
    END IF
  END SUBROUTINE attribute_init_field
 
  !> @brief Prepends the attribute value to an already existing attribute.  If the
  !!     attribute isn't yet defined, then creates a new attribute
  SUBROUTINE prepend_attribute_field(out_field, att_name, prepend_value, err_msg)
    TYPE(output_field_type), INTENT(inout) :: out_field !< output field that will get the attribute
    CHARACTER(len=*), INTENT(in) :: att_name !< Name of the attribute
    CHARACTER(len=*), INTENT(in) :: prepend_value !< Value to prepend
    CHARACTER(len=*), INTENT(out) , OPTIONAL :: err_msg !< Error message, passed back to calling routine
 
    INTEGER :: length, i, this_attribute
    CHARACTER(len=512) :: err_msg_local
 
    ! Initialize string characters
    err_msg_local=''
    IF ( PRESENT(err_msg) ) err_msg = ''
 
    ! Make sure the attributes for this out field have been initialized
    CALL attribute_init_field(out_field, err_msg_local)
    IF ( trim(err_msg_local) .NE. '' ) THEN
       IF ( fms_error_handler('diag_util_mod::prepend_attribute_field', trim(err_msg_local), err_msg) ) THEN
          RETURN
       END IF
    END IF
 
    ! Find if attribute exists
    this_attribute = 0
    DO i=1, out_field%num_attributes
       IF ( trim(out_field%attributes(i)%name) .EQ. trim(att_name) ) THEN
          this_attribute = i
          EXIT
       END IF
    END DO
 
    IF ( this_attribute > 0 ) THEN
       IF ( out_field%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_field', &
               & '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 field attributes
       this_attribute = out_field%num_attributes + 1
       IF ( this_attribute .GT. max_field_attributes ) THEN
          ! <ERROR STATUS="FATAL">
          !   Number of attributes exceeds max_field_attributes for attribute <name>.
          !   Increase diag_manager_nml:max_field_attributes.
          ! </ERROR>
          IF ( fms_error_handler('diag_util_mod::prepend_attribute_field',&
               & 'Number of attributes exceeds max_field_attributes for attribute "'&
               & //trim(att_name)//'".  Increase diag_manager_nml:max_field_attributes.',&
               & err_msg) ) THEN
             RETURN
          END IF
       ELSE
          out_field%num_attributes = this_attribute
          ! Set name and type
          out_field%attributes(this_attribute)%name = att_name
          out_field%attributes(this_attribute)%type = nf90_char
          ! Initialize catt to a blank string, as len_trim doesn't always work on an uninitialized string
          out_field%attributes(this_attribute)%catt = ''
       END IF
    END IF
 
    ! Check if string is already included, and return if found
    IF ( index(trim(out_field%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_field%attributes(this_attribute)%catt))
       IF ( length.GT.len(out_field%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_field',&
               & 'Prepend length for attribute "'//trim(att_name)//'" is longer than allowed.',&
               & err_msg) ) THEN
             RETURN
          END IF
       END IF
       ! Set fields
       out_field%attributes(this_attribute)%catt =&
            & trim(prepend_value)//' '//trim(out_field%attributes(this_attribute)%catt)
       out_field%attributes(this_attribute)%len = length
    END IF
  END SUBROUTINE prepend_attribute_field
  !> @brief Allocates the atttype in out_file
  SUBROUTINE attribute_init_file(out_file, err_msg)
    TYPE(file_type), INTENT(inout) :: out_file !< output file to allocate memory for attribute
    CHARACTER(LEN=*), INTENT(out), OPTIONAL :: err_msg !< Error message, passed back to calling function
 
    INTEGER :: istat
 
    ! Initialize err_msg
    IF ( PRESENT(err_msg) ) err_msg = ''
 
    ! Allocate memory for the attributes
    IF ( .NOT.allocated(out_file%attributes) ) THEN
       ALLOCATE(out_file%attributes(max_field_attributes), stat=istat)
       IF ( istat.NE.0 ) THEN
          ! <ERROR STATUS="FATAL">
          !   Unable to allocate memory for file attributes
          ! </ERROR>
          IF ( fms_error_handler('diag_util_mod::attribute_init_file', &
             &  'Unable to allocate memory for file attributes', err_msg) ) THEN
             RETURN
          END IF
       ELSE
          ! Set equal to 0.  It will be increased when attributes added
          out_file%num_attributes = 0
       END IF
    END IF
  END SUBROUTINE attribute_init_file
 
  !> @brief Prepends the attribute value to an already existing attribute.  If the
  !!     attribute isn't yet defined, then creates a new attribute
  SUBROUTINE prepend_attribute_file(out_file, att_name, prepend_value, err_msg)
    TYPE(file_type), INTENT(inout) :: out_file !< output file that will get the attribute
    CHARACTER(len=*), INTENT(in) :: att_name !< Name of the attribute
    CHARACTER(len=*), INTENT(in) :: prepend_value !< Value to prepend
    CHARACTER(len=*), INTENT(out) , OPTIONAL :: err_msg !< Error message, passed back to calling routine
 
    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_file(out_file, err_msg_local)
    IF ( trim(err_msg_local) .NE. '' ) THEN
       IF ( fms_error_handler('diag_util_mod::prepend_attribute_file', trim(err_msg_local), err_msg) ) THEN
          RETURN
       END IF
    END IF
 
    ! Find if attribute exists
    this_attribute = 0
    DO i=1, out_file%num_attributes
       IF ( trim(out_file%attributes(i)%name) .EQ. trim(att_name) ) THEN
          this_attribute = i
          EXIT
       END IF
    END DO
 
    IF ( this_attribute > 0 ) THEN
       IF ( out_file%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_file',&
               & '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_file%num_attributes + 1
       IF ( this_attribute .GT. max_file_attributes ) THEN
          ! <ERROR STATUS="FATAL">
          !   Number of attributes exceeds max_file_attributes for attribute <name>.
          !   Increase diag_manager_nml:max_file_attributes.
          ! </ERROR>
          IF ( fms_error_handler('diag_util_mod::prepend_attribute_file',&
               & 'Number of attributes exceeds max_file_attributes for attribute "'&
               &//trim(att_name)//'".  Increase diag_manager_nml:max_file_attributes.',&
               & err_msg) ) THEN
             RETURN
          END IF
       ELSE
          out_file%num_attributes = this_attribute
          ! Set name and type
          out_file%attributes(this_attribute)%name = att_name
          out_file%attributes(this_attribute)%type = nf90_char
          ! Initialize catt to a blank string, as len_trim doesn't always work on an uninitialized string
          out_file%attributes(this_attribute)%catt = ''
       END IF
    END IF
 
    ! Only add string only if not already defined
    IF ( index(trim(out_file%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_file%attributes(this_attribute)%catt))
       IF ( length.GT.len(out_file%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_file%attributes(this_attribute)%catt =&
            & trim(prepend_value)//' '//trim(out_file%attributes(this_attribute)%catt)
       out_file%attributes(this_attribute)%len = length
    END IF
  END SUBROUTINE prepend_attribute_file
 
  !> @brief Get the a diag_file's start_time as it is defined in the diag_table
  !! @return the start_time for the file
  function get_file_start_time(file_num) &
   result(start_time)
   integer,         intent(in)  :: file_num   !< File number of the file to get the start_time from
 
   TYPE(time_type) :: start_time !< The start_time to return
 
   start_time = files(file_num)%start_time
  end function get_file_start_time
END MODULE diag_util_mod
!> @}
! close documentation grouping