mpp_io_write.inc

Go to the documentation of this file.

! -*-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/>.
!***********************************************************************

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                            !
!                             MPP_WRITE_META                                 !
!                                                                            !
! This series of routines is used to describe the contents of the file       !
! being written on <unit>. Each file can contain any number of fields,       !
! which can be functions of 0-3 spatial axes and 0-1 time axes. Axis         !
! descriptors are stored in the <axistype> structure and field               !
! descriptors in the <fieldtype> structure.                                  !
!                                                                            !
!  type, public :: axistype                                                  !
!     sequence                                                               !
!     character(len=128) :: name                                             !
!     character(len=128) :: units                                            !
!     character(len=256) :: longname                                         !
!     integer :: sense           !+/-1, depth or height?                     !
!     type(domain1D) :: domain                                               !
!     real, pointer :: data(:) !axis values (not used if time axis)          !
!     integer :: id                                                          !
!  end type axistype                                                         !
!                                                                            !
!  type, public :: fieldtype                                                 !
!     sequence                                                               !
!     character(len=128) :: name                                             !
!     character(len=128) :: units                                            !
!     character(len=256) :: longname                                         !
!     character(len=256) :: standard_name  !CF standard name                 !
!     real :: min, max, missing, fill, scale, add                            !
!     type(axistype), pointer :: axis(:)                                     !
!     integer :: id                                                          !
!  end type fieldtype                                                        !
!                                                                            !
! The metadata contained in the type is always written for each axis and     !
! field. Any other metadata one wishes to attach to an axis or field         !
! can subsequently be passed to mpp_write_meta using the ID, as shown below. !
!                                                                            !
! mpp_write_meta can take several forms:                                     !
!                                                                            !
!  mpp_write_meta( unit, name, rval=rval, pack=pack )                        !
!  mpp_write_meta( unit, name, ival=ival )                                   !
!  mpp_write_meta( unit, name, cval=cval )                                   !
!      integer, intent(in) :: unit                                           !
!      character(len=*), intent(in) :: name                                  !
!      real, intent(in), optional :: rval(:)                                 !
!      integer, intent(in), optional :: ival(:)                              !
!      character(len=*), intent(in), optional :: cval                        !
!                                                                            !
! This form defines global metadata associated with the file as a            !
! whole. The attribute is named <name> and can take on a real, integer       !
! or character value. <rval> and <ival> can be scalar or 1D arrays.          !
!                                                                            !
!  mpp_write_meta( unit, id, name, rval=rval, pack=pack )                    !
!  mpp_write_meta( unit, id, name, ival=ival )                               !
!  mpp_write_meta( unit, id, name, cval=cval )                               !
!      integer, intent(in) :: unit, id                                       !
!      character(len=*), intent(in) :: name                                  !
!      real, intent(in), optional :: rval(:)                                 !
!      integer, intent(in), optional :: ival(:)                              !
!      character(len=*), intent(in), optional :: cval                        !
!                                                                            !
! This form defines metadata associated with a previously defined            !
! axis or field, identified to mpp_write_meta by its unique ID <id>.         !
! The attribute is named <name> and can take on a real, integer              !
! or character value. <rval> and <ival> can be scalar or 1D arrays.          !
! This need not be called for attributes already contained in                !
! the type.                                                                  !
!                                                                            !
! PACK can take values 1,2,4,8. This only has meaning when writing           !
! floating point numbers. The value of PACK defines the number of words      !
! written into 8 bytes. For pack=4 and pack=8, an integer value is           !
! written: rval is assumed to have been scaled to the appropriate dynamic    !
! range.                                                                     !
! PACK currently only works for netCDF files, and is ignored otherwise.      !
!                                                                            !
!  subroutine mpp_write_meta_axis( unit, axis, name, units, longname, &      !
!        cartesian, sense, domain, data )                                    !
!     integer, intent(in) :: unit                                            !
!     type(axistype), intent(inout) :: axis                                  !
!     character(len=*), intent(in) :: name, units, longname                  !
!     character(len=*), intent(in), optional :: cartesian                    !
!     integer, intent(in), optional :: sense                                 !
!     type(domain1D), intent(in), optional :: domain                         !
!     real, intent(in), optional :: data(:)                                  !
!                                                                            !
! This form defines a time or space axis. Metadata corresponding to the      !
! type above are written to the file on <unit>. A unique ID for subsequent   !
! references to this axis is returned in axis%id. If the <domain>            !
! element is present, this is recognized as a distributed data axis          !
! and domain decomposition information is also written if required (the      !
! domain decomposition info is required for multi-fileset multi-threaded     !
! I/O). If the <datLINK> element is allocated, it is considered to be a      !
! space axis, otherwise it is a time axis with an unlimited dimension.       !
! Only one time axis is allowed per file.                                    !
!                                                                            !
!  subroutine mpp_write_meta_field( unit, field, axes, name, units, longname !
!       stanadard_name, min, max, missing, fill, scale, add, pack)           !
!     integer, intent(in) :: unit                                            !
!     type(fieldtype), intent(out) :: field                                  !
!     type(axistype), intent(in) :: axes(:)                                  !
!     character(len=*), intent(in) :: name, units, longname, standard_name   !
!     real, intent(in), optional :: min, max, missing, fill, scale, add      !
!     integer, intent(in), optional :: pack                                  !
!                                                                            !
! This form defines a field. Metadata corresponding to the type              !
! above are written to the file on <unit>. A unique ID for subsequent        !
! references to this field is returned in field%id. At least one axis        !
! must be associated, 0D variables are not considered. mpp_write_meta        !
! must previously have been called on all axes associated with this          !
! field.                                                                     !
!                                                                            !
! The mpp_write_meta package also includes subroutines write_attribute and   !
! write_attribute_netcdf, that are private to this module.                   !
!                                                                            !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine mpp_write_meta_global( unit, name, rval, ival, cval, pack)
!writes a global metadata attribute to unit <unit>
!attribute <name> can be an real, integer or character
!one and only one of rval, ival, and cval should be present
!the first found will be used
!for a non-netCDF file, it is encoded into a string "GLOBAL <name> <val>"
      integer, intent(in) :: unit
      character(len=*), intent(in) :: name
      real,             intent(in), optional :: rval(:)
      integer,          intent(in), optional :: ival(:)
      character(len=*), intent(in), optional :: cval
      integer, intent(in), optional :: pack

!      call mpp_clock_begin(mpp_write_clock)
      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe) then
!         call mpp_clock_end(mpp_write_clock)
         return
      endif
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized ) &
           call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )

      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
#ifdef use_netCDF
          call write_attribute_netcdf( unit, NF_GLOBAL, name, rval, ival, cval, pack )
#endif
      else
          call write_attribute( unit, 'GLOBAL '//trim(name), rval, ival, cval, pack )
      end if
!      call mpp_clock_end(mpp_write_clock)

      return
    end subroutine mpp_write_meta_global

!versions of above to support <rval> and <ival> as scalars (because of f90 strict rank matching)
    subroutine mpp_write_meta_global_scalar_r( unit, name, rval, pack )
      integer, intent(in) :: unit
      character(len=*), intent(in) :: name
      real, intent(in) :: rval
      integer, intent(in), optional :: pack

      call mpp_write_meta_global( unit, name, rval=(/rval/), pack=pack )
      return
    end subroutine mpp_write_meta_global_scalar_r

    subroutine mpp_write_meta_global_scalar_i( unit, name, ival, pack )
      integer, intent(in) :: unit
      character(len=*), intent(in) :: name
      integer, intent(in) :: ival
      integer, intent(in), optional :: pack

      call mpp_write_meta_global( unit, name, ival=(/ival/), pack=pack )
      return
    end subroutine mpp_write_meta_global_scalar_i

    subroutine mpp_write_meta_var( unit, id, name, rval, ival, cval, pack)
!writes a metadata attribute for variable <id> to unit <unit>
!attribute <name> can be an real, integer or character
!one and only one of rval, ival, and cval should be present
!the first found will be used
!for a non-netCDF file, it is encoded into a string "<id> <name> <val>"
      integer, intent(in) :: unit, id
      character(len=*), intent(in) :: name
      real,             intent(in), optional :: rval(:)
      integer,          intent(in), optional :: ival(:)
      character(len=*), intent(in), optional :: cval
      integer,          intent(in), optional :: pack

      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe) then
         return
      endif
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized ) &
           call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )

      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
          call write_attribute_netcdf( unit, id, name, rval, ival, cval, pack )
      else
          write( text, '(a,i4,a)' )'VARIABLE ', id, ' '//name
          call write_attribute( unit, trim(text), rval, ival, cval, pack )
      end if

      return
    end subroutine mpp_write_meta_var

!versions of above to support <rval> and <ival> as scalar (because of f90 strict rank matching)
    subroutine mpp_write_meta_scalar_r( unit, id, name, rval, pack )
      integer, intent(in) :: unit, id
      character(len=*), intent(in) :: name
      real, intent(in) :: rval
      integer, intent(in), optional :: pack

      call mpp_write_meta( unit, id, name, rval=(/rval/), pack=pack )
      return
    end subroutine mpp_write_meta_scalar_r

    subroutine mpp_write_meta_scalar_i( unit, id, name, ival,pack )
      integer, intent(in) :: unit, id
      character(len=*), intent(in) :: name
      integer, intent(in) :: ival
      integer, intent(in), optional :: pack

      call mpp_write_meta( unit, id, name, ival=(/ival/),pack=pack )
      return
    end subroutine mpp_write_meta_scalar_i


    subroutine mpp_write_axis_data (unit, axes )
      integer, intent(in) :: unit
      type(axistype), dimension(:), intent(in) :: axes

      integer :: naxis

      naxis = size (axes)
      allocate (mpp_file(unit)%axis(naxis))
      mpp_file(unit)%axis(1:naxis) = axes(1:naxis)
#ifdef use_netCDF
      if( mpp_file(unit)%action.EQ.MPP_WRONLY )then
         if(header_buffer_val>0) then
            error = NF__ENDDEF(mpp_file(unit)%ncid,header_buffer_val,4,0,4)
         else
            error = NF_ENDDEF(mpp_file(unit)%ncid)
         endif
      endif
#endif
    end subroutine mpp_write_axis_data

    subroutine mpp_def_dim_nodata(unit,name,size)
      integer, intent(in) :: unit
      character(len=*), intent(in) :: name
      integer, intent(in) :: size
      integer :: error,did

    ! This routine assumes the file is in define mode
      if(.NOT. mpp_file(unit)%write_on_this_pe) return
#ifdef use_netCDF
      error = NF_DEF_DIM(mpp_file(unit)%ncid,name,size,did)
      call netcdf_err(error, mpp_file(unit),string='Axis='//trim(name))
#endif
    end subroutine mpp_def_dim_nodata

    subroutine mpp_def_dim_int(unit,name,dsize,longname,data)
      integer, intent(in) :: unit
      character(len=*), intent(in) :: name
      integer, intent(in) :: dsize
      character(len=*), intent(in) :: longname
      integer, intent(in) :: data(:)
      integer :: error,did,id

    ! This routine assumes the file is in define mode
#ifdef use_netCDF
      if(.NOT. mpp_file(unit)%write_on_this_pe) return
      error = NF_DEF_DIM(mpp_file(unit)%ncid,name,dsize,did)
      call netcdf_err(error, mpp_file(unit),string='Axis='//trim(name))

    ! Write dimension data.
      error = NF_DEF_VAR( mpp_file(unit)%ncid, name, NF_INT, 1, did, id )
      call netcdf_err( error, mpp_file(unit), string=' axis varable '//trim(name))

      error = NF_PUT_ATT_TEXT( mpp_file(unit)%ncid, id, 'long_name', len_trim(longname), longname )
      call netcdf_err( error, mpp_file(unit), string=' Attribute=long_name' )

      if( mpp_file(unit)%action.EQ.MPP_WRONLY )then
         if(header_buffer_val>0) then
            error = NF__ENDDEF(mpp_file(unit)%ncid,header_buffer_val,4,0,4)
         else
            error = NF_ENDDEF(mpp_file(unit)%ncid)
         endif
      endif
      call netcdf_err( error, mpp_file(unit), string=' subroutine mpp_def_dim')
      error = NF_PUT_VARA_INT ( mpp_file(unit)%ncid, id, 1, size(data), data )
      call netcdf_err( error, mpp_file(unit), string=' axis varable '//trim(name))
      error = NF_REDEF(mpp_file(unit)%ncid)
      call netcdf_err( error, mpp_file(unit), string=' subroutine mpp_def_dim')
#endif
    return
  end subroutine mpp_def_dim_int

    subroutine mpp_def_dim_real(unit,name,dsize,longname,data)
      integer, intent(in) :: unit
      character(len=*), intent(in) :: name
      integer, intent(in) :: dsize
      character(len=*), intent(in) :: longname
      real, intent(in) :: data(:)
      integer :: error,did,id

    ! This routine assumes the file is in define mode
#ifdef use_netCDF
      if(.NOT. mpp_file(unit)%write_on_this_pe) return
      error = NF_DEF_DIM(mpp_file(unit)%ncid,name,dsize,did)
      call netcdf_err(error, mpp_file(unit),string='Axis='//trim(name))

    ! Write dimension data.
      error = NF_DEF_VAR( mpp_file(unit)%ncid, name, NF_INT, 1, did, id )
      call netcdf_err( error, mpp_file(unit), string=' axis varable '//trim(name))

      error = NF_PUT_ATT_TEXT( mpp_file(unit)%ncid, id, 'long_name', len_trim(longname), longname )
      call netcdf_err( error, mpp_file(unit), string=' Attribute=long_name' )

      if( mpp_file(unit)%action.EQ.MPP_WRONLY )then
         if(header_buffer_val>0) then
            error = NF__ENDDEF(mpp_file(unit)%ncid,header_buffer_val,4,0,4)
         else
            error = NF_ENDDEF(mpp_file(unit)%ncid)
         endif
      endif
      call netcdf_err( error, mpp_file(unit), string=' subroutine mpp_def_dim')
      error = NF_PUT_VARA_INT ( mpp_file(unit)%ncid, id, 1, size(data), data )
      call netcdf_err( error, mpp_file(unit), string=' axis varable '//trim(name))
      error = NF_REDEF(mpp_file(unit)%ncid)
      call netcdf_err( error, mpp_file(unit), string=' subroutine mpp_def_dim')
#endif
    return
  end subroutine mpp_def_dim_real



    subroutine mpp_write_meta_axis_r1d( unit, axis, name, units, longname, cartesian, sense, domain, data, min, calendar)
!load the values in an axistype (still need to call mpp_write)
!write metadata attributes for axis
!it is declared intent(inout) so you can nullify pointers in the incoming object if needed
!the f90 standard doesn't guarantee that intent(out) on a type guarantees that its pointer components will be unassociated
      integer,          intent(in)           :: unit
      type(axistype),   intent(inout)        :: axis
      character(len=*), intent(in)           :: name, units, longname
      character(len=*), intent(in), optional :: cartesian
      integer,          intent(in), optional :: sense
      type(domain1D),   intent(in), optional :: domain
      real,             intent(in), optional :: data(:)
      real,             intent(in), optional :: min
      character(len=*), intent(in), optional :: calendar

      integer :: is, ie, isg, ieg
      integer :: istat
      logical :: domain_exist
      type(domain2d), pointer :: io_domain => NULL()

!      call mpp_clock_begin(mpp_write_clock)
      !--- the shift and cartesian information is needed in mpp_write_meta_field from all the pe.
      !--- we may revise this in the future.
      axis%cartesian = 'N'
      if( PRESENT(cartesian) )axis%cartesian = cartesian

      domain_exist = .false.

     if( PRESENT(domain) ) then
         domain_exist = .true.
         call mpp_get_global_domain( domain, isg, ieg )
         if(mpp_file(unit)%io_domain_exist) then
            io_domain => mpp_get_io_domain(mpp_file(unit)%domain)
            if(axis%cartesian=='X') then
               call mpp_get_global_domain( io_domain, xbegin=is, xend=ie)
            else if(axis%cartesian=='Y') then
               call mpp_get_global_domain( io_domain, ybegin=is, yend=ie)
            endif
         else
            call mpp_get_compute_domain( domain, is, ie )
         endif
      else if( PRESENT(data) )then
         isg=1; ieg=size(data(:)); is=isg; ie=ieg
      endif

      axis%shift = 0
      if( PRESENT(data) .AND. domain_exist ) then
         if( size(data(:)) == ieg-isg+2 ) then
            axis%shift = 1
            ie = ie   + 1
            ieg = ieg + 1
         endif
      endif

      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe) then
!         call mpp_clock_end(mpp_write_clock)
         return
      endif
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized ) &
           call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )

!pre-existing pointers need to be nullified
      if( ASSOCIATED(axis%data) ) then
         DEALLOCATE(axis%data, stat=istat)
      endif
!load axistype
      axis%name     = name
      axis%units    = units
      axis%longname = longname
      if( PRESENT(calendar)  ) axis%calendar  = calendar
      if( PRESENT(sense)     ) axis%sense     = sense
      if( PRESENT(data) )then
         if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. domain_exist ) then
            axis%len = ie - is + 1
            allocate(axis%data(axis%len))
            axis%data = data(is-isg+1:ie-isg+1)
         else
            axis%len = size(data(:))
            allocate(axis%data(axis%len))
            axis%data = data
         endif
      endif
!write metadata
      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
#ifdef use_netCDF
!write axis def
!space axes are always floats, time axis is always double
          if( ASSOCIATED(axis%data) )then !space axis
              error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, axis%len, axis%did )
              call netcdf_err( error, mpp_file(unit), axis )
              if(pack_size == 1) then
                 error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_DOUBLE, 1, axis%did, axis%id )
              else ! pack_size == 2
                 error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_FLOAT, 1, axis%did, axis%id )
              endif
              call netcdf_err( error, mpp_file(unit), axis )
          else                            !time axis
              if( mpp_file(unit)%id.NE.-1 ) &
                   call mpp_error( FATAL, 'MPP_WRITE_META_AXIS: There is already a time axis for this file.' )
              error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, NF_UNLIMITED, axis%did )
              call netcdf_err( error, mpp_file(unit), axis )
              if(pack_size == 1) then
                 error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_DOUBLE, 1, axis%did, axis%id )
              else ! pack_size == 2
                 error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_FLOAT, 1, axis%did, axis%id )
              endif
              call netcdf_err( error, mpp_file(unit), axis )
              mpp_file(unit)%id = axis%id !file ID is the same as time axis varID
          end if
#endif
      else
          varnum = varnum + 1
          axis%id = varnum
          axis%did = varnum
!write axis def
          write( text, '(a,i4,a)' )'AXIS ', axis%id, ' name'
          call write_attribute( unit, trim(text), cval=axis%name )
          write( text, '(a,i4,a)' )'AXIS ', axis%id, ' size'
          if( ASSOCIATED(axis%data) )then !space axis
!              if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
!                  call write_attribute( unit, trim(text), ival=(/ie-is+1/) )
!              else
                  call write_attribute( unit, trim(text), ival=(/size(axis%data(:))/) )
!              end if
          else                            !time axis
              if( mpp_file(unit)%id.NE.-1 ) &
                   call mpp_error( FATAL, 'MPP_WRITE_META_AXIS: There is already a time axis for this file.' )
              call write_attribute( unit, trim(text), ival=(/0/) ) !a size of 0 indicates time axis
              mpp_file(unit)%id = axis%id
          end if
      end if
!write axis attributes
      call mpp_write_meta( unit, axis%id, 'long_name', cval=axis%longname) ; axis%natt = axis%natt + 1
      if (lowercase(trim(axis%units)).ne.'none' .OR. .NOT.cf_compliance) then
        call mpp_write_meta( unit, axis%id, 'units',     cval=axis%units) ; axis%natt = axis%natt + 1
      endif
      if( PRESENT(calendar) ) then
        if (.NOT.cf_compliance) then
          call mpp_write_meta( unit, axis%id, 'calendar', cval=axis%calendar)
        else
          call mpp_write_meta( unit, axis%id, 'calendar', cval=lowercase(axis%calendar))
        endif
        axis%natt = axis%natt + 1
      endif
      if( PRESENT(cartesian) ) then
        if (.NOT.cf_compliance) then
          call mpp_write_meta( unit, axis%id, 'cartesian_axis', cval=axis%cartesian)
          axis%natt = axis%natt + 1
        else
          if (trim(axis%cartesian).ne.'N') then
            call mpp_write_meta( unit, axis%id, 'axis', cval=axis%cartesian)
            axis%natt = axis%natt + 1
          endif
        endif
      endif
      if( PRESENT(sense) )then
          if( sense.EQ.-1 )then
              call mpp_write_meta( unit, axis%id, 'positive', cval='down')
              axis%natt = axis%natt + 1
          else if( sense.EQ.1 )then
              call mpp_write_meta( unit, axis%id, 'positive', cval='up')
              axis%natt = axis%natt + 1
          else
             ! silently ignore values of sense other than +/-1.
          end if
      end if
      if( PRESENT(min) ) then
        call mpp_write_meta( unit, axis%id, 'valid_min', rval=min)
        axis%natt = axis%natt + 1
      endif
      if( mpp_file(unit)%threading.EQ.MPP_MULTI .AND. mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. domain_exist )then
          call mpp_write_meta( unit, axis%id, 'domain_decomposition', ival=(/isg,ieg,is,ie/))
          axis%natt = axis%natt + 1
      end if
      if( verbose )print '(a,2i6,x,a,2i3)', 'MPP_WRITE_META: Wrote axis metadata, pe, unit, axis%name, axis%id, axis%did=', &
           pe, unit, trim(axis%name), axis%id, axis%did

      mpp_file(unit)%ndim = max(1,mpp_file(unit)%ndim + 1)

!      call mpp_clock_end(mpp_write_clock)
      return
    end subroutine mpp_write_meta_axis_r1d

    subroutine mpp_write_meta_axis_i1d(unit, axis, name, units, longname, data, min, compressed)
!load the values in an axistype (still need to call mpp_write)
!write metadata attributes for axis
!it is declared intent(inout) so you can nullify pointers in the incoming object if needed
!the f90 standard doesn't guarantee that intent(out) on a type guarantees that its pointer components will be unassociated
      integer,          intent(in)           :: unit
      type(axistype),   intent(inout)        :: axis
      character(len=*), intent(in)           :: name, units, longname
      integer,          intent(in)           :: data(:)
      integer,          intent(in), optional :: min
      character(len=*), intent(in), optional :: compressed

      integer :: istat
      logical :: domain_exist
      type(domain2d), pointer :: io_domain => NULL()

!      call mpp_clock_begin(mpp_write_clock)
      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META_I1D: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe) then
!         call mpp_clock_end(mpp_write_clock)
         return
      endif
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized ) &
           call mpp_error( FATAL, 'MPP_WRITE_META_I1D: cannot write metadata to file after an mpp_write.' )

!pre-existing pointers need to be nullified
      if( ASSOCIATED(axis%idata) ) then
         DEALLOCATE(axis%idata, stat=istat)
      endif
!load axistype
      axis%name     = name
      axis%units    = units
      axis%longname = longname
      if( PRESENT(compressed)) axis%compressed = trim(compressed)
      axis%len = size(data(:))
      allocate(axis%idata(axis%len))
       axis%idata = data
!write metadata
#ifdef use_netCDF
      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
          error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, axis%len, axis%did )
          call netcdf_err( error, mpp_file(unit), axis )
          error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_INT, 1, axis%did, axis%id )
          call netcdf_err( error, mpp_file(unit), axis )
      else
          call mpp_error( FATAL, 'MPP_WRITE_META_AXIS_I1D: Only netCDF format is currently supported.' )
      end if
#endif
!write axis attributes
      call mpp_write_meta( unit, axis%id, 'long_name', cval=axis%longname) ; axis%natt = axis%natt + 1
      if (lowercase(trim(axis%units)).ne.'none' .OR. .NOT.cf_compliance) then
        call mpp_write_meta( unit, axis%id, 'units',     cval=axis%units) ; axis%natt = axis%natt + 1
      endif
      if( PRESENT(compressed) ) then
        call mpp_write_meta( unit, axis%id, 'compress', cval=axis%compressed)
        axis%natt = axis%natt + 1
      endif
      if( PRESENT(min) ) then
        call mpp_write_meta( unit, axis%id, 'valid_min', ival=min)
        axis%natt = axis%natt + 1
      endif
      if( verbose )print '(a,2i6,x,a,2i3)', 'MPP_WRITE_META: Wrote axis metadata, pe, unit, axis%name, axis%id, axis%did=', &
           pe, unit, trim(axis%name), axis%id, axis%did

      mpp_file(unit)%ndim = max(1,mpp_file(unit)%ndim + 1)

!      call mpp_clock_end(mpp_write_clock)
      return
    end subroutine mpp_write_meta_axis_i1d


    subroutine mpp_write_meta_axis_unlimited(unit, axis, name, data, unlimited, units, longname)
!load the values in an axistype (still need to call mpp_write)
!write metadata attributes for axis
!it is declared intent(inout) so you can nullify pointers in the incoming object if needed
!the f90 standard doesn't guarantee that intent(out) on a type guarantees that its pointer components will be unassociated
      integer,          intent(in)           :: unit
      type(axistype),   intent(inout)        :: axis
      character(len=*), intent(in)           :: name
      integer,          intent(in)           :: data  ! Number of elements to be written
      logical,          intent(in)           :: unlimited  ! Provides unique arg signature
      character(len=*), intent(in), optional :: units, longname

      integer :: istat
      logical :: domain_exist
      type(domain2d), pointer :: io_domain => NULL()

!      call mpp_clock_begin(mpp_write_clock)
      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META_I1D: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe) then
!         call mpp_clock_end(mpp_write_clock)
         return
      endif
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized ) &
           call mpp_error( FATAL, 'MPP_WRITE_META_I1D: cannot write metadata to file after an mpp_write.' )

!load axistype
      axis%name     = name
      if(present(units)) axis%units    = units
      if(present(longname)) axis%longname = longname
      axis%len = 1
      allocate(axis%idata(1))
      axis%idata = data
!write metadata
#ifdef use_netCDF
      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
          error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, NF_UNLIMITED, axis%did )
          call netcdf_err( error, mpp_file(unit), axis )
          error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_INT, 0, axis%did, axis%id )
          call netcdf_err( error, mpp_file(unit), axis )
      else
          call mpp_error( FATAL, 'MPP_WRITE_META_AXIS_UNLIMITED: Only netCDF format is currently supported.' )
      end if
#endif
!write axis attributes
      if(present(longname)) then
         call mpp_write_meta(unit,axis%id,'long_name',cval=axis%longname); axis%natt=axis%natt+1
      endif
      if(present(units)) then
         if (lowercase(trim(axis%units)).ne.'none' .OR. .NOT.cf_compliance) then
           call mpp_write_meta(unit,axis%id,'units',     cval=axis%units); axis%natt=axis%natt+1
         endif
      endif
      if( verbose )print '(a,2i6,x,a,2i3)', &
          'MPP_WRITE_META_UNLIMITED: Wrote axis metadata, pe, unit, axis%name, axis%id, axis%did=', &
           pe, unit, trim(axis%name), axis%id, axis%did

      mpp_file(unit)%ndim = max(1,mpp_file(unit)%ndim + 1)

!      call mpp_clock_end(mpp_write_clock)
      return
    end subroutine mpp_write_meta_axis_unlimited


    subroutine mpp_write_meta_field( unit, field, axes, name, units, longname,&
         min, max, missing, fill, scale, add, pack, time_method, standard_name, checksum)
!define field: must have already called mpp_write_meta(axis) for each axis
      integer, intent(in) :: unit
      type(fieldtype), intent(inout) :: field
      type(axistype), intent(in) :: axes(:)
      character(len=*), intent(in) :: name, units, longname
      real, intent(in), optional :: min, max, missing, fill, scale, add
      integer, intent(in), optional :: pack
      character(len=*), intent(in), optional :: time_method
      character(len=*), intent(in), optional :: standard_name
      integer(LONG_KIND), dimension(:), intent(in), optional :: checksum
!this array is required because of f77 binding on netCDF interface
      integer, allocatable :: axis_id(:)
      real :: a, b
      integer :: i, istat, ishift, jshift
      character(len=64) :: checksum_char

!      call mpp_clock_begin(mpp_write_clock)

      !--- figure out the location of data, this is needed in mpp_write.
      !--- for NON-symmetry domain, the position is not an issue.
      !--- we may need to rethink how to address the symmetric issue.
      ishift = 0; jshift = 0
      do i = 1, size(axes(:))
         select case ( lowercase( axes(i)%cartesian ) )
         case ( 'x' )
            ishift = axes(i)%shift
         case ( 'y' )
            jshift = axes(i)%shift
         end select
      end do

      field%position = CENTER
      if(ishift == 1 .AND. jshift == 1) then
         field%position = CORNER
      else if(ishift == 1) then
         field%position = EAST
      else if(jshift == 1) then
         field%position = NORTH
      endif

      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )

      if( .NOT.mpp_file(unit)%write_on_this_pe) then
         if( .NOT. ASSOCIATED(field%axes) )allocate(field%axes(1)) !temporary fix
!         call mpp_clock_end(mpp_write_clock)
         return
      endif
      if( .NOT.mpp_file(unit)%opened ) call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized )  then
!     File has already been written to and needs to be returned to define mode.
#ifdef use_netCDF
        error = NF_REDEF(mpp_file(unit)%ncid)
#endif
        mpp_file(unit)%initialized = .false.
      endif
!           call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )

!pre-existing pointers need to be nullified
      if( ASSOCIATED(field%axes) ) DEALLOCATE(field%axes, stat=istat)
      if( ASSOCIATED(field%size) ) DEALLOCATE(field%size, stat=istat)
!fill in field metadata
      field%name = name
      field%units = units
      field%longname = longname
      allocate( field%axes(size(axes(:))) )
      field%axes = axes
      field%ndim = size(axes(:))
      field%time_axis_index = -1 !this value will never match any axis index
!size is buffer area for the corresponding axis info: it is required to buffer this info in the fieldtype
!because axis might be reused in different files
      allocate( field%size(size(axes(:))) )
      do i = 1,size(axes(:))
         if( ASSOCIATED(axes(i)%data) )then !space axis
             field%size(i) = size(axes(i)%data(:))
         else               !time
             field%size(i) = 1
             field%time_axis_index = i
         end if
      end do
!attributes
      if( PRESENT(min) )          field%min           = min
      if( PRESENT(max) )          field%max           = max
      if( PRESENT(scale) )        field%scale         = scale
      if( PRESENT(add) )          field%add           = add
      if( PRESENT(standard_name)) field%standard_name = standard_name
      if( PRESENT(missing) )      field%missing       = missing
      if( PRESENT(fill) )         field%fill          = fill
      field%checksum      = 0
      if( PRESENT(checksum) )     field%checksum(1:size(checksum)) = checksum(:)

      ! Issue warning if fill and missing are different
      if ( (present(fill).and.present(missing)) .and. (field%missing .ne. field%fill) ) then
         call mpp_error(WARNING, 'MPP_WRITE_META: NetCDF attributes _FillValue and missing_value should be equal.')
      end if
!pack is currently used only for netCDF
      field%pack = 2        !default write 32-bit floats
      if( PRESENT(pack) )field%pack = pack
      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
#ifdef use_netCDF
          allocate( axis_id(size(field%axes(:))) )
          do i = 1,size(field%axes(:))
             axis_id(i) = field%axes(i)%did
          end do
!write field def
          select case (field%pack)
              case(0)
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_INT, size(field%axes(:)), axis_id, field%id )
              case(1)
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_DOUBLE, size(field%axes(:)), axis_id, field%id )
              case(2)
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_FLOAT,  size(field%axes(:)), axis_id, field%id )
              case(4)
                  if( .NOT.PRESENT(scale) .OR. .NOT.PRESENT(add) ) &
                       call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: scale and add must be supplied when pack=4.' )
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_SHORT,  size(field%axes(:)), axis_id, field%id )
              case(8)
                  if( .NOT.PRESENT(scale) .OR. .NOT.PRESENT(add) ) &
                       call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: scale and add must be supplied when pack=8.' )
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_BYTE,   size(field%axes(:)), axis_id, field%id )
              case default
                  call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: only legal packing values are 1,2,4,8.' )
          end select
          call netcdf_err( error, mpp_file(unit), field=field )
          deallocate(axis_id)
#ifndef use_netCDF3
          if(shuffle .NE. 0 .OR. deflate .NE. 0) then
             error = NF_DEF_VAR_DEFLATE(mpp_file(unit)%ncid, field%id, shuffle, deflate, deflate_level)
             call netcdf_err( error, mpp_file(unit), field=field )
          endif
#endif
#endif
      else
          varnum = varnum + 1
          field%id = varnum
          if( PRESENT(pack) )call mpp_error( WARNING, 'MPP_WRITE_META: Packing is currently available only on netCDF files.' )
!write field def
          write( text, '(a,i4,a)' )'FIELD ', field%id, ' name'
          call write_attribute( unit, trim(text), cval=field%name )
          write( text, '(a,i4,a)' )'FIELD ', field%id, ' axes'
          call write_attribute( unit, trim(text), ival=field%axes(:)%did )
      end if
!write field attributes: these names follow netCDF conventions
      call mpp_write_meta( unit, field%id, 'long_name', cval=field%longname)
      if (lowercase(trim(field%units)).ne.'none' .OR. .NOT.cf_compliance) then
        call mpp_write_meta( unit, field%id, 'units',     cval=field%units)
      endif
!all real attributes must be written as packed
      if( PRESENT(min) .AND. PRESENT(max) )then
          if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
              call mpp_write_meta( unit, field%id, 'valid_range', rval=(/min,max/), pack=pack )
          else
              a = nint((min-add)/scale)
              b = nint((max-add)/scale)
              call mpp_write_meta( unit, field%id, 'valid_range', rval=(/a,  b  /), pack=pack )
          end if
      else if( PRESENT(min) )then
          if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
              call mpp_write_meta( unit, field%id, 'valid_min', rval=field%min, pack=pack )
          else
              a = nint((min-add)/scale)
              call mpp_write_meta( unit, field%id, 'valid_min', rval=a, pack=pack )
          end if
      else if( PRESENT(max) )then
          if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
              call mpp_write_meta( unit, field%id, 'valid_max', rval=field%max, pack=pack )
          else
              a = nint((max-add)/scale)
              call mpp_write_meta( unit, field%id, 'valid_max', rval=a, pack=pack )
          end if
      end if
! write missing_value
      if ( present(missing) ) then
         if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
            call mpp_write_meta( unit, field%id, 'missing_value', rval=field%missing, pack=pack )
         else
            a = nint((missing-add)/scale)
            call mpp_write_meta( unit, field%id, 'missing_value', rval=a, pack=pack )
         end if
      end if
! write _FillValue
      if ( present(fill) ) then
         if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
            call mpp_write_meta( unit, field%id, '_FillValue', rval=field%fill, pack=pack )
         else if (field%pack==0) then ! some safety checks for integer fills
            if ( present(scale).OR.present(add) ) then
               call mpp_error(FATAL,"add,scale not currently implimented for pack=0 int handling, try reals instead.")
            else
               ! Trust No One
               call mpp_write_meta( unit, field%id, '_FillValue', ival=MPP_FILL_INT, pack=pack )
            end if
         else
            a = nint((fill-add)/scale)
            call mpp_write_meta( unit, field%id, '_FillValue', rval=a, pack=pack )
         end if
      end if

      if( field%pack.NE.1 .AND. field%pack.NE.2 )then
          call mpp_write_meta( unit, field%id, 'packing', ival=field%pack )
          if( PRESENT(scale) )call mpp_write_meta( unit, field%id, 'scale_factor',  rval=field%scale )
          if( PRESENT(add)   )call mpp_write_meta( unit, field%id, 'add_offset',    rval=field%add   )
      end if

      if( present(checksum) )then
          write (checksum_char,'(Z16)') field%checksum(1)
          do i = 2,size(checksum)
            write (checksum_char,'(a,Z16)') trim(checksum_char)//",",checksum(i)
          enddo
          call mpp_write_meta( unit, field%id, 'checksum', cval=checksum_char )
      end if

      if ( PRESENT(time_method) ) then
          call mpp_write_meta(unit,field%id, 'cell_methods',cval='time: '//trim(time_method))
      endif
      if ( PRESENT(standard_name)) &
           call mpp_write_meta(unit,field%id,'standard_name ', cval=field%standard_name)

      if( verbose )print '(a,2i6,x,a,i3)', 'MPP_WRITE_META: Wrote field metadata: pe, unit, field%name, field%id=', &
           pe, unit, trim(field%name), field%id

!      call mpp_clock_end(mpp_write_clock)
      return
    end subroutine mpp_write_meta_field

    subroutine write_attribute( unit, name, rval, ival, cval, pack )
!called to write metadata for non-netCDF I/O
      integer, intent(in) :: unit
      character(len=*), intent(in) :: name
      real, intent(in), optional :: rval(:)
      integer, intent(in), optional :: ival(:)
      character(len=*), intent(in), optional :: cval
!pack is currently ignored in this routine: only used by netCDF I/O
      integer, intent(in), optional :: pack

      if( mpp_file(unit)%nohdrs )return
!encode text string
      if( PRESENT(rval) )then
          write( text,* )trim(name)//'=', rval
      else if( PRESENT(ival) )then
          write( text,* )trim(name)//'=', ival
      else if( PRESENT(cval) )then
          text = ' '//trim(name)//'='//trim(cval)
      else
          call mpp_error( FATAL, 'WRITE_ATTRIBUTE: one of rval, ival, cval must be present.' )
      end if
      if( mpp_file(unit)%format.EQ.MPP_ASCII )then
!implies sequential access
          write( unit,fmt='(a)' )trim(text)//char(10)
      else                      !MPP_IEEE32 or MPP_NATIVE
          if( mpp_file(unit)%access.EQ.MPP_SEQUENTIAL )then
              write(unit)trim(text)//char(10)
          else                  !MPP_DIRECT
              write( unit,rec=mpp_file(unit)%record )trim(text)//char(10)
              if( verbose )print '(a,i6,a,i3)', 'WRITE_ATTRIBUTE: PE=', pe, ' wrote record ', mpp_file(unit)%record
              mpp_file(unit)%record = mpp_file(unit)%record + 1
          end if
      end if
      return
    end subroutine write_attribute

    subroutine write_attribute_netcdf( unit, id, name, rval, ival, cval, pack )
!called to write metadata for netCDF I/O
      integer, intent(in) :: unit
      integer, intent(in) :: id
      character(len=*), intent(in) :: name
      real,             intent(in), optional :: rval(:)
      integer,          intent(in), optional :: ival(:)
      character(len=*), intent(in), optional :: cval
      integer, intent(in), optional :: pack
      integer, allocatable :: rval_i(:)
#ifdef use_netCDF
      if( PRESENT(rval) )then
!pack was only meaningful for FP numbers, but is now extended by the ival branch of this routine
          if( PRESENT(pack) )then
             if( pack== 0 ) then !! here be dragons, use ival branch!...
                if( KIND(rval).EQ.DOUBLE_KIND )then
                   call mpp_error( FATAL, &
                        'WRITE_ATTRIBUTE_NETCDF: attempting to write internal NF_INT, currently int32, as double.' )
                   error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_DOUBLE, size(rval(:)), rval )
                else if( KIND(rval).EQ.FLOAT_KIND )then
                   call mpp_error( FATAL, &
                        'WRITE_ATTRIBUTE_NETCDF: attempting to write internal NF_INT, currently int32, as float.' )
                   error = NF_PUT_ATT_REAL  ( mpp_file(unit)%ncid, id, name, NF_DOUBLE, size(rval(:)), rval )
                end if
                call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
             else if( pack.EQ.1 )then
                  if( KIND(rval).EQ.DOUBLE_KIND )then
                      error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_DOUBLE, size(rval(:)), rval )
                  else if( KIND(rval).EQ.FLOAT_KIND )then
                      call mpp_error( WARNING, &
                           'WRITE_ATTRIBUTE_NETCDF: attempting to write internal 32-bit real as external 64-bit.' )
                      error = NF_PUT_ATT_REAL  ( mpp_file(unit)%ncid, id, name, NF_DOUBLE, size(rval(:)), rval )
                  end if
                  call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
              else if( pack.EQ.2 )then
                  if( KIND(rval).EQ.DOUBLE_KIND )then
                      error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_FLOAT,  size(rval(:)), rval )
                  else if( KIND(rval).EQ.FLOAT_KIND )then
                      error = NF_PUT_ATT_REAL  ( mpp_file(unit)%ncid, id, name, NF_FLOAT,  size(rval(:)), rval )
                  end if
                  call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
              else if( pack.EQ.4 )then
                  allocate( rval_i(size(rval(:))) )
                  rval_i = rval
                  if( KIND(rval).EQ.DOUBLE_KIND )then
                      error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_SHORT,  size(rval_i(:)), rval )
                  else if( KIND(rval).EQ.FLOAT_KIND )then
                      error = NF_PUT_ATT_REAL  ( mpp_file(unit)%ncid, id, name, NF_SHORT,  size(rval_i(:)), rval )
                  end if
                  call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
                  deallocate(rval_i)
              else if( pack.EQ.8 )then
                  allocate( rval_i(size(rval(:))) )
                  rval_i = rval
                  if( KIND(rval).EQ.DOUBLE_KIND )then
                      error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_BYTE,   size(rval_i(:)), rval )
                  else if( KIND(rval).EQ.FLOAT_KIND )then
                      error = NF_PUT_ATT_REAL  ( mpp_file(unit)%ncid, id, name, NF_BYTE,   size(rval_i(:)), rval )
                  end if
                  call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
                  deallocate(rval_i)
              else
                  call mpp_error( FATAL, 'WRITE_ATTRIBUTE_NETCDF: only legal packing values are 1,2,4,8.' )
              end if
          else
!default is to write FLOATs (32-bit)
              if( KIND(rval).EQ.DOUBLE_KIND )then
                  error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_FLOAT,  size(rval(:)), rval )
              else if( KIND(rval).EQ.FLOAT_KIND )then
                  error = NF_PUT_ATT_REAL  ( mpp_file(unit)%ncid, id, name, NF_FLOAT,  size(rval(:)), rval )
              end if
              call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
          end if
      else if( PRESENT(ival) )then
          if( PRESENT(pack) ) then
             if (pack ==0) then
                if (KIND(ival).EQ.LONG_KIND ) then
                   call mpp_error(FATAL,'only use NF_INTs with pack=0 for now')
                end if
                error = NF_PUT_ATT_INT( mpp_file(unit)%ncid, id, name, NF_INT, size(ival(:)), ival ) !!XXX int32_t..
                call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
             else
                call mpp_error( FATAL, 'WRITE_ATTRIBUTE_NETCDF: only implimented ints when pack=0, else use reals.' )
             endif
          else
          error = NF_PUT_ATT_INT ( mpp_file(unit)%ncid, id, name, NF_INT, size(ival(:)), ival )
          call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
          end if
      else if( present(cval) )then
          if (.NOT.cf_compliance .or. trim(name).NE.'calendar') then
            error = NF_PUT_ATT_TEXT( mpp_file(unit)%ncid, id, name, len_trim(cval), cval )
          else
            error = NF_PUT_ATT_TEXT( mpp_file(unit)%ncid, id, name, len_trim(cval), lowercase(cval) )
          endif
          call netcdf_err( error, mpp_file(unit), string=' Attribute='//name )
      else
          call mpp_error( FATAL, 'WRITE_ATTRIBUTE_NETCDF: one of rval, ival, cval must be present.' )
      end if
#endif /* use_netCDF */
      return
    end subroutine write_attribute_netcdf

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                      !
!                             MPP_WRITE                                !
!                                                                      !
! mpp_write is used to write data to the file on <unit> using the      !
! file parameters supplied by mpp_open(). Axis and field definitions   !
! must have previously been written to the file using mpp_write_meta.  !
!                                                                      !
! mpp_write can take 2 forms, one for distributed data and one for     !
! non-distributed data. Distributed data refer to arrays whose two     !
! fastest-varying indices are domain-decomposed. Distributed data      !
! must be 2D or 3D (in space). Non-distributed data can be 0-3D.       !
!                                                                      !
! In all calls to mpp_write, tstamp is an optional argument. It is to  !
! be omitted if the field was defined not to be a function of time.    !
! Results are unpredictable if the argument is supplied for a time-    !
! independent field, or omitted for a time-dependent field. Repeated   !
! writes of a time-independent field are also not recommended. One     !
! time level of one field is written per call.                         !
!                                                                      !
!                                                                      !
! For non-distributed data, use                                        !
!                                                                      !
!  mpp_write( unit, field, data, tstamp )                              !
!     integer, intent(in) :: unit                                      !
!     type(fieldtype), intent(in) :: field                             !
!     real(DOUBLE_KIND), optional :: tstamp                            !
!     data is real and can be scalar or of rank 1-3.                   !
!                                                                      !
! For distributed data, use                                            !
!                                                                      !
!  mpp_write( unit, field, domain, data, tstamp )                      !
!     integer, intent(in) :: unit                                      !
!     type(fieldtype), intent(in) :: field                             !
!     type(domain2D), intent(in) :: domain                             !
!     real(DOUBLE_KIND), optional :: tstamp                            !
!     data is real and can be of rank 2 or 3.                          !
!                                                                      !
!  mpp_write( unit, axis )                                             !
!     integer, intent(in) :: unit                                      !
!     type(axistype), intent(in) :: axis                               !
!                                                                      !
! This call writes the actual co-ordinate values along each space      !
! axis. It must be called once for each space axis after all other     !
! metadata has been written.                                           !
!                                                                      !
! The mpp_write package also includes the routine write_record which   !
! performs the actual write. This routine is private to this module.   !
!                                                                      !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#undef WRITE_RECORD_
#define WRITE_RECORD_ write_record_default
#undef MPP_WRITE_2DDECOMP_2D_
#define MPP_WRITE_2DDECOMP_2D_ mpp_write_2ddecomp_r2d
#undef MPP_WRITE_2DDECOMP_3D_
#define MPP_WRITE_2DDECOMP_3D_ mpp_write_2ddecomp_r3d
#undef MPP_WRITE_2DDECOMP_4D_
#define MPP_WRITE_2DDECOMP_4D_ mpp_write_2ddecomp_r4d
#undef MPP_TYPE_
#define MPP_TYPE_ real
#include <mpp_write_2Ddecomp.h>

#ifdef OVERLOAD_R8
#define WRITE_RECORD_ write_record_r8
#undef MPP_WRITE_2DDECOMP_2D_
#undef MPP_WRITE_2DDECOMP_2D_
#define MPP_WRITE_2DDECOMP_2D_ mpp_write_2ddecomp_r2d_r8
#undef MPP_WRITE_2DDECOMP_3D_
#define MPP_WRITE_2DDECOMP_3D_ mpp_write_2ddecomp_r3d_r8
#undef MPP_WRITE_2DDECOMP_4D_
#define MPP_WRITE_2DDECOMP_4D_ mpp_write_2ddecomp_r4d_r8
#undef MPP_TYPE_
#define MPP_TYPE_ real(DOUBLE_KIND)
#include <mpp_write_2Ddecomp.h>
#endif

#undef MPP_WRITE_COMPRESSED_1D_
#define MPP_WRITE_COMPRESSED_1D_ mpp_write_compressed_r1d
#undef MPP_WRITE_COMPRESSED_2D_
#define MPP_WRITE_COMPRESSED_2D_ mpp_write_compressed_r2d
#undef MPP_WRITE_COMPRESSED_3D_
#define MPP_WRITE_COMPRESSED_3D_ mpp_write_compressed_r3d
#undef MPP_TYPE_
#define MPP_TYPE_ real
#include <mpp_write_compressed.h>

#undef MPP_WRITE_UNLIMITED_AXIS_1D_
#define MPP_WRITE_UNLIMITED_AXIS_1D_ mpp_write_unlimited_axis_r1d
#undef MPP_TYPE_
#define MPP_TYPE_ real
#include <mpp_write_unlimited_axis.h>

#undef MPP_WRITE_
#define MPP_WRITE_ mpp_write_r0D
#undef MPP_TYPE_
#define MPP_TYPE_ real
#undef MPP_RANK_
#define MPP_RANK_ !
#undef MPP_WRITE_RECORD_
#define MPP_WRITE_RECORD_ call write_record_default( unit, field, 1, (/data/), tstamp)
#include <mpp_write.h>

#undef MPP_WRITE_
#define MPP_WRITE_ mpp_write_r1D
#undef MPP_TYPE_
#define MPP_TYPE_ real
#undef MPP_WRITE_RECORD_
#define MPP_WRITE_RECORD_ call write_record_default( unit, field, size(data(:)), data, tstamp)
#undef MPP_RANK_
#define MPP_RANK_ (:)
#include <mpp_write.h>

#undef MPP_WRITE_
#define MPP_WRITE_ mpp_write_r2D
#undef MPP_TYPE_
#define MPP_TYPE_ real
#undef MPP_WRITE_RECORD_
#define MPP_WRITE_RECORD_ call write_record_default( unit, field, size(data(:,:)), data, tstamp )
#undef MPP_RANK_
#define MPP_RANK_ (:,:)
#include <mpp_write.h>

#undef MPP_WRITE_
#define MPP_WRITE_ mpp_write_r3D
#undef MPP_TYPE_
#define MPP_TYPE_ real
#undef MPP_WRITE_RECORD_
#define MPP_WRITE_RECORD_ call write_record_default( unit, field, size(data(:,:,:)), data, tstamp)
#undef MPP_RANK_
#define MPP_RANK_ (:,:,:)
#include <mpp_write.h>

#undef MPP_WRITE_
#define MPP_WRITE_ mpp_write_r4D
#undef MPP_TYPE_
#define MPP_TYPE_ real
#undef MPP_WRITE_RECORD_
#define MPP_WRITE_RECORD_ call write_record_default( unit, field, size(data(:,:,:,:)), data, tstamp)
#undef MPP_RANK_
#define MPP_RANK_ (:,:,:,:)
#include <mpp_write.h>

    subroutine mpp_write_axis( unit, axis )
      integer, intent(in) :: unit
      type(axistype), intent(in) :: axis
      type(fieldtype) :: field

      call mpp_clock_begin(mpp_write_clock)
      if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe ) then
         call mpp_clock_end(mpp_write_clock)
         return
      endif
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE: invalid unit number.' )
!we convert axis to type(fieldtype) in order to call write_record
      field = default_field
      allocate( field%axes(1) )
      field%axes(1) = axis
      allocate( field%size(1) )
      field%size(1) = axis%len
      field%id = axis%id

      field%name = axis%name
      field%longname = axis%longname
      field%units = axis%units

      if(ASSOCIATED(axis%data))then
        allocate( field%axes(1)%data(size(axis%data) ))
        field%axes(1)%data = axis%data
        call write_record( unit, field, axis%len, axis%data )
      elseif(ASSOCIATED(axis%idata))then
        allocate( field%axes(1)%data(size(axis%idata) ))
        field%axes(1)%data = REAL(axis%idata)
        field%pack=4
        call write_record( unit, field, axis%len, REAL(axis%idata) )
      else
        call mpp_error( FATAL, 'MPP_WRITE_AXIS: No data associated with axis.' )
      endif

      deallocate(field%axes(1)%data)
      deallocate(field%axes,field%size)

      call mpp_clock_end(mpp_write_clock)
      return
    end subroutine mpp_write_axis

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                      !
!                          MPP_COPY_META                               !
!                                                                      !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine mpp_copy_meta_global( unit, gatt )
!writes a global metadata attribute to unit <unit>
!attribute <name> can be an real, integer or character
!one and only one of rval, ival, and cval should be present
!the first found will be used
!for a non-netCDF file, it is encoded into a string "GLOBAL <name> <val>"
      integer, intent(in) :: unit
      type(atttype), intent(in) :: gatt
      integer :: len, error

      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe )return
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized ) then
!     File has already been written to and needs to be returned to define mode.
#ifdef use_netCDF
        error = NF_REDEF(mpp_file(unit)%ncid)
#endif
        mpp_file(unit)%initialized = .false.
      endif
!           call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )
#ifdef use_netCDF
      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
         if( gatt%type.EQ.NF_CHAR )then
            len = gatt%len
            call write_attribute_netcdf( unit, NF_GLOBAL, gatt%name, cval=gatt%catt(1:len) )
         else
            call write_attribute_netcdf( unit, NF_GLOBAL, gatt%name, rval=gatt%fatt )
         endif
      else
         if( gatt%type.EQ.NF_CHAR )then
            len=gatt%len
            call write_attribute( unit, 'GLOBAL '//trim(gatt%name), cval=gatt%catt(1:len) )
         else
            call write_attribute( unit, 'GLOBAL '//trim(gatt%name), rval=gatt%fatt )
         endif
     end if
#else
     call mpp_error( FATAL, 'MPP_READ currently requires use_netCDF option' )
#endif
      return
    end subroutine mpp_copy_meta_global

    subroutine mpp_copy_meta_axis( unit, axis, domain )
!load the values in an axistype (still need to call mpp_write)
!write metadata attributes for axis.  axis is declared inout
!because the variable and dimension ids are altered

      integer, intent(in) :: unit
      type(axistype), intent(inout) :: axis
      type(domain1D), intent(in), optional :: domain
      character(len=512) :: text
      integer :: i, len, is, ie, isg, ieg, error

!      call mpp_clock_begin(mpp_write_clock)
      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe ) then
!         call mpp_clock_end(mpp_write_clock)
         return
      endif
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized )  then
!     File has already been written to and needs to be returned to define mode.
#ifdef use_netCDF
        error = NF_REDEF(mpp_file(unit)%ncid)
#endif
        mpp_file(unit)%initialized = .false.
      endif
!           call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )

! redefine domain if present
      if( PRESENT(domain) )then
          axis%domain = domain
      else
          axis%domain = NULL_DOMAIN1D
      end if

#ifdef use_netCDF
!write metadata
      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then

!write axis def
          if( ASSOCIATED(axis%data) )then !space axis
              if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
                  call mpp_get_compute_domain( axis%domain, is, ie )
                  call mpp_get_global_domain( axis%domain, isg, ieg )
                  ie  = ie + axis%shift
                  ieg = ieg + axis%shift
                  error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, ie-is+1, axis%did )
              else
                  error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, size(axis%data(:)),          axis%did )
              end if
              call netcdf_err( error, mpp_file(unit), axis )
              error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_FLOAT, 1, axis%did, axis%id )
              call netcdf_err( error, mpp_file(unit), axis )
          else                            !time axis
              error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, NF_UNLIMITED, axis%did )
              call netcdf_err( error, mpp_file(unit), axis )
              error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_DOUBLE, 1, axis%did, axis%id )
              call netcdf_err( error, mpp_file(unit), axis )
              mpp_file(unit)%id = axis%id !file ID is the same as time axis varID
              mpp_file(unit)%recdimid = axis%did ! record dimension id
          end if
      else
          varnum = varnum + 1
          axis%id = varnum
          axis%did = varnum
!write axis def
          write( text, '(a,i4,a)' )'AXIS ', axis%id, ' name'
          call write_attribute( unit, trim(text), cval=axis%name )
          write( text, '(a,i4,a)' )'AXIS ', axis%id, ' size'
          if( ASSOCIATED(axis%data) )then !space axis
              if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
                  call mpp_get_compute_domain(axis%domain, is, ie)
                  call write_attribute( unit, trim(text), ival=(/ie-is+1/) ) ! ??? is, ie is not initialized
              else
                  call write_attribute( unit, trim(text), ival=(/size(axis%data(:))/) )
              end if
          else                            !time axis
              if( mpp_file(unit)%id.NE.-1 ) &
                   call mpp_error( FATAL, 'MPP_WRITE_META_AXIS: There is already a time axis for this file.' )
              call write_attribute( unit, trim(text), ival=(/0/) ) !a size of 0 indicates time axis
              mpp_file(unit)%id = axis%id
          end if
      end if
!write axis attributes

      do i=1,axis%natt
         if( axis%Att(i)%name.NE.default_att%name )then
            if( axis%Att(i)%type.EQ.NF_CHAR )then
               len = axis%Att(i)%len
               call mpp_write_meta( unit, axis%id, axis%Att(i)%name, cval=axis%Att(i)%catt(1:len) )
            else
               call mpp_write_meta( unit, axis%id, axis%Att(i)%name, rval=axis%Att(i)%fatt)
            endif
         endif
      enddo

      if( mpp_file(unit)%threading.EQ.MPP_MULTI .AND. mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
          call mpp_write_meta( unit, axis%id, 'domain_decomposition', ival=(/isg,ieg,is,ie/) )
      end if
      if( verbose )print '(a,2i6,x,a,2i3)', 'MPP_WRITE_META: Wrote axis metadata, pe, unit, axis%name, axis%id, axis%did=', &
           pe, unit, trim(axis%name), axis%id, axis%did
#else
      call mpp_error( FATAL, 'MPP_READ currently requires use_netCDF option' )
#endif
!      call mpp_clock_end(mpp_write_clock)
      return
    end subroutine mpp_copy_meta_axis

    subroutine mpp_copy_meta_field( unit, field, axes )
!useful for copying field metadata from a previous call to mpp_read_meta
!define field: must have already called mpp_write_meta(axis) for each axis
      integer, intent(in) :: unit
      type(fieldtype), intent(inout) :: field
      type(axistype), intent(in), optional :: axes(:)
!this array is required because of f77 binding on netCDF interface
      integer, allocatable :: axis_id(:)
      real :: a, b
      integer :: i, error

!      call mpp_clock_begin(mpp_write_clock)
      if( .NOT.module_is_initialized    )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
      if( .NOT. mpp_file(unit)%write_on_this_pe ) then
!         call mpp_clock_end(mpp_write_clock)
         return
      endif
      if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
      if( mpp_file(unit)%initialized )  then
!     File has already been written to and needs to be returned to define mode.
#ifdef use_netCDF
        error = NF_REDEF(mpp_file(unit)%ncid)
#endif
        mpp_file(unit)%initialized = .false.
      endif
!           call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )

       if( field%pack.NE.1 .AND. field%pack.NE.2 )then
            if( field%pack.NE.4 .AND. field%pack.NE.8 ) &
               call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: only legal packing values are 1,2,4,8.' )
      end if

      if (PRESENT(axes)) then
         deallocate(field%axes)
         deallocate(field%size)
         allocate(field%axes(size(axes(:))))
         allocate(field%size(size(axes(:))))
         field%axes = axes
         do i=1,size(axes(:))
            if (ASSOCIATED(axes(i)%data)) then
               field%size(i) = size(axes(i)%data(:))
            else
               field%size(i) = 1
               field%time_axis_index = i
            endif
         enddo
      endif

      if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
#ifdef use_netCDF
          allocate( axis_id(size(field%axes(:))) )
          do i = 1,size(field%axes(:))
             axis_id(i) = field%axes(i)%did
          end do
!write field def
          select case (field%pack)
              case(1)
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_DOUBLE, size(field%axes(:)), axis_id, field%id )
              case(2)
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_FLOAT,  size(field%axes(:)), axis_id, field%id )
              case(4)
!                 if( field%scale.EQ.default_field%scale .OR. field%add.EQ.default_field%add ) &
!                      call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: scale and add must be supplied when pack=4.' )
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_SHORT,  size(field%axes(:)), axis_id, field%id )
              case(8)
!                 if( field%scale.EQ.default_field%scale .OR. field%add.EQ.default_field%add ) &
!                      call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: scale and add must be supplied when pack=8.' )
                  error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_BYTE,   size(field%axes(:)), axis_id, field%id )
              case default
                  call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: only legal packing values are 1,2,4,8.' )
          end select
          deallocate( axis_id )
#endif
      else
          varnum = varnum + 1
          field%id = varnum
          if( field%pack.NE.default_field%pack ) &
           call mpp_error( WARNING, 'MPP_WRITE_META: Packing is currently available only on netCDF files.' )
!write field def
          write( text, '(a,i4,a)' )'FIELD ', field%id, ' name'
          call write_attribute( unit, trim(text), cval=field%name )
          write( text, '(a,i4,a)' )'FIELD ', field%id, ' axes'
          call write_attribute( unit, trim(text), ival=field%axes(:)%did )
      end if
!write field attributes: these names follow netCDF conventions
      call mpp_write_meta( unit, field%id, 'long_name', cval=field%longname )
      if (lowercase(trim(field%units)).ne.'none' .OR. .NOT.cf_compliance) then
        call mpp_write_meta( unit, field%id, 'units',     cval=field%units    )
      endif
!all real attributes must be written as packed
      if( (field%min.NE.default_field%min) .AND. (field%max.NE.default_field%max) )then
          if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
              call mpp_write_meta( unit, field%id, 'valid_range', rval=(/field%min,field%max/), pack=field%pack )
          else
              a = nint((field%min-field%add)/field%scale)
              b = nint((field%max-field%add)/field%scale)
              call mpp_write_meta( unit, field%id, 'valid_range', rval=(/a,  b  /), pack=field%pack )
          end if
      else if( field%min.NE.default_field%min )then
          if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
              call mpp_write_meta( unit, field%id, 'valid_min', rval=field%min, pack=field%pack )
          else
              a = nint((field%min-field%add)/field%scale)
              call mpp_write_meta( unit, field%id, 'valid_min', rval=a, pack=field%pack )
          end if
      else if( field%max.NE.default_field%max )then
          if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
              call mpp_write_meta( unit, field%id, 'valid_max', rval=field%max, pack=field%pack )
          else
              a = nint((field%max-field%add)/field%scale)
              call mpp_write_meta( unit, field%id, 'valid_max', rval=a, pack=field%pack )
          end if
      end if
      if( field%missing.NE.default_field%missing )then
          if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
              call mpp_write_meta( unit, field%id, 'missing_value', rval=field%missing, pack=field%pack )
          else
              a = nint((field%missing-field%add)/field%scale)
              call mpp_write_meta( unit, field%id, 'missing_value', rval=a, pack=field%pack )
          end if
      end if
      if( field%fill.NE.default_field%fill )then
          if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
              call mpp_write_meta( unit, field%id, '_FillValue', rval=field%missing, pack=field%pack )
          else
              a = nint((field%fill-field%add)/field%scale)
              call mpp_write_meta( unit, field%id, '_FillValue', rval=a, pack=field%pack )
          end if
      end if
      if( field%pack.NE.1 .AND. field%pack.NE.2 )then
          call mpp_write_meta( unit, field%id, 'packing', ival=field%pack )
          if( field%scale.NE.default_field%scale )call mpp_write_meta( unit, field%id, 'scale_factor',  rval=field%scale )
          if( field%add.NE.default_field%add   )call mpp_write_meta( unit, field%id, 'add_offset',    rval=field%add   )
      end if
      if( verbose )print '(a,2i6,x,a,i3)', 'MPP_WRITE_META: Wrote field metadata: pe, unit, field%name, field%id=', &
           pe, unit, trim(field%name), field%id

!      call mpp_clock_end(mpp_write_clock)
      return
    end subroutine mpp_copy_meta_field

    subroutine mpp_modify_axis_meta( axis, name, units, longname, cartesian, data )

      type(axistype), intent(inout) :: axis
      character(len=*), intent(in), optional :: name, units, longname, cartesian
      real, dimension(:), intent(in), optional :: data

      if (PRESENT(name)) axis%name = trim(name)
      if (PRESENT(units)) axis%units = trim(units)
      if (PRESENT(longname)) axis%longname = trim(longname)
      if (PRESENT(cartesian)) axis%cartesian = trim(cartesian)
      if (PRESENT(data)) then
         axis%len = size(data(:))
         if (ASSOCIATED(axis%data)) deallocate(axis%data)
         allocate(axis%data(axis%len))
         axis%data = data
      endif

      return
    end subroutine mpp_modify_axis_meta

    subroutine mpp_modify_field_meta( field, name, units, longname, min, max, missing, axes )

      type(fieldtype), intent(inout) :: field
      character(len=*), intent(in), optional :: name, units, longname
      real, intent(in), optional :: min, max, missing
      type(axistype), dimension(:), intent(inout), optional :: axes

      if (PRESENT(name)) field%name = trim(name)
      if (PRESENT(units)) field%units = trim(units)
      if (PRESENT(longname)) field%longname = trim(longname)
      if (PRESENT(min)) field%min = min
      if (PRESENT(max)) field%max = max
      if (PRESENT(missing)) field%missing = missing
!      if (PRESENT(axes)) then
!         axis%len = size(data(:))
!         deallocate(axis%data)
!         allocate(axis%data(axis%len))
!         axis%data = data
!      endif

      return
    end subroutine mpp_modify_field_meta