module ncdio_pio

  !-----------------------------------------------------------------------
  !BOP
  !
  ! !MODULE: ncdio_pioMod
  !
  ! !DESCRIPTION:
  ! Generic interfaces to write fields to netcdf files for CLM
  !
  ! !USES:
  use shr_kind_mod   , only : r8 => shr_kind_r8, i4=>shr_kind_i4, shr_kind_cl
  use shr_infnan_mod , only : nan => shr_infnan_nan,  isnan => shr_infnan_isnan, assignment(=)
  use shr_sys_mod    , only : shr_sys_abort
  use shr_file_mod   , only : shr_file_getunit, shr_file_freeunit
  use shr_string_mod , only : shr_string_toUpper
  use shr_log_mod    , only : errMsg => shr_log_errMsg
  use spmdMod        , only : masterproc, mpicom, iam, npes
  use spmdMod        , only : MPI_REAL8, MPI_INTEGER, MPI_LOGICAL
  use elm_varcon     , only : grlnd, nameg, namet, namel, namec, namep
  use elm_varcon     , only : spval,ispval
  use elm_varctl     , only : single_column, scm_multcols, scm_nx, scm_ny, iulog
  use shr_sys_mod    , only : shr_sys_flush
  use decompMod      , only : get_elmlevel_gsize,get_elmlevel_gsmap
  use perf_mod       , only : t_startf, t_stopf
  use fileutils      , only : getavu, relavu
  use mct_mod        , only : mct_gsMap, mct_gsMap_lsize, mct_gsMap_gsize, mct_gsMap_orderedPoints
  use pio            , only : file_desc_t, io_desc_t, iosystem_desc_t, pio_64bit_offset
  use pio            , only : pio_bcast_error, pio_char, pio_clobber, pio_closefile, pio_createfile, pio_def_dim
  use pio            , only : pio_def_var, pio_double, pio_enddef, pio_get_att, pio_get_var, pio_global, pio_initdecomp
  use pio            , only : pio_inq_att, pio_inq_dimid, pio_inq_dimlen, pio_inq_dimname, pio_inq_vardimid, pio_inq_varid
  use pio            , only : pio_inq_varname, pio_inq_varndims, pio_inquire, pio_int, pio_internal_error
  use pio            , only : pio_noclobber, pio_noerr, pio_nofill, pio_nowrite, pio_offset_kind, pio_openfile
  use pio            , only : pio_put_att, pio_put_var, pio_read_darray, pio_real, pio_seterrorhandling
  use pio            , only : pio_setframe, pio_unlimited, pio_write, pio_write_darray, var_desc_t, pio_iotask_rank
  !
  ! !PUBLIC TYPES:
  implicit none
  private
  save
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: check_var          ! determine if variable is on netcdf file
  public :: check_att          ! check if attribute is on file
  public :: check_dim          ! validity check on dimension
  public :: ncd_pio_openfile   ! open a file
  public :: ncd_pio_createfile ! create a new file
  public :: ncd_pio_closefile  ! close a file
  public :: ncd_pio_init       ! called from clm_comp
  public :: ncd_enddef         ! end define mode
  public :: ncd_putatt         ! put attribute
  public :: ncd_getatt         ! get attribute
  public :: ncd_defdim         ! define dimension
  public :: ncd_inqdid         ! inquire dimension id
  public :: ncd_inqdname       ! inquire dimension name
  public :: ncd_inqdlen        ! inquire dimension length
  public :: ncd_inqfdims       ! inquire file dimnesions 
  public :: ncd_defvar         ! define variables
  public :: ncd_inqvid         ! inquire variable id
  public :: ncd_inqvname       ! inquire variable name
  public :: ncd_inqvdims       ! inquire variable ndims
  public :: ncd_inqvdids       ! inquire variable dimids
  public :: ncd_inqvdlen       ! inquire variable dimension size
  public :: ncd_inqvdname      ! inquire variable dimension name
  public :: ncd_io             ! write local data

  integer,parameter,public :: ncd_int       = pio_int
  integer,parameter,public :: ncd_log       =-pio_int
  integer,parameter,public :: ncd_float     = pio_real
  integer,parameter,public :: ncd_double    = pio_double
  integer,parameter,public :: ncd_char      = pio_char
  integer,parameter,public :: ncd_global    = pio_global
  integer,parameter,public :: ncd_write     = pio_write
  integer,parameter,public :: ncd_nowrite   = pio_nowrite
  integer,parameter,public :: ncd_clobber   = pio_clobber
  integer,parameter,public :: ncd_noclobber = pio_noclobber
  integer,parameter,public :: ncd_nofill    = pio_nofill
  integer,parameter,public :: ncd_unlimited = pio_unlimited

  ! PIO types needed for ncdio_pio interface calls
  public file_desc_t
  public var_desc_t

  !
  ! !PRIVATE MEMBER FUNCTIONS:
  !
  interface ncd_defvar
     module procedure ncd_defvar_bynf
     module procedure ncd_defvar_bygrid
  end interface

  interface ncd_putatt
     module procedure ncd_putatt_int
     module procedure ncd_putatt_real
     module procedure ncd_putatt_char
  end interface

  interface ncd_getatt
     module procedure ncd_getatt_char
  end interface ncd_getatt

  interface ncd_io 
     module procedure ncd_io_char_var0_start_glob

     !DIMS 0,1
     module procedure ncd_io_{DIMS}d_log_glob
     
     !DIMS 0,1,2,3,4  
     !TYPE int,double
     module procedure ncd_io_{DIMS}d_{TYPE}_glob

     !DIMS 0,1,2
     !TYPE text
     module procedure ncd_io_{DIMS}d_{TYPE}_glob

     !TYPE int,double
     !DIMS 1,2,3,4
     module procedure ncd_io_{DIMS}d_{TYPE}

     !TYPE logical
     !DIMS 1
     module procedure ncd_io_{DIMS}d_{TYPE}
  end interface

  interface ncd_inqvdlen
     module procedure ncd_inqvdlen_byDesc
     module procedure ncd_inqvdlen_byName
  end interface

  interface ncd_inqvdname
     module procedure ncd_inqvdname_byDesc
     module procedure ncd_inqvdname_byName
  end interface ncd_inqvdname

  private :: ncd_getiodesc      ! obtain iodesc
  private :: scam_field_offsets ! get offset to proper lat/lon gridcell for SCAM

  integer,parameter,private :: debug = 0             ! local debug level

  integer , parameter  , public  :: max_string_len = 256     ! length of strings
  real(r8), parameter  , public  :: fillvalue = 1.e36_r8     ! fill value for netcdf fields

  integer, public :: io_type

  type(iosystem_desc_t), pointer, public  :: pio_subsystem

  type iodesc_plus_type
     character(len=64) :: name
     type(IO_desc_t)   :: iodesc
     integer           :: type
     integer           :: ndims
     integer           :: dims(5)
     integer           :: dimids(5) 
  end type iodesc_plus_type
  integer,parameter      ,private :: max_iodesc = 100
  integer                ,private :: num_iodesc = 0
  type(iodesc_plus_type) ,private, target :: iodesc_list(max_iodesc)
  !-----------------------------------------------------------------------

contains

  !-----------------------------------------------------------------------
  subroutine ncd_pio_init()
    !
    ! !DESCRIPTION:
    ! Initial PIO
    !
    ! !USES:
    use shr_pio_mod     , only : shr_pio_getiosys, shr_pio_getiotype
    use elm_varctl      , only : inst_name
    !-----------------------------------------------------------------------

    PIO_subsystem => shr_pio_getiosys(inst_name)
    io_type       =  shr_pio_getiotype(inst_name)

  end subroutine ncd_pio_init

  !-----------------------------------------------------------------------
  subroutine ncd_pio_openfile(file, fname, mode)
    !
    ! !DESCRIPTION:
    ! Open a NetCDF PIO file
    !
    ! !ARGUMENTS:
    class(file_desc_t) , intent(inout) :: file   ! Output PIO file handle
    character(len=*)   , intent(in)    :: fname  ! Input filename to open
    integer            , intent(in)    :: mode   ! file mode
    !
    ! !LOCAL VARIABLES:
    integer :: ierr
    !-----------------------------------------------------------------------

    ierr = pio_openfile(pio_subsystem, file, io_type, fname, mode)

    if(ierr/= PIO_NOERR) then
       call shr_sys_abort('ncd_pio_openfile ERROR: Failed to open file')
    else if(pio_iotask_rank(pio_subsystem)==0) then
       write(iulog,*) 'Opened existing file ', trim(fname), file%fh
    end if

  end subroutine ncd_pio_openfile

  !-----------------------------------------------------------------------
  subroutine ncd_pio_closefile(file)
    !
    ! !DESCRIPTION:
    ! Close a NetCDF PIO file
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout) :: file   ! PIO file handle to close
    !-----------------------------------------------------------------------

    call pio_closefile(file)

  end subroutine ncd_pio_closefile

  !-----------------------------------------------------------------------
  subroutine ncd_pio_createfile(file, fname, avoid_pnetcdf)
    !
    ! !DESCRIPTION:
    ! Create a new NetCDF file with PIO
    !
    ! !USES:
    use pio, only : pio_iotype_pnetcdf, pio_iotype_netcdf
    use elm_instance,   only: lnd_id
    use shr_pio_mod,    only: shr_pio_getioformat
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout) :: file    ! PIO file descriptor
    character(len=*) ,  intent(in)    :: fname   ! File name to create

    ! BUG(wjs, 2014-10-20, bugz 1730) Workaround for
    ! http://bugs.cgd.ucar.edu/show_bug.cgi?id=1730
    logical, intent(in), optional     :: avoid_pnetcdf
    !
    ! !LOCAL VARIABLES:
    logical :: l_avoid_pnetcdf  ! local version of avoid_pnetcdf
    integer :: my_io_type
    integer :: ierr
    integer :: mode             ! IO mode according to netcdf format used
    !-----------------------------------------------------------------------

    l_avoid_pnetcdf = .false.
    if (present(avoid_pnetcdf)) then
       l_avoid_pnetcdf = avoid_pnetcdf
    end if

    my_io_type = io_type
    if (l_avoid_pnetcdf) then
       if (my_io_type == pio_iotype_pnetcdf) then
          my_io_type = pio_iotype_netcdf
          if (pio_iotask_rank(pio_subsystem)==0) then
             write(iulog,*) 'Workaround for bugz 1730: creating'
             write(iulog,*) trim(fname)
             write(iulog,*) 'with type netcdf instead of pnetcdf'
          end if
       end if
    end if

    mode = shr_pio_getioformat(lnd_id)
    ierr = pio_createfile(pio_subsystem, file, my_io_type, fname, ior(PIO_CLOBBER,mode))

    if(ierr/= PIO_NOERR) then
       call shr_sys_abort( ' ncd_pio_createfile ERROR: Failed to open file to write: '//trim(fname))
    else if(pio_iotask_rank(pio_subsystem)==0) then
       write(iulog,*) 'Opened file ', trim(fname),  ' to write', file%fh
    end if

  end subroutine ncd_pio_createfile

  !-----------------------------------------------------------------------
  subroutine check_var(ncid, varname, vardesc, readvar, print_err )
    !
    ! !DESCRIPTION:
    ! Check if variable is on netcdf file
    !
    ! !ARGUMENTS:
    class(file_desc_t) , intent(inout)  :: ncid      ! PIO file descriptor
    character(len=*)   , intent(in)     :: varname   ! Varible name to check
    type(Var_desc_t)   , intent(out)    :: vardesc   ! Output variable descriptor
    logical            , intent(out)    :: readvar   ! If variable exists or not
    logical, optional  , intent(in)     :: print_err ! If should print about error
    !
    ! !LOCAL VARIABLES:
    integer :: ret     ! return value
    logical :: log_err ! if should log error
    character(len=*),parameter :: subname='check_var' ! subroutine name
    !-----------------------------------------------------------------------


    if ( present(print_err) )then
       log_err = print_err
    else
       log_err = .true.
    end if
    readvar = .true.
    call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
    ret = PIO_inq_varid (ncid, varname, vardesc)
    if (ret /= PIO_noerr) then
       readvar = .false.
       if (masterproc .and. log_err) &
            write(iulog,*) subname//': variable ',trim(varname),' is not on dataset'
    end if
    call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)

  end subroutine check_var

  !-----------------------------------------------------------------------
  subroutine check_att(ncid, varid, attrib, att_found)
    !
    ! !DESCRIPTION:
    ! Check if attribute is on file
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    integer           ,intent(in)    :: varid     ! netcdf var id
    character(len=*)  ,intent(in)    :: attrib    ! netcdf attrib
    logical           ,intent(out)   :: att_found ! true if the attribute was found
    !
    ! !LOCAL VARIABLES:
    integer :: att_type  ! attribute type
    integer(PIO_OFFSET_KIND) :: att_len   ! attribute length
    integer :: status

    character(len=*), parameter :: subname = 'check_att'
    !-----------------------------------------------------------------------
    
    att_found = .true.
    call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
    status = PIO_inq_att(ncid, varid, trim(attrib), att_type, att_len)
    if (status /= PIO_noerr) then
       att_found = .false.
    end if
    call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)

  end subroutine check_att

  !-----------------------------------------------------------------------
  subroutine check_dim(ncid, dimname, value)
    !
    ! !DESCRIPTION:
    ! Validity check on dimension
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(in) :: ncid      ! PIO file handle
    character(len=*) , intent(in) :: dimname   ! Dimension name
    integer, intent(in)           :: value     ! Expected dimension size
    !
    ! !LOCAL VARIABLES:
    integer :: dimid, dimlen    ! temporaries
    integer :: status           ! error code      
    character(len=*),parameter :: subname='check_dim' ! subroutine name
    !-----------------------------------------------------------------------

    status = pio_inq_dimid (ncid, trim(dimname), dimid)
    status = pio_inq_dimlen (ncid, dimid, dimlen)
    if (dimlen /= value) then
       write(iulog,*) subname//' ERROR: mismatch of input dimension ',dimlen, &
            ' with expected value ',value,' for variable ',trim(dimname)
       call shr_sys_abort(errMsg(__FILE__, __LINE__))
    end if

  end subroutine check_dim

  !-----------------------------------------------------------------------
  subroutine ncd_enddef(ncid)
    !
    ! !DESCRIPTION:
    ! enddef netcdf file
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    !
    ! !LOCAL VARIABLES:
    integer :: status   ! error status
    !-----------------------------------------------------------------------

    status = PIO_enddef(ncid)

  end subroutine ncd_enddef

  !-----------------------------------------------------------------------
  subroutine ncd_inqdid(ncid,name,dimid,dimexist)
    !
    ! !DESCRIPTION:
    ! inquire on a dimension id
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid   ! netcdf file id
    character(len=*) , intent(in) :: name      ! dimension name
    integer          , intent(out):: dimid     ! dimension id
    logical,optional , intent(out):: dimexist  ! if this dimension exists or not
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    if ( present(dimexist) )then
       call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
    end if
    status = PIO_inq_dimid(ncid,name,dimid)
    if ( present(dimexist) )then
       if ( status == PIO_NOERR)then
          dimexist = .true.
       else
          dimexist = .false.
       end if
       call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
    end if

  end subroutine ncd_inqdid

  !-----------------------------------------------------------------------
  subroutine ncd_inqdlen(ncid,dimid,len,name)
    !
    ! !DESCRIPTION:
    ! enddef netcdf file
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout) :: ncid       ! netcdf file id
    integer           , intent(inout) :: dimid      ! dimension id
    integer           , intent(out)   :: len        ! dimension len
    character(len=*), optional, intent(in) :: name  ! dimension name
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    if ( present(name) )then
       call ncd_inqdid(ncid,name,dimid)
    end if
    len = -1
    status = PIO_inq_dimlen(ncid,dimid,len)

  end subroutine ncd_inqdlen

  !-----------------------------------------------------------------------
  subroutine ncd_inqdname(ncid,dimid,dname)
    !
    ! !DESCRIPTION:
    ! inquire dim name
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(in) :: ncid      ! netcdf file id
    integer           , intent(in) :: dimid     ! dimension id
    character(len=*)  , intent(out):: dname     ! dimension name
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    status = PIO_inq_dimname(ncid,dimid,dname)

  end subroutine ncd_inqdname

  !-----------------------------------------------------------------------
  subroutine ncd_inqfdims(ncid, isgrid2d, ni, nj, ns)
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout):: ncid
    logical           , intent(out)  :: isgrid2d
    integer           , intent(out)  :: ni
    integer           , intent(out)  :: nj
    integer           , intent(out)  :: ns
    !
    ! !LOCAL VARIABLES:
    integer  :: dimid                                ! netCDF id
    integer  :: ier                                  ! error status 
    character(len=32) :: subname = 'ncd_inqfdims' ! subroutine name
    !-----------------------------------------------------------------------

    if (single_column) then
       ni = 1
       nj = 1
       ns = 1
       isgrid2d = .true.
       RETURN
    end if

    if (scm_multcols) then
       ni = scm_nx
       nj = scm_ny
       ns = ni * nj
       isgrid2d = .true.
       RETURN
    end if

    ni = 0
    nj = 0

    call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
    ier = pio_inq_dimid (ncid, 'lon', dimid)
    if (ier == PIO_NOERR) ier = pio_inq_dimlen(ncid, dimid, ni)
    ier = pio_inq_dimid (ncid, 'lat', dimid)
    if (ier == PIO_NOERR) ier = pio_inq_dimlen(ncid, dimid, nj)

    ier = pio_inq_dimid (ncid, 'lsmlon', dimid)
    if (ier == PIO_NOERR) ier = pio_inq_dimlen(ncid, dimid, ni)
    ier = pio_inq_dimid (ncid, 'lsmlat', dimid)
    if (ier == PIO_NOERR) ier = pio_inq_dimlen(ncid, dimid, nj)

    ier = pio_inq_dimid (ncid, 'ni', dimid)
    if (ier == PIO_NOERR) ier = pio_inq_dimlen(ncid, dimid, ni)
    ier = pio_inq_dimid (ncid, 'nj', dimid)
    if (ier == PIO_NOERR) ier = pio_inq_dimlen(ncid, dimid, nj)

    ier = pio_inq_dimid (ncid, 'gridcell', dimid)
    if (ier == PIO_NOERR) then
       ier = pio_inq_dimlen(ncid, dimid, ni)
       if (ier == PIO_NOERR) nj = 1
    end if
   
    call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)

    if (ni == 0 .or. nj == 0) then
       write(iulog,*) trim(subname),' ERROR: ni,nj = ',ni,nj,' cannot be zero '
       call shr_sys_abort(errMsg(__FILE__, __LINE__))
    end if

    if (nj == 1) then
       isgrid2d = .false.
    else
       isgrid2d = .true.
    end if

    ns = ni*nj

  end subroutine ncd_inqfdims

  !-----------------------------------------------------------------------
  subroutine ncd_inqvid(ncid,name,varid,vardesc,readvar)
    !
    ! !DESCRIPTION:
    ! Inquire on a variable ID
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout) :: ncid      ! netcdf file id
    character(len=*)  , intent(in)    :: name      ! variable name
    integer           , intent(out)   :: varid     ! variable id
    type(Var_desc_t)  , intent(out)   :: vardesc   ! variable descriptor
    logical, optional , intent(out)   :: readvar   ! does variable exist
    !
    ! !LOCAL VARIABLES:
    integer :: ret               ! return code
    character(len=*),parameter :: subname='ncd_inqvid' ! subroutine name
    !-----------------------------------------------------------------------

    if (present(readvar)) then
       readvar = .false.
       call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
       ret = PIO_inq_varid(ncid,name,vardesc)
       if (ret /= PIO_noerr) then
          if (masterproc) write(iulog,*) subname//': variable ',trim(name),' is not on dataset'
          readvar = .false.
       else
          readvar = .true.
       end if
       call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
    else
       ret = PIO_inq_varid(ncid,name,vardesc)
    endif
    varid = vardesc%varid

  end subroutine ncd_inqvid

  !-----------------------------------------------------------------------
  subroutine ncd_inqvdims(ncid,ndims,vardesc)
    !
    ! !DESCRIPTION:
    ! inquire variable dimensions
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(in)   :: ncid      ! netcdf file id
    integer           , intent(out)  :: ndims     ! variable ndims
    type(Var_desc_t)  , intent(inout):: vardesc   ! variable descriptor
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    ndims = -1
    status = PIO_inq_varndims(ncid,vardesc,ndims)

  end subroutine ncd_inqvdims

  !-----------------------------------------------------------------------
  subroutine ncd_inqvname(ncid,varid,vname,vardesc)
    !
    ! !DESCRIPTION:
    ! inquire variable name
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(in)   :: ncid      ! netcdf file id
    integer           , intent(in)   :: varid     ! variable id
    character(len=*)  , intent(out)  :: vname     ! variable vname
    type(Var_desc_t)  , intent(inout):: vardesc   ! variable descriptor
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    vname = ''
    status = PIO_inq_varname(ncid,vardesc,vname)

  end subroutine ncd_inqvname

  !-----------------------------------------------------------------------
  subroutine ncd_inqvdids(ncid,dids,vardesc)
    !
    ! !DESCRIPTION:
    ! inquire variable dimension ids
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(in)  :: ncid      ! netcdf file id
    integer          ,intent(out)  :: dids(:)   ! variable dids
    type(Var_desc_t) ,intent(inout):: vardesc   ! variable descriptor
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    dids = -1
    status = PIO_inq_vardimid(ncid,vardesc,dids)

  end subroutine ncd_inqvdids

  !-----------------------------------------------------------------------
  subroutine ncd_inqvdlen_byDesc(ncid,vardesc,dimnum,dlen,err_code)
    !
    ! !DESCRIPTION:
    ! inquire size of one of a variable's dimensions, given a vardesc
    !
    ! If the variable has n dimensions, then dimnum should be between 1 and n; this routine
    ! returns the size of the dimnum'th dimension.
    !
    ! If there is an error condition, dlen will be -1, and err_code will hold the error
    ! code; possible error codes are:
    ! 0: no error
    ! 1: dimnum out of range
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    type(Var_desc_t)  ,intent(inout) :: vardesc   ! variable descriptor
    integer           ,intent(in)    :: dimnum    ! dimension number to query
    integer           ,intent(out)   :: dlen      ! length of the dimension
    integer           ,intent(out)   :: err_code  ! error code (0 means no error)
    !
    ! !LOCAL VARIABLES:
    integer              :: ndims      ! number of dimensions
    integer, allocatable :: dimids(:)  ! dimension IDs

    integer, parameter   :: dlen_invalid = -1
    integer, parameter   :: error_none = 0
    integer, parameter   :: error_dimnum_out_of_range = 1
    !-----------------------------------------------------------------------

    err_code = error_none

    call ncd_inqvdims(ncid, ndims, vardesc)

    if (dimnum > 0 .and. dimnum <= ndims) then
       allocate(dimids(ndims))
       call ncd_inqvdids(ncid, dimids, vardesc)
       call ncd_inqdlen(ncid, dimids(dimnum), dlen)
       deallocate(dimids)
    else
       dlen = dlen_invalid
       err_code = error_dimnum_out_of_range
    end if

  end subroutine ncd_inqvdlen_byDesc


  !-------------------------------------------------------------------------
  subroutine ncd_inqvdlen_byName(ncid,varname,dimnum,dlen,err_code)
    !
    ! !DESCRIPTION:
    ! inquire size of one of a variable's dimensions, given a variable name
    !
    ! If the variable has n dimensions, then dimnum should be between 1 and n; this routine
    ! returns the size of the dimnum'th dimension.
    !
    ! If there is an error condition, dlen will be -1, and err_code will hold the error
    ! code; possible error codes are:
    ! 0: no error
    ! 1: dimnum out of range
    ! 11: variable not found
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    character(len=*)  ,intent(in)    :: varname   ! variable name
    integer           ,intent(in)    :: dimnum    ! dimension number to query
    integer           ,intent(out)   :: dlen      ! length of the dimension
    integer           ,intent(out)   :: err_code  ! error code (0 means no error)
    !
    ! !LOCAL VARIABLES:
    type(Var_desc_t) :: vardesc        ! variable descriptor
    logical          :: readvar        ! whether the variable was found
    integer, parameter :: dlen_invalid = -1
    integer, parameter :: error_variable_not_found = 11
    !-----------------------------------------------------------------------

    call check_var(ncid, varname, vardesc, readvar)
    if (readvar) then
       call ncd_inqvdlen_byDesc(ncid, vardesc, dimnum, dlen, err_code)
    else
       dlen = dlen_invalid
       err_code = error_variable_not_found
    end if

 end subroutine ncd_inqvdlen_byName


 !-----------------------------------------------------------------------
 subroutine ncd_inqvdname_byDesc(ncid,vardesc,dimnum,dname,err_code)
    !
    ! !DESCRIPTION:
    ! Inquire name of one of a variable's dimensions, given a vardesc
    !
    ! If the variable has n dimensions, then dimnum should be between 1 and n; this
    ! routine returns the name of the dimnum'th dimension.
    !
    ! If there is an error condition, dname will be ' ', and err_code will hold the error
    ! code; possible error codes are:
    ! 0: no error
    ! 1: dimnum out of range
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    type(Var_desc_t)  ,intent(inout) :: vardesc   ! variable descriptor
    integer           ,intent(in)    :: dimnum    ! dimension number to query
    character(len=*)  ,intent(out)   :: dname     ! name of the dimension
    integer           ,intent(out)   :: err_code  ! error code (0 means no error)
    !
    ! !LOCAL VARIABLES:
    integer              :: ndims      ! number of dimensions
    integer, allocatable :: dimids(:)  ! dimension IDs

    character(len=*), parameter :: dname_invalid = ' '
    integer, parameter :: error_none = 0
    integer, parameter :: error_dimnum_out_of_range = 1

    character(len=*), parameter :: subname = 'ncd_inqvdname_byDesc'
    !-----------------------------------------------------------------------
    err_code = error_none

    call ncd_inqvdims(ncid, ndims, vardesc)

    if (dimnum > 0 .and. dimnum <= ndims) then
       allocate(dimids(ndims))
       call ncd_inqvdids(ncid, dimids, vardesc)
       call ncd_inqdname(ncid, dimids(dimnum), dname)
       deallocate(dimids)
    else
       dname = dname_invalid
       err_code = error_dimnum_out_of_range
    end if

 end subroutine ncd_inqvdname_byDesc

 !-----------------------------------------------------------------------

 subroutine ncd_inqvdname_byName(ncid,varname,dimnum,dname,err_code)
    !
    ! !DESCRIPTION:
    ! Inquire name of one of a variable's dimensions, given a variable name
    !
    ! If the variable has n dimensions, then dimnum should be between 1 and n; this
    ! routine returns the name of the dimnum'th dimension.
    !
    ! If there is an error condition, dname will be ' ', and err_code will hold the error
    ! code; possible error codes are:
    ! 0: no error
    ! 1: dimnum out of range
    ! 11: variable not found
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    character(len=*)  ,intent(in)    :: varname   ! variable name
    integer           ,intent(in)    :: dimnum    ! dimension number to query
    character(len=*)  ,intent(out)   :: dname     ! name of the dimension
    integer           ,intent(out)   :: err_code  ! error code (0 means no error)
    !
    ! !LOCAL VARIABLES:
    type(Var_desc_t) :: vardesc        ! variable descriptor
    logical          :: readvar        ! whether the variable was found
    character(len=*), parameter :: dname_invalid = ' '
    integer, parameter :: error_variable_not_found = 11

    character(len=*), parameter :: subname = 'ncd_inqvdname_byName'
    !-----------------------------------------------------------------------

    call check_var(ncid, varname, vardesc, readvar)
    if (readvar) then
       call ncd_inqvdname_byDesc(ncid, vardesc, dimnum, dname, err_code)
    else
       dname = dname_invalid
       err_code = error_variable_not_found
    end if

 end subroutine ncd_inqvdname_byName

  !-----------------------------------------------------------------------

  subroutine ncd_putatt_int(ncid,varid,attrib,value,xtype)
    !
    ! !DESCRIPTION:
    ! put integer attributes
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    integer           ,intent(in)    :: varid     ! netcdf var id
    character(len=*)  ,intent(in)    :: attrib    ! netcdf attrib
    integer           ,intent(in)    :: value     ! netcdf attrib value
    integer,optional  ,intent(in)    :: xtype     ! netcdf data type
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    status = PIO_put_att(ncid,varid,trim(attrib),value)

  end subroutine ncd_putatt_int

  !-----------------------------------------------------------------------
  subroutine ncd_putatt_char(ncid,varid,attrib,value,xtype)
    !
    ! !DESCRIPTION:
    ! put character attributes
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    integer           ,intent(in)    :: varid     ! netcdf var id
    character(len=*)  ,intent(in)    :: attrib    ! netcdf attrib
    character(len=*)  ,intent(in)    :: value     ! netcdf attrib value
    integer,optional  ,intent(in)    :: xtype     ! netcdf data type
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    status = PIO_put_att(ncid,varid,trim(attrib),value)

  end subroutine ncd_putatt_char

  !-----------------------------------------------------------------------
  subroutine ncd_putatt_real(ncid,varid,attrib,value,xtype)
    !
    ! !DESCRIPTION:
    ! put real attributes
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    integer           ,intent(in)    :: varid     ! netcdf var id
    character(len=*)  ,intent(in)    :: attrib    ! netcdf attrib
    real(r8)          ,intent(in)    :: value     ! netcdf attrib value
    integer           ,intent(in)    :: xtype     ! netcdf data type
    !
    ! !LOCAL VARIABLES:
    integer :: status
    real*4  :: value4
    !-----------------------------------------------------------------------

    value4 = value

    if (xtype == pio_double) then
       status = PIO_put_att(ncid,varid,trim(attrib),value)
    else
       status = PIO_put_att(ncid,varid,trim(attrib),value4)
    endif

  end subroutine ncd_putatt_real

  !-----------------------------------------------------------------------
  subroutine ncd_getatt_char(ncid,varid,attrib,value)
    !
    ! !DESCRIPTION:
    ! get a character attribute
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid      ! netcdf file id
    integer           ,intent(in)    :: varid     ! netcdf var id
    character(len=*)  ,intent(in)    :: attrib    ! netcdf attrib
    character(len=*)  ,intent(out)   :: value     ! netcdf attrib value
    !
    ! !LOCAL VARIABLES:
    integer :: status
    
    character(len=*), parameter :: subname = 'ncd_getatt_char'
    !-----------------------------------------------------------------------
    
    status = PIO_get_att(ncid,varid,trim(attrib),value)

  end subroutine ncd_getatt_char


  !-----------------------------------------------------------------------
  subroutine ncd_defdim(ncid,attrib,value,dimid)
    !
    ! !DESCRIPTION:
    ! define dimension
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(in) :: ncid      ! netcdf file id
    character(len=*)  , intent(in) :: attrib    ! netcdf attrib
    integer           , intent(in) :: value     ! netcdf attrib value
    integer           , intent(out):: dimid     ! netcdf dimension id
    !
    ! !LOCAL VARIABLES:
    integer :: status
    !-----------------------------------------------------------------------

    status = pio_def_dim(ncid,attrib,value,dimid)

  end subroutine ncd_defdim

  !-----------------------------------------------------------------------
  subroutine ncd_defvar_bynf(ncid, varname, xtype, ndims, dimid, varid, &
       long_name, standard_name, units, cell_method, missing_value, fill_value, &
       imissing_value, ifill_value, comment, flag_meanings, &
       flag_values, nvalid_range )
    !
    ! !DESCRIPTION:
    !  Define a netcdf variable
    !
    ! !ARGUMENTS:
    class(file_desc_t) , intent(inout) :: ncid                  ! netcdf file id
    character(len=*)   , intent(in)  :: varname                 ! variable name
    integer            , intent(in)  :: xtype                   ! external type
    integer            , intent(in)  :: ndims                   ! number of dims
    integer            , intent(inout) :: varid                 ! returned var id
    integer            , intent(in), optional :: dimid(:)       ! dimids
    character(len=*)   , intent(in), optional :: long_name      ! attribute
    character(len=*)   , intent(in), optional :: standard_name  ! attribute
    character(len=*)   , intent(in), optional :: units          ! attribute
    character(len=*)   , intent(in), optional :: cell_method    ! attribute
    character(len=*)   , intent(in), optional :: comment        ! attribute
    character(len=*)   , intent(in), optional :: flag_meanings(:) ! attribute
    real(r8)           , intent(in), optional :: missing_value  ! attribute for real
    real(r8)           , intent(in), optional :: fill_value     ! attribute for real
    integer            , intent(in), optional :: imissing_value ! attribute for int
    integer            , intent(in), optional :: ifill_value    ! attribute for int
    integer            , intent(in), optional :: flag_values(:)  ! attribute for int
    integer            , intent(in), optional :: nvalid_range(2)  ! attribute for int
    !
    ! !LOCAL VARIABLES:
    integer :: n                   ! indices
    integer :: ldimid(4)           ! local dimid 
    integer :: dimid0(1)           ! local dimid
    integer :: status              ! error status 
    integer :: lxtype              ! local external type (in case logical variable)
    type(var_desc_t)   :: vardesc  ! local vardesc
    character(len=128) :: dimname  ! temporary
    character(len=256) :: str      ! temporary
    character(len=*),parameter :: subname='ncd_defvar_bynf' ! subroutine name
    !-----------------------------------------------------------------------

    varid = -1

    dimid0 = 0
    ldimid = 0
    if (present(dimid)) then
       ldimid(1:ndims) = dimid(1:ndims)
    else   ! ndims must be zero if dimid not present
       if (ndims /= 0) then
          write(iulog,*) subname//' ERROR: dimid not supplied and ndims ne 0 ',trim(varname),ndims
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       endif
    endif

    if ( xtype == ncd_log )then
       lxtype = ncd_int
    else
       lxtype = xtype
    end if
    if (masterproc .and. debug > 1) then
       write(iulog,*) 'Error in defining variable = ', trim(varname)
       write(iulog,*) subname//' ',trim(varname),lxtype,ndims,ldimid(1:ndims)
    endif

    if (ndims >  0) then 
       status = pio_inq_dimname(ncid,ldimid(ndims),dimname)
    end if

    ! Define variable
    if (present(dimid)) then
       status = PIO_def_var(ncid,trim(varname),lxtype,dimid(1:ndims),vardesc)
    else
       status = PIO_def_var(ncid,trim(varname),lxtype,dimid0        ,vardesc)
    endif
    varid = vardesc%varid

    !
    ! Add attributes
    !
    if (present(long_name)) then
       call ncd_putatt(ncid, varid, 'long_name', trim(long_name))
    end if
    if (present(standard_name)) then
       if (standard_name /= ' ') then
          call ncd_putatt(ncid, varid, 'standard_name', trim(standard_name))
       end if
    end if
    if (present(flag_values)) then
       status = PIO_put_att(ncid,varid,'flag_values',flag_values)
       if ( .not. present(flag_meanings)) then
          write(iulog,*) 'Error in defining variable = ', trim(varname)
          call shr_sys_abort(" ERROR:: flag_values set -- but not flag_meanings"//errMsg(__FILE__, __LINE__))
       end if
    end if
    if (present(flag_meanings)) then
       if ( .not. present(flag_values)) then
          write(iulog,*) 'Error in defining variable = ', trim(varname)
          call shr_sys_abort(" ERROR:: flag_meanings set -- but not flag_values"//errMsg(__FILE__, __LINE__) )
       end if
       if ( size(flag_values) /= size(flag_meanings) ) then
          write(iulog,*) 'Error in defining variable = ', trim(varname)
          call shr_sys_abort(" ERROR:: flag_meanings and flag_values dimension different"//errMsg(__FILE__, __LINE__))
       end if
       str = flag_meanings(1)
       do n = 1, size(flag_meanings)
          if ( index(flag_meanings(n), ' ') /= 0 )then
             write(iulog,*) 'Error in defining variable = ', trim(varname)
             call shr_sys_abort(" ERROR:: flag_meanings has an invalid space in it"//errMsg(__FILE__, __LINE__) )
          end if
          if ( n > 1 ) str = trim(str)//" "//flag_meanings(n)
       end do
       status = PIO_put_att(ncid,varid,'flag_meanings', trim(str) )
    end if
    if (present(comment)) then
       call ncd_putatt(ncid, varid, 'comment', trim(comment))
    end if
    if (present(units)) then
       call ncd_putatt(ncid, varid, 'units', trim(units))
    end if
    if (present(cell_method)) then
       str = 'time: ' // trim(cell_method)
       call ncd_putatt(ncid, varid, 'cell_methods', trim(str))
    end if
    if (present(fill_value)) then
       call ncd_putatt(ncid, varid, '_FillValue', fill_value, lxtype)
    end if
    if (present(missing_value)) then
       call ncd_putatt(ncid, varid, 'missing_value', missing_value, lxtype)
    end if
    if (present(ifill_value)) then
       call ncd_putatt(ncid, varid, '_FillValue', ifill_value, lxtype)
    end if
    if (present(imissing_value)) then
       call ncd_putatt(ncid, varid, 'missing_value', imissing_value, lxtype)
    end if
    if (present(nvalid_range)) then
       status = PIO_put_att(ncid,varid,'valid_range', nvalid_range )
    end if
    if ( xtype == ncd_log )then
       status = PIO_put_att(ncid,varid,'flag_values',     (/0, 1/) )
       status = PIO_put_att(ncid,varid,'flag_meanings',  "FALSE TRUE" )
       status = PIO_put_att(ncid,varid,'valid_range',    (/0, 1/) )
    end if

  end subroutine ncd_defvar_bynf

  !-----------------------------------------------------------------------
  subroutine ncd_defvar_bygrid(ncid, varname, xtype, &
       dim1name, dim2name, dim3name, dim4name, dim5name, &
       long_name, standard_name, units, cell_method, missing_value, fill_value, &
       imissing_value, ifill_value, switchdim, comment, &
       flag_meanings, flag_values, nvalid_range )
    !
    ! !DESCRIPTION:
    !  Define a netcdf variable
    !
    ! !ARGUMENTS:
    class(file_desc_t) , intent(inout) :: ncid                 ! netcdf file id
    character(len=*)   , intent(in)  :: varname                 ! variable name
    integer            , intent(in)  :: xtype                   ! external type
    character(len=*)   , intent(in), optional :: dim1name       ! dimension name
    character(len=*)   , intent(in), optional :: dim2name       ! dimension name
    character(len=*)   , intent(in), optional :: dim3name       ! dimension name
    character(len=*)   , intent(in), optional :: dim4name       ! dimension name
    character(len=*)   , intent(in), optional :: dim5name       ! dimension name
    character(len=*)   , intent(in), optional :: long_name      ! attribute
    character(len=*)   , intent(in), optional :: standard_name  ! attribute
    character(len=*)   , intent(in), optional :: units          ! attribute
    character(len=*)   , intent(in), optional :: cell_method    ! attribute
    character(len=*)   , intent(in), optional :: comment        ! attribute
    character(len=*)   , intent(in), optional :: flag_meanings(:) ! attribute
    real(r8)           , intent(in), optional :: missing_value  ! attribute for real
    real(r8)           , intent(in), optional :: fill_value     ! attribute for real
    integer            , intent(in), optional :: imissing_value ! attribute for int
    integer            , intent(in), optional :: ifill_value    ! attribute for int
    logical            , intent(in), optional :: switchdim      ! true=> permute dim1 and dim2 for output
    integer            , intent(in), optional :: flag_values(:)  ! attribute for int
    integer            , intent(in), optional :: nvalid_range(2)  ! attribute for int
    !
    ! !LOCAL VARIABLES:
    integer :: n              ! indices
    integer :: ndims          ! dimension counter
    integer :: dimid(5)       ! dimension ids
    integer :: varid          ! variable id
    integer :: itmp           ! temporary
    character(len=256) :: str ! temporary
    character(len=*),parameter :: subname='ncd_defvar_bygrid' ! subroutine name
    !-----------------------------------------------------------------------

    dimid(:) = 0

    ! Determine dimension ids for variable

    if (present(dim1name)) call ncd_inqdid(ncid, dim1name, dimid(1))
    if (present(dim2name)) call ncd_inqdid(ncid, dim2name, dimid(2))
    if (present(dim3name)) call ncd_inqdid(ncid, dim3name, dimid(3))
    if (present(dim4name)) call ncd_inqdid(ncid, dim4name, dimid(4))
    if (present(dim5name)) call ncd_inqdid(ncid, dim5name, dimid(5))

    ! Permute dim1 and dim2 if necessary

    if (present(switchdim)) then
       itmp = dimid(2)
       dimid(2) = dimid(1)
       dimid(1) = itmp
    end if

    ! Define variable

    ndims = 0
    if (present(dim1name)) then
       do n = 1, size(dimid)
          if (dimid(n) /= 0) ndims = ndims + 1
       end do
    end if

    call ncd_defvar_bynf(ncid,varname,xtype,ndims,dimid,varid, &
         long_name=long_name, standard_name=standard_name,units=units, cell_method=cell_method, &
         missing_value=missing_value, fill_value=fill_value, &
         imissing_value=imissing_value, ifill_value=ifill_value, &
         comment=comment, flag_meanings=flag_meanings, &
         flag_values=flag_values, nvalid_range=nvalid_range )

  end subroutine ncd_defvar_bygrid

  !------------------------------------------------------------------------
  subroutine ncd_io_char_var0_start_glob(vardesc, data, flag, ncid, start )
    !
    ! !DESCRIPTION:
    ! netcdf I/O of global character array with start indices input
    !
    ! !ARGUMENTS:
    class(file_desc_t),intent(inout) :: ncid             ! netcdf file id
    character(len=*) , intent(in)    :: flag             ! 'read' or 'write'
    type(var_desc_t) , intent(in)    :: vardesc          ! local vardesc pointer
    character(len=*) , intent(inout) :: data             ! raw data for this index
    integer          , intent(in)    :: start(:)         ! output bounds
    !
    ! !LOCAL VARIABLES:
    integer :: status               ! error code
    character(len=*),parameter :: subname='ncd_io_char_var0_start_glob'
    !-----------------------------------------------------------------------

    if (flag == 'read') then

       status = pio_get_var(ncid, vardesc, start, data )

    elseif (flag == 'write') then

       status = pio_put_var(ncid, vardesc, start, data )

    endif

  end subroutine ncd_io_char_var0_start_glob

  !------------------------------------------------------------------------
  !DIMS 0,1
  subroutine ncd_io_{DIMS}d_log_glob(varname, data, flag, ncid, readvar, nt, posNOTonfile)
    !
    ! !DESCRIPTION:
    ! netcdf I/O of global integer variable
    !
    ! !ARGUMENTS:
    class(file_desc_t) , intent(inout) :: ncid      ! netcdf file id
    character(len=*)   , intent(in)    :: flag      ! 'read' or 'write'
    character(len=*)   , intent(in)    :: varname   ! variable name
    logical            , intent(inout) :: data{DIMSTR} ! raw data
    logical, optional  , intent(out)   :: readvar   ! was var read?
    integer, optional  , intent(in)    :: nt        ! time sample index
    logical            , optional, intent(in) :: posNOTonfile ! position is NOT on this file
    !
    ! !LOCAL VARIABLES:
    integer           :: varid              ! netCDF variable id
    integer           :: start(2), count(2) ! output bounds
    integer           :: status             ! error code
    logical           :: varpresent         ! if true, variable is on tape
    integer           :: idata
    integer, pointer  :: idata1d(:)         ! Temporary integer data to send to file
    character(len=32) :: vname              ! variable error checking
    type(var_desc_t)  :: vardesc            ! local vardesc pointer
    character(len=*),parameter :: subname='ncd_io_{DIMS}d_log_glob'
    !-----------------------------------------------------------------------

    start(:) = 0
    count(:) = 0

    if (flag == 'read') then

       call ncd_inqvid(ncid, varname, varid, vardesc, readvar=varpresent)
       if (varpresent) then
          if (single_column .and. present(posNOTonfile) ) then
             if ( .not. posNOTonfile )then
                call shr_sys_abort(' ERROR: scalar var is NOT compatable with posNOTonfile = .false.'//&
                     errMsg(__FILE__, __LINE__))
             end if
          endif
#if ({DIMS}==0) 
          status = pio_get_var(ncid, varid, idata)
          if ( idata == 0 )then
             data = .false.
          else if ( idata == 1 )then
             data = .true.
          else
             call shr_sys_abort(' ERROR: bad integer value for logical data'//errMsg(__FILE__, __LINE__))
          end if
#else
          allocate(idata1d(size(data))) 
          data = (idata == 1)
          if ( any(idata1d /= 0 .and. idata1d /= 1) )then
             call shr_sys_abort(' ERROR: read in bad integer value(s) for logical data'//errMsg(__FILE__, __LINE__))
          end if
          deallocate(idata1d)
#endif
       endif
       if (present(readvar)) readvar = varpresent

    elseif (flag == 'write') then

#if ({DIMS}==0) 
       start(1) = 1 ; count(1) = 1
       if (present(nt)) start(1) = nt
       call ncd_inqvid  (ncid, varname, varid, vardesc)
       allocate(idata1d(1))
       if ( data )then
          idata1d(1) = 1
       else
          idata1d(1) = 0
       end if
       status = pio_put_var(ncid, varid, start, count, idata1d)
       deallocate(idata1d)
#else
       start(1) = 1 ; count(1) = size(data)
       start(2) = 1 ; count(2) = 1
       if (present(nt)) start(2) = nt
       allocate(idata1d(size(data))) 
       where( data )
          idata1d = 1
       elsewhere
          idata1d = 0
       end where
       call ncd_inqvid  (ncid, varname, varid, vardesc)
       status = pio_put_var(ncid, varid, start, count, idata1d)
       deallocate( idata1d )
#endif

    endif   ! flag

  end subroutine ncd_io_{DIMS}d_log_glob

  !------------------------------------------------------------------------
  !DIMS 0,1,2,3,4
  !TYPE int,double
  subroutine ncd_io_{DIMS}d_{TYPE}_glob(varname, data, flag, ncid, readvar, nt, posNOTonfile) 
    !
    ! !DESCRIPTION:
    ! netcdf I/O of global variable
    !
    ! !ARGUMENTS:
    class(file_desc_t),         intent(inout) :: ncid         ! netcdf file id
    character(len=*),           intent(in)    :: flag         ! 'read' or 'write'
    character(len=*),           intent(in)    :: varname      ! variable name
    {VTYPE}         ,           intent(inout) :: data{DIMSTR} ! raw data
    logical         , optional, intent(out)   :: readvar      ! was var read?
    integer         , optional, intent(in)    :: nt           ! time sample index
    logical         , optional, intent(in)    :: posNOTonfile ! position is NOT on this file
    !
    ! !LOCAL VARIABLES:
    integer           :: m
    integer           :: varid              ! netCDF variable id
    integer           :: start(5), count(5) ! output bounds
    integer           :: status             ! error code
    logical           :: varpresent         ! if true, variable is on tape
    logical           :: found              ! if true, found lat/lon dims on file
    logical           :: scm_restart_file   ! is this a restart file in SCM mode?
    character(len=32) :: vname              ! variable error checking
    type(var_desc_t)  :: vardesc            ! local vardesc pointer
    {VTYPE} :: temp(1)
    character(len=*),parameter :: subname='ncd_io_{DIMS}d_{TYPE}_glob'

    {VTYPE}, allocatable :: iop_data(:,:)
    !-----------------------------------------------------------------------

    start(:) = 0
    count(:) = 0

    if (flag == 'read') then

       ! If running SCM mode then check to see if this is a restart file
       !   or not.  If it is a restart file then the standard PIO routines
       !   will be activated, if not then the specialized routines specific
       !   for SCM and SCM mult cols will be used.
       if (scm_multcols .or. single_column) then
         call is_restart_file_scmmode(ncid, scm_restart_file)
       else
         scm_restart_file = .false.
       endif

       call ncd_inqvid(ncid, varname, varid, vardesc, readvar=varpresent)

#if ({DIMS}==0) 
       if (varpresent) then
          status = pio_get_var(ncid, vardesc, data)
          if (single_column .and. present(posNOTonfile) ) then
             if ( .not. posNOTonfile )then
                call shr_sys_abort(' ERROR: scalar var is NOT compatable with posNOTonfile = .false.'//&
                     errMsg(__FILE__, __LINE__))
             end if
          endif
       end if
#else
       if (varpresent) then
          if ((scm_multcols .or. single_column) .and. .not. scm_restart_file) then
             call scam_field_offsets(ncid,'undefined', vardesc,&
                  start, count, found=found, posNOTonfile=posNOTonfile)
             if ( found )then
                if (single_column) status = pio_get_var(ncid, varid, start, count, data)
                if (scm_multcols) then
                  allocate(iop_data(1:count(1),1:count(2)))
                  status = pio_get_var(ncid, varid, start(1:2), count(1:2), iop_data)
#if ({DIMS}==1)
                  data(:) = iop_data(1,1) ! Use data at IOP lat/lon for whole globe
#endif
                  deallocate(iop_data)
                endif
             else
                status = pio_get_var(ncid, varid, data)
             end if
          else
             status = pio_get_var(ncid, varid, data)
          endif
       endif
#endif
       if (present(readvar)) readvar = varpresent

    elseif (flag == 'write') then

       call ncd_inqvid  (ncid, varname, varid, vardesc)

       if (present(nt)) then
#if ({DIMS}==0) 
          start(1) = nt  ; count(1) = 1
          temp(1) = data
          status = pio_put_var(ncid, varid, start, count, temp)
#elif ({DIMS}==1) 
          start(1) = 1 ; count(1) = size(data) 
          start(2) = nt ; count(2) = 1
          status = pio_put_var(ncid, varid, start, count, data)
#elif ({DIMS}==2) 
          start(1) = 1  ; count(1) = size(data, dim=1)
          start(2) = 1  ; count(2) = size(data, dim=2)
          start(3) = nt  ; count(3) = 1
          status = pio_put_var(ncid, varid, start, count, data)
#elif ({DIMS}==3) 
          start(1) = 1  ; count(1) = size(data,dim=1)
          start(2) = 1  ; count(2) = size(data,dim=2)
          start(3) = 1  ; count(3) = size(data,dim=3)
          start(4) = nt ; count(4) = 1
          status = pio_put_var(ncid, varid, start, count, data)
#elif ({DIMS}==4) 
          start(1) = 1  ; count(1) = size(data,dim=1)
          start(2) = 1  ; count(2) = size(data,dim=2)
          start(3) = 1  ; count(3) = size(data,dim=3)
          start(4) = 1  ; count(4) = size(data,dim=4)
          start(5) = nt ; count(5) = 1
          status = pio_put_var(ncid, varid, start, count, data)
#endif
       else
          status = pio_put_var(ncid, varid, data)
       end if

    endif

end subroutine ncd_io_{DIMS}d_{TYPE}_glob

  !------------------------------------------------------------------------
  !DIMS 0,1,2
  !TYPE text
  subroutine ncd_io_{DIMS}d_{TYPE}_glob(varname, data, flag, ncid, readvar, nt, posNOTonfile)
    !
    ! !DESCRIPTION:
    ! netcdf I/O of global variable
    !
    ! !ARGUMENTS:
    class(file_desc_t),         intent(inout) :: ncid         ! netcdf file id
    character(len=*),           intent(in)    :: flag         ! 'read' or 'write'
    character(len=*),           intent(in)    :: varname      ! variable name
    {VTYPE}         ,           intent(inout) :: data{DIMSTR} ! raw data
    logical         , optional, intent(out)   :: readvar      ! was var read?
    integer         , optional, intent(in)    :: nt           ! time sample index
    logical         , optional, intent(in)    :: posNOTonfile ! position is NOT on this file
    !
    ! !LOCAL VARIABLES:
    integer           :: m
    integer           :: varid              ! netCDF variable id
    integer           :: start(4), count(4) ! output bounds
    integer           :: status             ! error code
    logical           :: varpresent         ! if true, variable is on tape
    character(len=1)  :: tmpString(max_string_len)     ! temp for manipulating output string
    type(var_desc_t)  :: vardesc            ! local vardesc pointer
    character(len=*),parameter :: subname='ncd_io_{DIMS}d_{TYPE}_glob'
    !-----------------------------------------------------------------------

    start(:) = 0
    count(:) = 0

    if (flag == 'read') then

       call ncd_inqvid(ncid, varname, varid, vardesc, readvar=varpresent)

       if (varpresent) then
          data   = ' '
          status = pio_get_var(ncid, varid, data)
       end if
       if (present(readvar)) readvar = varpresent

    elseif (flag == 'write') then

       call ncd_inqvid  (ncid, varname, varid, vardesc)

#if ({DIMS}==0) 
       if (present(nt)) then
          do m = 1,len(data)
             tmpString(m:m) = data(m:m)
          end do
          start(1) = 1 ; count(1) = len(data)
          start(2) = nt; count(2) = 1
          if ( count(1) > size(tmpString) )then
             write(iulog,*) subname//' ERROR: input string size is too large:'
          end if
          status = pio_put_var(ncid, varid, start, count, tmpString(1:count(1)))
       else
          status = pio_put_var(ncid, varid, data )
       end if
#elif ({DIMS}==1) 
       if (present(nt)) then
          start(1) = 1 ; count(1) = len(data)
          start(2) = 1 ; count(2) = size(data)
          start(3) = nt; count(3) = 1
          status = pio_put_var(ncid, varid, start, count, data)
       else
          status = pio_put_var(ncid, varid, data)
       end if
#elif ({DIMS}==2) 
       if (present(nt)) then
          start(1) = 1  ; count(1) = len(data)
          start(2) = 1  ; count(2) = size(data,dim=1) 
          start(3) = 1  ; count(3) = size(data,dim=2)
          start(4) = nt ; count(4) = 1
          status = pio_put_var(ncid, varid, start, count, data)
       else
          status = pio_put_var(ncid, varid, data)
       end if
#endif

    endif

  end subroutine ncd_io_{DIMS}d_{TYPE}_glob 

  !-----------------------------------------------------------------------

  !TYPE int,double,logical
  subroutine ncd_io_1d_{TYPE}(varname, data, dim1name, flag, ncid, nt, readvar, cnvrtnan2fill)
    !
    ! !DESCRIPTION:
    ! netcdf I/O for 1d 
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout)        :: ncid          ! netcdf file id
    character(len=*) , intent(in)            :: flag          ! 'read' or 'write'
    character(len=*) , intent(in)            :: varname       ! variable name
    {VTYPE}          , pointer               :: data(:)       ! local decomposition data
    character(len=*) , intent(in)            :: dim1name      ! dimension name
    integer          , optional, intent(in)  :: nt            ! time sample index
    logical          , optional, intent(out) :: readvar       ! true => variable is on initial dataset (read only)
    logical          , optional, intent(in)  :: cnvrtnan2fill ! true => convert any NaN's to _FillValue (spval)
    !
    ! Local Variables

    {VTYPE}, allocatable :: iop_data2d(:,:)
    {VTYPE}, allocatable :: iop_data1d(:)

    character(len=8)                 :: elmlevel   ! elmlevel
    character(len=32)                :: dimname    ! temporary
    integer                          :: n          ! index      
    integer                          :: n_time_ind !time slice index of the "start" array
    integer                          :: iodnum     ! iodesc num in list
    integer                          :: varid      ! varid
    integer                          :: ndims      ! ndims for var
    integer                          :: ndims_iod  ! ndims iodesc for var
    integer                          :: dims(4)    ! dim sizes       
    integer                          :: dids(4)    ! dim ids
    integer                          :: start(3)   ! netcdf start index
    integer                          :: count(3)   ! netcdf count index
    integer                          :: start_iop2d(2)   ! netcdf start index
    integer                          :: count_iop2d(2)   ! netcdf count index
    integer                          :: start_iop1d(1)   ! netcdf start index
    integer                          :: count_iop1d(1)   ! netcdf count index
    integer                          :: status     ! error code  
    logical                          :: varpresent ! if true, variable is on tape
    logical                          :: scm_restart_file ! restart file in SCM mode?
    integer                , pointer :: idata(:)   ! Temporary integer data to send to file
    integer                , pointer :: compDOF(:)
    type(iodesc_plus_type) , pointer :: iodesc_plus
    type(var_desc_t)                 :: vardesc
    character(len=*),parameter       :: subname='ncd_io_1d_{TYPE}' ! subroutine name
    !-----------------------------------------------------------------------

    start(:) = 0
    count(:) = 0

    elmlevel = dim1name

    if (scm_multcols .or. single_column) then
      call is_restart_file_scmmode(ncid, scm_restart_file)
    else
      scm_restart_file = .false.
    endif

    if (masterproc .and. debug > 1) then
       write(iulog,*) subname//' ',trim(flag),' ',trim(varname),' ',trim(elmlevel)
    end if

#if ({ITYPE}==TYPEDOUBLE) 
    if ( present(cnvrtnan2fill) )then
       if (.not. cnvrtnan2fill) then
          call shr_sys_abort(' ERROR: cnvrtnan2fill present but NOT set to true -- MUST set it to TRUE if used'//&
               errMsg(__FILE__, __LINE__))
       endif
    end if
#endif

    if (flag == 'read') then

       call ncd_inqvid(ncid, varname, varid, vardesc, readvar=varpresent)
       if (varpresent) then
          if ((scm_multcols .or. single_column) .and. .not. scm_restart_file) then
             start(:) = 1 ; count(:) = 1
             call scam_field_offsets(ncid,elmlevel,vardesc,start,count)
             if (present(nt)) then
                !-------------------------------------------------------------
                !Find the last non-zero occurrence in the "start" array
                !-------------------------------------------------------------

                !Assuming that the last non-zero index corresponds to the
                !time dimension. We need to read the "nt" th time dimension here

                n_time_ind = 3

                do while(start(n_time_ind) == 0)
                   n_time_ind = n_time_ind - 1 ! time index for the start array
                   if(n_time_ind == 1) exit    ! do not test for start(1)
                enddo
                start(n_time_ind) = nt ; count(n_time_ind) = 1
             endif

#if ({ITYPE}==TYPELOGICAL)
             allocate(idata(size(data))) 
             status = pio_get_var(ncid, varid, start, count, idata)
             data = (idata == 1)
             if ( any(idata /= 0 .and. idata /= 1) )then
                call shr_sys_abort(' ERROR: read in bad integer value(s) for logical data'//errMsg(__FILE__, __LINE__))
             end if
             deallocate( idata )
#else
             if (single_column) status=pio_get_var(ncid,varid,start,count,data)

             ! If in SCM multiple column mode, then extract data for the column closest
             !   to the IOP lat/lon from initial condition file and propogate it to
             !   all the columns in the domain.
             if (scm_multcols) then

               if (count(2) .gt. 0) then

                 start_iop2d(:) = start(1:2)
                 count_iop2d(:) = count(1:2)

                 allocate(iop_data2d(count(1),count(2)))
                 status = pio_get_var(ncid, varid, start_iop2d, count_iop2d, iop_data2d)
                 data(:) = iop_data2d(1,1)
                 deallocate(iop_data2d)

               else

                 start_iop1d(:) = start(1)
                 count_iop1d(:) = count(1)

                 allocate(iop_data1d(count(1)))
                 status = pio_get_var(ncid, varid, start_iop1d, count_iop1d, iop_data1d)
                 data(:) = iop_data1d(1)
                 deallocate(iop_data1d)

               endif

             endif ! scm_multcols
#endif
          else
             status = pio_inq_varndims(ncid, vardesc, ndims)
             status = pio_inq_vardimid(ncid, vardesc, dids)
             status = pio_inq_dimname(ncid,dids(ndims),dimname)
             if ('time' == trim(dimname)) then
                ndims_iod = ndims - 1
             else
                ndims_iod = ndims
             end if
             do n = 1,ndims_iod
                status = pio_inq_dimlen(ncid,dids(n),dims(n))
             enddo
#if ({ITYPE}==TYPELOGICAL)
             call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
                  PIO_INT, iodnum)
#else
             call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
                  PIO_{TYPE}, iodnum)
#endif
             iodesc_plus => iodesc_list(iodnum)
             if (present(nt)) then
                call pio_setframe(ncid,vardesc, int(nt,kind=PIO_OFFSET_KIND))
             end if
#if ({ITYPE}==TYPELOGICAL)
             allocate(idata(size(data))) 
             call pio_read_darray(ncid, vardesc, iodesc_plus%iodesc, idata, status)
             data = (idata == 1)
             if ( any(idata /= 0 .and. idata /= 1) )then
                call shr_sys_abort(' ERROR: read in bad integer value(s) for logical data'//errMsg(__FILE__, __LINE__))
             end if
             deallocate( idata )
#else
             call pio_read_darray(ncid, vardesc, iodesc_plus%iodesc, data, status)
#endif
          end if
       end if
       if (present(readvar)) readvar = varpresent

    elseif (flag == 'write') then

       call ncd_inqvid(ncid, varname ,varid, vardesc)
       status = pio_inq_varndims(ncid, vardesc, ndims)
       status = pio_inq_vardimid(ncid, vardesc, dids)
       status = pio_inq_dimname(ncid,dids(ndims),dimname)
       if ('time' == trim(dimname)) then
          ndims_iod = ndims - 1
       else
          ndims_iod = ndims
       end if
       do n = 1,ndims_iod
          status = pio_inq_dimlen(ncid,dids(n),dims(n))
       enddo
#if ({ITYPE}==TYPELOGICAL)
       call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
            PIO_INT, iodnum)
#else
       call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
            PIO_{TYPE}, iodnum)
#endif
       iodesc_plus => iodesc_list(iodnum)
       if (present(nt)) then
          call pio_setframe(ncid, vardesc, int(nt,kind=PIO_OFFSET_KIND))
       end if
#if ({ITYPE}==TYPELOGICAL)
       allocate( idata(size(data)) ) 
       where( data )
          idata = 1
       elsewhere
          idata = 0
       end where
       call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, idata, status, fillval=0)
       deallocate( idata )
#elif ({ITYPE}==TYPEINT)
       call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, data, status, fillval=0)
#elif ({ITYPE}==TYPEDOUBLE)
       call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, data, status, fillval=spval)
#endif
    else

       if (masterproc) then
          write(iulog,*) subname//' ERROR: unsupported flag ',trim(flag)
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       endif

    endif

  end subroutine ncd_io_1d_{TYPE}

  !-----------------------------------------------------------------------

  !TYPE int,double
  subroutine ncd_io_2d_{TYPE}(varname, data, dim1name, lowerb2, upperb2, &
       flag, ncid, nt, readvar, switchdim, cnvrtnan2fill)
    !
    ! !DESCRIPTION:
    ! Netcdf i/o of 2d 
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout) :: ncid         ! netcdf file id
    character(len=*) , intent(in)  :: flag            ! 'read' or 'write'
    character(len=*) , intent(in)  :: varname         ! variable name
    {VTYPE}          , pointer     :: data(:,:)       ! local decomposition input data
    character(len=*) , intent(in)  :: dim1name        ! dimension 1 name
    integer, optional, intent(in)  :: nt              ! time sample index
    integer, optional, intent(in)  :: lowerb2,upperb2 ! lower and upper bounds of second dimension
    logical, optional, intent(out) :: readvar         ! true => variable is on initial dataset (read only)
    logical, optional, intent(in)  :: switchdim       ! true=> permute dim1 and dim2 for output
    logical, optional, intent(in)  :: cnvrtnan2fill   ! true => convert any NaN's to _FillValue (spval)
    !
    ! !LOCAL VARIABLES:
#if ({ITYPE}==TYPEINT) 
    integer , pointer :: temp(:,:)
#else
    real(r8), pointer :: temp(:,:)
#endif

    {VTYPE}, allocatable :: iop_data2d(:,:)
    {VTYPE}, allocatable :: iop_data1d(:)

    integer           :: ndim1,ndim2       
    character(len=8)  :: elmlevel   ! elmlevel
    character(len=32) :: dimname    ! temporary      
    integer           :: status     ! error status
    integer           :: ndims      ! ndims total for var
    integer           :: ndims_iod  ! ndims iodesc for var
    integer           :: varid      ! varid
    integer           :: n,i,j      ! indices
    integer           :: n_time_ind !time slice index of the "start" array
    integer           :: dims(4)    ! dim sizes       
    integer           :: dids(4)    ! dim ids
    integer           :: iodnum     ! iodesc num in list
    integer           :: start(4)   ! netcdf start index
    integer           :: count(4)   ! netcdf count index
    integer           :: start_iop1d(1)   ! netcdf start index
    integer           :: count_iop1d(1)   ! netcdf count index
    integer           :: start_iop2d(2)   ! netcdf start index
    integer           :: count_iop2d(2)   ! netcdf count index
    logical           :: varpresent ! if true, variable is on tape
    logical           :: scm_restart_file ! restart file in SCM mode?
    integer           :: lb1,lb2
    integer           :: ub1,ub2
    integer           :: offset, ii, jj
    type(iodesc_plus_type) , pointer  :: iodesc_plus
    type(var_desc_t)                  :: vardesc
    character(len=*),parameter :: subname='ncd_io_2d_{TYPE}' ! subroutine name
    !-----------------------------------------------------------------------

    start(:)=0
    count(:)=0

    elmlevel = dim1name

    if (scm_multcols .or. single_column) then
      call is_restart_file_scmmode(ncid, scm_restart_file)
    else
      scm_restart_file = .false.
    endif

    if (masterproc .and. debug > 1) then
       write(iulog,*) trim(subname),' ',trim(flag),' ',trim(varname),' ',trim(elmlevel)
    end if

#if ({ITYPE}==TYPEDOUBLE) 
    if ( present(cnvrtnan2fill) )then
       if (.not. cnvrtnan2fill) then
          call shr_sys_abort( ' ERROR: cnvrtnan2fill present but NOT set to true -- MUST set it to TRUE if used'//&
               errMsg(__FILE__, __LINE__))
       endif
    end if
#endif

    lb1 = lbound(data, dim=1)
    ub1 = ubound(data, dim=1)
    lb2 = lbound(data, dim=2)
    ub2 = ubound(data, dim=2)

    if (present(switchdim)) then
       if (present(lowerb2)) lb2 = lowerb2
       if (present(upperb2)) ub2 = upperb2
       allocate(temp(lb2:ub2,lb1:ub1))
    end if

    if (flag == 'read') then

       call ncd_inqvid(ncid, varname, varid, vardesc, readvar=varpresent)
       if (varpresent) then
          if ((scm_multcols .or. single_column) .and. .not. scm_restart_file) then
             start(:) = 1 ; count(:) = 1
             call scam_field_offsets(ncid, elmlevel, vardesc, start, count)
                if (present(nt)) then

                   !-------------------------------------------------------------
                   !Find the last non-zero occurrence in the "start" array
                   !-------------------------------------------------------------

                   !Assuming that the last non-zero index corresponds to the
                   !time dimension. We need to read the "nt" th time dimension here

                   n_time_ind = 4
                   do while(start(n_time_ind) == 0)
                      n_time_ind = n_time_ind - 1 ! time index for the start array
                      if(n_time_ind == 1) exit    ! do not test for start(1)
                   enddo
                   start(n_time_ind) = nt ; count(n_time_ind) = 1
                endif
             if (present(switchdim) .and. .not. scm_multcols) then
                status = pio_get_var(ncid, vardesc, start, count, temp)
                do j = lb2,ub2
                   do i = lb1,ub1
                      data(i,j) = temp(j,i) 
                   end do
                end do
             else
                if (single_column) status = pio_get_var(ncid, vardesc, start, count, data)

                ! If in SCM multiple column mode, then extract data for the column closest
                !   to the IOP lat/lon from initial condition file and propogate it to
                !   all the columns in the domain.
                if (scm_multcols) then

                  if (count(2) .gt. 0) then

                    start_iop2d(:) = start(1:2)
                    count_iop2d(:) = count(1:2)

                    allocate(iop_data2d(count(1),count(2)))
                    status = pio_get_var(ncid, varid, start_iop2d, count_iop2d, iop_data2d)
                    data(:,:) = iop_data2d(1,1)
                    deallocate(iop_data2d)

                  else
                    start_iop1d(:) = start(1)
                    count_iop1d(:) = count(1)

                    allocate(iop_data1d(count(1)))
                    status = pio_get_var(ncid, varid, start_iop1d, count_iop1d, iop_data1d)
                    data(:,:) = iop_data1d(1)
                    deallocate(iop_data1d)

                  endif

                endif ! if scm_multcols

             endif ! if scm_multcols or single_column
          else
             status = pio_inq_varndims(ncid, vardesc, ndims)
             status = pio_inq_vardimid(ncid,vardesc, dids)
             status = pio_inq_dimname(ncid, dids(ndims), dimname)
             if (ndims == 0) then
                write(iulog,*) trim(subname),' ERROR: ndims must be greater than 0'
                call shr_sys_abort(errMsg(__FILE__, __LINE__))
             end if
             if ('time' == trim(dimname)) then
                ndims_iod = ndims - 1
             else
                ndims_iod = ndims
             end if
             do n = 1,ndims_iod
                status = pio_inq_dimlen(ncid,dids(n),dims(n))
             enddo
             if (present(switchdim)) then
                call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
                     PIO_{TYPE}, iodnum, switchdim=.true.)
             else
                call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
                     PIO_{TYPE}, iodnum)
             end if
             iodesc_plus => iodesc_list(iodnum)
             if (present(nt)) then
                call pio_setframe(ncid, vardesc, int(nt,kind=PIO_OFFSET_KIND))
             end if
             if (present(switchdim)) then
                call pio_read_darray(ncid, vardesc, iodesc_plus%iodesc, temp, status)
                do j = lb2,ub2
                do i = lb1,ub1
                   data(i,j) = temp(j,i) 
                end do
                end do
             else
                call pio_read_darray(ncid, vardesc, iodesc_plus%iodesc, data, status)
             end if
          end if
#if ({ITYPE}!=TYPEINT)
          if ( present(cnvrtnan2fill) )then
             do j = lb2,ub2
             do i = lb1,ub1
                if ( data(i,j) == spval )then
                   data(i,j) = nan
                end if
             end do
             end do
          end if
#endif
       end if
       if (present(readvar)) readvar = varpresent

    else if (flag == 'write') then

       call ncd_inqvid(ncid, varname ,varid, vardesc)
       status = pio_inq_varndims(ncid, vardesc, ndims)
       status = pio_inq_vardimid(ncid, vardesc , dids)
       if (ndims == 0) then
          write(iulog,*) trim(subname),' ERROR: ndims must be greater than 0'
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       end if
       status = pio_inq_dimname(ncid,dids(ndims),dimname)
       if ('time' == trim(dimname)) then
          ndims_iod = ndims - 1
       else
          ndims_iod = ndims
       end if
       do n = 1,ndims_iod
          status = pio_inq_dimlen(ncid,dids(n),dims(n))
       enddo
       if (present(switchdim)) then
          call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
               PIO_{TYPE}, iodnum, switchdim=.true.)
       else
          call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
               PIO_{TYPE}, iodnum)
       end if
       iodesc_plus => iodesc_list(iodnum)
       if (present(nt)) then
          call pio_setframe(ncid, vardesc, int(nt,kind=PIO_OFFSET_KIND))
       end if
       if (present(switchdim)) then
          do j = lb2,ub2
             do i = lb1,ub1
                temp(j,i) = data(i,j)
             end do
          end do
       end if
#if ({ITYPE}==TYPEINT)
       if (present(switchdim)) then
          call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, temp, status, fillval=0)
       else
          call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, data, status, fillval=0)
       end if
#else
       if (present(switchdim)) then
          call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, temp, status, fillval=spval)
       else
          call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, data, status, fillval=spval)
       end if
       if ( present(cnvrtnan2fill) )then
          do j = lb2,ub2
             do i = lb1,ub1
                if ( isnan(data(i,j)) )then
                   data(i,j) = spval
                end if
             end do
          end do
       end if
#endif

    else

       if (masterproc) then
          write(iulog,*) subname,' error: unsupported flag ',trim(flag)
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       endif

    endif

    if (present(switchdim)) then
       deallocate(temp)
    end if

  end subroutine ncd_io_2d_{TYPE}

  !-----------------------------------------------------------------------

  !TYPE int,double
  subroutine ncd_io_3d_{TYPE}(varname, data, dim1name, flag, ncid, nt, readvar)
    !
    ! !DESCRIPTION:
    ! Netcdf i/o of 3d 
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout) :: ncid         ! netcdf file id
    character(len=*) , intent(in)  :: flag            ! 'read' or 'write'
    character(len=*) , intent(in)  :: varname         ! variable name
    {VTYPE}          , pointer     :: data(:,:,:)     ! local decomposition input data
    character(len=*) , intent(in)  :: dim1name        ! dimension 1 name
    integer, optional, intent(in)  :: nt              ! time sample index
    logical, optional, intent(out) :: readvar         ! true => variable is on initial dataset (read only)
    !

    {VTYPE}, allocatable :: iop_data3d(:,:,:)
    {VTYPE}, allocatable :: iop_data2d(:,:)

    ! !LOCAL VARIABLES:
    integer, parameter               :: max_ind = 5   !index for arrays which read data from the netcdf files
    integer                          :: n_time_ind !time slice index of the "start" array
    integer                          :: i, j, ii, jj
    integer                          :: ndim1,ndim2
    character(len=8)                 :: elmlevel   ! elmlevel
    character(len=32)                :: dimname    ! temporary      
    integer                          :: status     ! error status
    integer                          :: ndims      ! ndims total for var
    integer                          :: ndims_iod  ! ndims iodesc for var
    integer                          :: varid      ! varid
    integer                          :: ind          ! index      
    integer                          :: dims(max_ind)    ! dim sizes       
    integer                          :: dids(max_ind)    ! dim ids
    integer                          :: iodnum     ! iodesc num in list
    integer                          :: start(max_ind)   ! netcdf start index
    integer                          :: count(max_ind)   ! netcdf count index
    integer                          :: start_iop2d(2)   ! netcdf start index
    integer                          :: count_iop2d(2)   ! netcdf count index
    integer                          :: start_iop3d(3)   ! netcdf start index
    integer                          :: count_iop3d(3)   ! netcdf count index
    logical                          :: varpresent ! if true, variable is on tape
    logical                          :: scm_restart_file ! restart file in SCM mode
    type(iodesc_plus_type) , pointer :: iodesc_plus
    type(var_desc_t)                 :: vardesc
    character(len=*),parameter :: subname='ncd_io_3d_{TYPE}' ! subroutine name
    !-----------------------------------------------------------------------

    elmlevel = dim1name

    if (scm_multcols .or. single_column) then
      call is_restart_file_scmmode(ncid, scm_restart_file)
    else
      scm_restart_file = .false.
    endif

    if (masterproc .and. debug > 1) then
       write(iulog,*) trim(subname),' ',trim(flag),' ',trim(varname),' ',trim(elmlevel)
    end if

    if (flag == 'read') then

       call ncd_inqvid(ncid, varname, varid, vardesc, readvar=varpresent)
       if (varpresent) then
          if ((scm_multcols .or. single_column) .and. .not. scm_restart_file) then
             start(:) = 1
             count(:) = 1
             call scam_field_offsets(ncid, elmlevel, vardesc, start, count)
             if (present(nt)) then
                !-------------------------------------------------------------
                !Find the last  non-zero occurrence in the "start" array
                !-------------------------------------------------------------

                !Assuming that the last non-zero index corresponds to the
                !time dimension. We need to read the "nt" th time dimension here

                n_time_ind = max_ind
                do while(start(n_time_ind) == 0)
                   n_time_ind = n_time_ind -1 ! time index for the start array
                   if(n_time_ind == 1) exit   ! do not test for start(1)
                enddo
                start(n_time_ind) = nt ; count(n_time_ind) = 1

             endif

             if (single_column) status = pio_get_var(ncid, vardesc, start, count, data)

             ! If in SCM multiple column mode, then extract data for the column closest
             !   to the IOP lat/lon from initial condition file and propogate it to
             !   all the columns in the domain.
             if (scm_multcols) then

               if (start(3) .eq. 0) then

                 count_iop2d(1:2) = count(1:2)
                 start_iop2d(1:2) = start(1:2)

                 allocate(iop_data2d(count_iop2d(1),count_iop2d(2)))

                 status = pio_get_var(ncid, vardesc, start_iop2d, count_iop2d, iop_data2d)

                 do i=lbound(data,dim=1),ubound(data,dim=1)
                   jj=1
                   do j=lbound(data,dim=3),ubound(data,dim=3)
                     data(i,:,j) = iop_data2d(1,jj)
                     jj=jj+1
                   enddo
                 enddo

                 deallocate(iop_data2d)

               else

                 count_iop3d(1:3) = count(1:3)
                 start_iop3d(1:3) = start(1:3)

                 allocate(iop_data3d(count_iop3d(1),count_iop3d(2),count_iop3d(3)))

                 status = pio_get_var(ncid, vardesc, start_iop3d, count_iop3d, iop_data3d)

                 do i=lbound(data,dim=1),ubound(data,dim=1)
                   jj=1
                   do j=lbound(data,dim=3),ubound(data,dim=3)
                     data(i,:,j) = iop_data3d(1,jj,1)
                     jj=jj+1
                   enddo
                 enddo

                 deallocate(iop_data3d)

               endif

            endif ! scm_multcols

          else
             status = pio_inq_varndims(ncid, vardesc, ndims)
             status = pio_inq_vardimid(ncid,vardesc, dids)
             status = pio_inq_dimname(ncid, dids(ndims), dimname)
             if (ndims == 0) then
                write(iulog,*) trim(subname),' ERROR: ndims must be greater than 0'
                call shr_sys_abort(errMsg(__FILE__, __LINE__))
             end if
             if ('time' == trim(dimname)) then
                ndims_iod = ndims - 1
             else
                ndims_iod = ndims
             end if
             do ind = 1,ndims_iod
                status = pio_inq_dimlen(ncid,dids(ind),dims(ind))
             enddo
             call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
                  PIO_{TYPE}, iodnum)
             iodesc_plus => iodesc_list(iodnum)
             if (present(nt)) then
                call pio_setframe(ncid, vardesc, int(nt,kind=PIO_OFFSET_KIND))
             end if
             call pio_read_darray(ncid, vardesc, iodesc_plus%iodesc, data, status)
          end if
       end if
       if (present(readvar)) readvar = varpresent

    else if (flag == 'write') then

       call ncd_inqvid(ncid, varname ,varid, vardesc)
       status = pio_inq_varndims(ncid, vardesc, ndims)
       status = pio_inq_vardimid(ncid, vardesc , dids)
       if (ndims == 0) then
          write(iulog,*) trim(subname),' ERROR: ndims must be greater than 0'
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       end if
       status = pio_inq_dimname(ncid,dids(ndims),dimname)
       if ('time' == trim(dimname)) then
          ndims_iod = ndims - 1
       else
          ndims_iod = ndims
       end if
       do ind = 1,ndims_iod
          status = pio_inq_dimlen(ncid,dids(ind),dims(ind))
       enddo
       call ncd_getiodesc(ncid, elmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
            PIO_{TYPE}, iodnum)
       iodesc_plus => iodesc_list(iodnum)
       if (present(nt)) then
          call pio_setframe(ncid, vardesc, int(nt,kind=PIO_OFFSET_KIND))
       end if
       call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, data, status)

    else

       if (masterproc) then
          write(iulog,*) subname,' error: unsupported flag ',trim(flag)
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       endif

    endif

  end subroutine ncd_io_3d_{TYPE}

  !-----------------------------------------------------------------------

  !TYPE int,double
  subroutine ncd_io_4d_{TYPE}(varname, data, dim1name, flag, ncid, nt, readvar)
    !
    ! !DESCRIPTION:
    ! Netcdf i/o of 4d 
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout) :: ncid         ! netcdf file id
    character(len=*) , intent(in)  :: flag            ! 'read' or 'write'
    character(len=*) , intent(in)  :: varname         ! variable name
    {VTYPE}          , pointer     :: data(:,:,:,:)     ! local decomposition input data
    character(len=*) , intent(in)  :: dim1name        ! dimension 1 name
    integer, optional, intent(in)  :: nt              ! time sample index
    logical, optional, intent(out) :: readvar         ! true => variable is on initial dataset (read only)

    {VTYPE}, allocatable :: iop_data4d(:,:,:,:)
    {VTYPE}, allocatable :: iop_data3d(:,:,:)

    !
    ! !LOCAL VARIABLES:
    integer                          :: i, j
    integer                          :: ndim1,ndim2       
    character(len=8)                 :: clmlevel   ! clmlevel
    character(len=32)                :: dimname    ! temporary      
    integer                          :: status     ! error status
    integer                          :: ndims      ! ndims total for var
    integer                          :: ndims_iod  ! ndims iodesc for var
    integer                          :: varid      ! varid
    integer                          :: n          ! index      
    integer                          :: n_time_ind !time slice index of the "start" array
    integer                          :: dims(5)    ! dim sizes       
    integer                          :: dids(5)    ! dim ids
    integer                          :: iodnum     ! iodesc num in list
    integer                          :: start(6)   ! netcdf start index
    integer                          :: count(6)   ! netcdf count index
    integer                          :: start_iop3d(3)   ! netcdf start index
    integer                          :: count_iop3d(3)   ! netcdf count index
    integer                          :: start_iop4d(4)   ! netcdf start index
    integer                          :: count_iop4d(4)   ! netcdf count index
    logical                          :: varpresent ! if true, variable is on tape
    logical                          :: scm_restart_file ! restart file in SCM mode
    type(iodesc_plus_type) , pointer :: iodesc_plus
    type(var_desc_t)                 :: vardesc
    character(len=*),parameter :: subname='ncd_io_4d_{TYPE}' ! subroutine name
    !-----------------------------------------------------------------------

    clmlevel = dim1name

    if (scm_multcols .or. single_column) then
      call is_restart_file_scmmode(ncid, scm_restart_file)
    else
      scm_restart_file = .false.
    endif

    if (masterproc .and. debug > 1) then
       write(iulog,*) trim(subname),' ',trim(flag),' ',trim(varname),' ',trim(clmlevel)
    end if

    if (flag == 'read') then
       call ncd_inqvid(ncid, varname, varid, vardesc, readvar=varpresent)
       if (varpresent) then
          if ((scm_multcols .or. single_column) .and. .not. scm_restart_file) then
             start(:) = 1; count(:) = 1
             call scam_field_offsets(ncid, clmlevel, vardesc, start, count)
             if (present(nt)) then

                !-------------------------------------------------------------
                !Find the last non-zero occurrence in the "start" array
                !-------------------------------------------------------------

                !Assuming that the last non-zero index corresponds to the
                !time dimension. We need to read the "nt" th time dimension here

                n_time_ind = 6
                do while(start(n_time_ind) == 0)
                   n_time_ind = n_time_ind -1 ! time index for the start array
                   if(n_time_ind == 1) exit   ! do not test for start(1)
                enddo
                start(n_time_ind) = nt ; count(n_time_ind) = 1
             endif

             if (single_column) status = pio_get_var(ncid, vardesc, start, count, data)
             if (scm_multcols) then
               if (start(4) .eq. 0) then

                 count_iop3d(1:3) = count(1:3)
                 start_iop3d(1:3) = start(1:3)

                 allocate(iop_data3d(count_iop3d(1),count_iop3d(2),count_iop3d(3)))

                 status = pio_get_var(ncid, vardesc, start_iop3d, count_iop3d, iop_data3d)

                 do i=1,count_iop3d(2)
                   do j=1,count_iop3d(3)
                      data(:,:,i,j) = iop_data3d(1,i,j)
                   enddo
                enddo

                 deallocate(iop_data3d)
               else
                 count_iop4d(1:4) = count(1:4)
                 start_iop4d(1:4) = start(1:4)

                 allocate(iop_data4d(count_iop4d(1),count_iop4d(2),count_iop4d(3),count_iop4d(4)))

                 status = pio_get_var(ncid, vardesc, start_iop4d, count_iop4d, iop_data4d)

                 do i=1,count_iop4d(2)
                   do j=1,count_iop4d(3)
                     data(:,:,i,j) = iop_data4d(1,i,j,1)
                   enddo
                 enddo

                 deallocate(iop_data4d)

               endif

            endif ! scm_multcols

          else
             status = pio_inq_varndims(ncid, vardesc, ndims)
             status = pio_inq_vardimid(ncid,vardesc, dids)
             status = pio_inq_dimname(ncid, dids(ndims), dimname)
             if (ndims == 0) then
                write(iulog,*) trim(subname),' ERROR: ndims must be greater than 0'
                call shr_sys_abort(errMsg(__FILE__, __LINE__))
             end if
             if ('time' == trim(dimname)) then
                ndims_iod = ndims - 1
             else
                ndims_iod = ndims
             end if
             do n = 1,ndims_iod
                status = pio_inq_dimlen(ncid,dids(n),dims(n))
             enddo
             call ncd_getiodesc(ncid, clmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
                  PIO_{TYPE}, iodnum)
             iodesc_plus => iodesc_list(iodnum)
             if (present(nt)) then
                call pio_setframe(ncid, vardesc, int(nt,kind=PIO_OFFSET_KIND))
             end if
             call pio_read_darray(ncid, vardesc, iodesc_plus%iodesc, data, status)
          end if
       end if
       if (present(readvar)) readvar = varpresent
    else if (flag == 'write') then

       call ncd_inqvid(ncid, varname ,varid, vardesc)
       status = pio_inq_varndims(ncid, vardesc, ndims)
       status = pio_inq_vardimid(ncid, vardesc , dids)
       if (ndims == 0) then
          write(iulog,*) trim(subname),' ERROR: ndims must be greater than 0'
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       end if
       status = pio_inq_dimname(ncid,dids(ndims),dimname)
       if ('time' == trim(dimname)) then
          ndims_iod = ndims - 1
       else
          ndims_iod = ndims
       end if
       do n = 1,ndims_iod
          status = pio_inq_dimlen(ncid,dids(n),dims(n))
       enddo
       call ncd_getiodesc(ncid, clmlevel, ndims_iod, dims(1:ndims_iod), dids(1:ndims_iod), &
            PIO_{TYPE}, iodnum)
       iodesc_plus => iodesc_list(iodnum)
       if (present(nt)) then
          call pio_setframe(ncid, vardesc, int(nt,kind=PIO_OFFSET_KIND))
       end if
       call pio_write_darray(ncid, vardesc, iodesc_plus%iodesc, data, status)

    else

       if (masterproc) then
          write(iulog,*) subname,' error: unsupported flag ',trim(flag)
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       endif

    endif

  end subroutine ncd_io_4d_{TYPE}

  !------------------------------------------------------------------------

  subroutine scam_field_offsets( ncid, dim1name, vardesc, start, count, &
       found, posNOTonfile)
    !
    ! !DESCRIPTION: 
    ! Read/Write initial data from/to netCDF instantaneous initial data file 
    !
    ! !USES:
    use elm_varctl,     only: scmlon,scmlat,single_column
    use shr_scam_mod,   only: shr_scam_getCloseLatLon
    use shr_string_mod, only: shr_string_toLower
    !
    ! !ARGUMENTS:
    class(file_desc_t), intent(inout) :: ncid         ! netcdf file id
    character(len=*)  , intent(in)    :: dim1name     ! dimension 1 name
    type(Var_desc_t)  , intent(inout) :: vardesc      ! variable descriptor
    integer           , intent(inout) :: start(:)     ! start index
    integer           , intent(inout) :: count(:)     ! count to retrieve
    logical, optional , intent(out)   :: found        ! if present return true if found
    ! dimensions on file else false if NOT present abort if can't find
    logical, optional , intent(in)    :: posNOTonfile ! Position is NOT on this file
    !
    ! !LOCAL VARIABLES:
    integer              :: cc,i,ii           ! index variable
    integer              :: data_offset       ! offset into land array 1st column 
    integer              :: ndata             ! number of column (or pft points to read)
    real(r8) , pointer   :: cols1dlon(:)      ! holds cols1d_ixy var
    real(r8) , pointer   :: cols1dlat(:)      ! holds cols1d_jxy var
    real(r8) , pointer   :: pfts1dlon(:)      ! holds pfts1d_ixy var
    real(r8) , pointer   :: pfts1dlat(:)      ! holds pfts1d_jxy var
    real(r8) , pointer   :: land1dlon(:)      ! holds land1d_ixy var
    real(r8) , pointer   :: land1dlat(:)      ! holds land1d_jxy var
    real(r8) , pointer   :: topo1dlon(:)      ! holds topo1d_ixy var
    real(r8) , pointer   :: topo1dlat(:)      ! holds topo1d_ixy var
    integer, allocatable :: cols(:)           ! grid cell columns for scam
    integer, allocatable :: pfts(:)           ! grid cell pfts for scam
    integer, allocatable :: landunits(:)      ! grid cell landunits for scam
    integer, allocatable :: topounits(:)      ! grid cell topounits for scam
    integer, allocatable :: dids(:)           ! dim ids
    integer              :: varid             ! netCDF variable id
    integer              :: status            ! return code
    integer              :: latidx,lonidx     ! latitude/longitude indices
    real(r8)             :: closelat,closelon ! closest latitude and longitude indices
    integer              :: ndims,dimlen      ! number of dimensions in desired variable
    character(len=32)    :: dimname           ! dimension name
    character(len=32)    :: subname = 'scam_field_offsets'
    !------------------------------------------------------------------------

    start(:)=0
    count(:)=0

    if ( present(posNOTonfile) )then
       if ( posNOTonfile )then
          if ( .not. present(found) )then
             call shr_sys_abort('ERROR: Bad subroutine calling structure posNOTonfile sent, but found was NOT!'//&
                  errMsg(__FILE__, __LINE__))
          end if
          found = .false.
          return
       end if
    end if

    ! find closest land grid cell for this point

    if ( present(found) )then
       call shr_scam_getCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,latidx,lonidx,found)
       if ( .not. found ) return
    else
       call shr_scam_getCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
    end if

    call ncd_inqvdims(ncid,ndims,vardesc)

    allocate(dids(ndims))
    status = pio_inq_vardimid(ncid, vardesc, dids)
    do i = 1,ndims
       status = pio_inq_dimname(ncid,dids(i),dimname)
       dimname=shr_string_toLower(dimname)
       status = pio_inq_dimlen(ncid,dids(i),dimlen)
       if ( trim(dimname)=='nj'.or. trim(dimname)=='lat'.or. trim(dimname)=='lsmlat') then
          start(i)=latidx
          count(i)=1
       !else if ( trim(dimname)=='ni'.or. trim(dimname)=='lon'.or. trim(dimname)=='lsmlon' .or. trim(dimname)=='gridcell' .or. trim(dimname)=='topounit') then
       else if ( trim(dimname)=='ni'.or. trim(dimname)=='lon'.or. trim(dimname)=='lsmlon' .or. trim(dimname)=='gridcell') then   
          start(i)=lonidx
          count(i)=1
       else if ( trim(dimname)=='column') then

          allocate (cols1dlon(dimlen))
          allocate (cols1dlat(dimlen))
          allocate (cols(dimlen))

          status = pio_inq_varid(ncid, 'cols1d_lon', varid)
          status = pio_get_var(ncid, varid, cols1dlon)
          status = pio_inq_varid(ncid, 'cols1d_lat', varid)
          status = pio_get_var(ncid, varid, cols1dlat)

          cols(:)     = huge(1)
          data_offset = huge(1)
          ii = 1
          ndata = 0
          do cc = 1, dimlen
             if (cols1dlon(cc) == closelon.and.cols1dlat(cc) == closelat) then
                cols(ii)=cc
                ndata  =ii
                ii=ii+1
             end if
          end do
          if (ndata == 0) then
             write(iulog,*)'couldnt find any columns for this latitude ',latidx,' and longitude ',lonidx
             call shr_sys_abort('ERROR:: no columns for this position'//errMsg(__FILE__, __LINE__))
          else
             data_offset=cols(1)
          end if

          deallocate (cols1dlon)
          deallocate (cols1dlat)
          deallocate (cols)

          start(i) = data_offset
          count(i) = ndata
       else if ( trim(dimname)=='pft') then

          allocate (pfts1dlon(dimlen))
          allocate (pfts1dlat(dimlen))
          allocate (pfts(dimlen))

          status = pio_inq_varid(ncid, 'pfts1d_lon', varid)
          status = pio_get_var(ncid, varid, pfts1dlon)

          status = pio_inq_varid(ncid, 'pfts1d_lat', varid)
          status = pio_get_var(ncid, varid, pfts1dlat)

          pfts(:)     = huge(1)
          data_offset = huge(1)
          ii     = 1
          ndata = 0
          do cc = 1, dimlen
             if (pfts1dlon(cc) == closelon.and.pfts1dlat(cc) == closelat) then
                pfts(ii)=cc
                ndata  =ii
                ii=ii+1
             end if
          end do
          if (ndata == 0) then
             write(iulog,*)'couldnt find any pfts for this latitude ',closelat,' and longitude ',closelon
             call shr_sys_abort('ERROR:: no PFTs for this position'//errMsg(__FILE__, __LINE__))
          else
             data_offset=pfts(1)
          end if

          deallocate (pfts1dlon)
          deallocate (pfts1dlat)
          deallocate (pfts)

          start(i) = data_offset
          count(i) = ndata
       else if ( trim(dimname)=='landunit') then

          allocate (land1dlon(dimlen))
          allocate (land1dlat(dimlen))
          allocate (landunits(dimlen))

          status = pio_inq_varid(ncid, 'land1d_lon', varid)
          status = pio_get_var(ncid, varid, land1dlon)

          status = pio_inq_varid(ncid, 'land1d_lat', varid)
          status = pio_get_var(ncid, varid, land1dlat)

          landunits(:) = huge(1)
          data_offset  = huge(1)
          ii     = 1
          ndata = 0
          do cc = 1, dimlen
             if (land1dlon(cc) == closelon.and.land1dlat(cc) == closelat) then
                landunits(ii)=cc
                ndata  =ii
                ii=ii+1
             end if
          end do
          if (ndata == 0) then
             write(iulog,*)'couldnt find any landunits for this latitude ',closelat,' and longitude ',closelon
             call shr_sys_abort('ERROR:: no landunits for this position'//errMsg(__FILE__, __LINE__))
          else
             data_offset=landunits(1)
          end if

          deallocate (land1dlon)
          deallocate (land1dlat)
          deallocate (landunits)

          start(i) = data_offset
          count(i) = ndata
       else if ( trim(dimname)=='topounit') then

          allocate (topo1dlon(dimlen))
          allocate (topo1dlat(dimlen))
          allocate (topounits(dimlen))

          status = pio_inq_varid(ncid, 'topo1d_lon', varid)
          status = pio_get_var(ncid, varid, topo1dlon)

          status = pio_inq_varid(ncid, 'topo1d_lat', varid)
          status = pio_get_var(ncid, varid, topo1dlat)

          topounits(:) = huge(1)
          data_offset  = huge(1)
          ii     = 1
          ndata = 0
          do cc = 1, dimlen
             if (topo1dlon(cc) == closelon.and.topo1dlat(cc) == closelat) then
                topounits(ii)=cc
                ndata  =ii
                ii=ii+1
             end if
          end do
          if (ndata == 0) then
             write(iulog,*)'couldnt find any topounits for this latitude',closelat,' and longitude ',closelon
             call shr_sys_abort('ERROR:: no topounits for this position'//errMsg(__FILE__, __LINE__))
          else
             data_offset=topounits(1)
          end if

          deallocate (topo1dlon)
          deallocate (topo1dlat)
          deallocate (topounits)

          start(i) = data_offset
          count(i) = ndata
       else
          start(i)=1
          count(i)=dimlen
       end if
    enddo

    deallocate(dids)

  end subroutine scam_field_offsets

  !------------------------------------------------------------------------

  subroutine is_restart_file_scmmode(ncid, restart_file)

  implicit none

  ! DESCRIPTION:
  ! determines if the input file is a restart file.  This is important to determine
  !  for doubly-periodic CRM mode (SCM mult columns) because we need to know if we are
  !  calling the standard PIO routines (for a restart file) or the specialized ones
  !  that are used for SCM mode.  Note that this routine is specific for determining
  !  restarts ONLY for DP CRM mode.

  ! ARGUMENTS:
  class(file_desc_t), intent(inout) :: ncid         ! netcdf file id
  logical, intent(out) :: restart_file

  ! LOCAL VARIABLES:
  integer :: ni, ier1, ier2, dimid
  character (len = 80) :: title
  !------------------------------------------------------------------------

  ! Initialize value to default
  restart_file = .false.

  ! Check to see if one of the dimensions is "gridcell"
  ier1 = pio_inq_dimid(ncid, 'gridcell', dimid)
  ! Check to see if the attribute "title" exists
  ier2 = pio_get_att(ncid, pio_global, 'title', title)

  ! A restart file needs to pass the above to requirements
  if (ier1 == PIO_NOERR .and. ier2 == PIO_NOERR) then
    ! Determine the number of columns in the restart file
    ier1 = pio_inq_dimlen(ncid, dimid, ni)
    ! Get the title from the attribute
    ier2 = pio_get_att(ncid, pio_global, 'title', title)

    ! Final check, for DP-CRM restart file then the number of columns should be
    !  identical to the number of SCM columns.  If not, this file is not a restart
    !  file.  For protection against happenstance, also make sure the attribute
    !  title is that which is typically wrote out for restarts.
    if (ni == scm_nx .and. (title == "ELM Restart information" .or. title == &
      "ELM Restart History information, required to continue a simulation")) then

      restart_file = .true.
    endif
  endif

  end subroutine is_restart_file_scmmode

  !------------------------------------------------------------------------

  subroutine ncd_getiodesc(ncid, elmlevel, ndims, dims, dimids, &
       xtype, iodnum, switchdim) 
    !
    ! !DESCRIPTION: 
    ! Returns an index to an io descriptor
    !
    ! !ARGUMENTS:
    class(file_desc_t) , intent(inout) :: ncid       ! PIO file descriptor
    character(len=8)   , intent(in)    :: elmlevel   ! elmlevel
    integer            , intent(in)    :: ndims      ! ndims for var      
    integer            , intent(in)    :: dims(:)    ! dim sizes
    integer            , intent(in)    :: dimids(:)  ! dim ids
    integer            , intent(in)    :: xtype      ! file external type
    integer            , intent(out)   :: iodnum     ! iodesc num in list
    logical,optional   , intent(in)    :: switchdim  ! switch level dimension and first dim 
    !
    ! !LOCAL VARIABLES:
    integer   :: k,m,n,cnt                     ! indices
    integer   :: basetype                      ! pio basetype
    integer   :: gsmap_lsize                   ! local size of gsmap
    integer*8 :: gsmap_gsize                   ! global size of gsmap
    integer*8 :: fullsize                      ! size of entire array on cdf
    integer*8 :: gsize                         ! global size of elmlevel
    integer*8 :: vsize                         ! other dimensions
    integer   :: vsize1, vsize2                ! other dimensions
    integer   :: status                        ! error status
    logical   :: found                         ! true => found created iodescriptor
    integer   :: ndims_file                    ! temporary
    character(len=64) dimname_file           ! dimension name on file
    character(len=64) dimname_iodesc         ! dimension name from io descriptor
    type(mct_gsMap),pointer       :: gsmap   ! global seg map
    integer, pointer,dimension(:) :: gsmOP   ! gsmap ordered points
    integer*8, pointer  :: compDOF(:)
    character(len=32) :: subname = 'ncd_getiodesc'
    !------------------------------------------------------------------------

    ! Determining if need to create a new io descriptor

    n = 1
    found = .false.
    do while (n <= num_iodesc .and. .not.found)
       if (ndims == iodesc_list(n)%ndims .and. xtype == iodesc_list(n)%type) then
          found = .true.
          ! First found implies that dimension sizes are the same 
          do m = 1,ndims
             if (dims(m) /= iodesc_list(n)%dims(m)) then
                found = .false.
             endif
          enddo
          ! If found - then also check that dimension names are equal - 
          ! dimension ids in iodescriptor are only used to query dimension
          ! names associated with that iodescriptor
          if (found) then
             do m = 1,ndims
                status = PIO_inq_dimname(ncid,dimids(m),dimname_file)
                status = PIO_inquire(ncid, ndimensions=ndims_file)
                if (iodesc_list(n)%dimids(m) > ndims_file) then 
                   found = .false.
                   exit
                else
                   status = PIO_inq_dimname(ncid,iodesc_list(n)%dimids(m),dimname_iodesc)
                   if (trim(dimname_file) /= trim(dimname_iodesc)) then
                      found = .false.
                      exit
                   end if
                end if
             end do
          end if
          if (found) then
             iodnum = n
             if (iodnum > num_iodesc) then
                write(iulog,*) trim(subname),' ERROR: iodnum out of range ',iodnum,num_iodesc
                call shr_sys_abort(errMsg(__FILE__, __LINE__))
             endif
             RETURN
          endif
       endif
       n = n + 1
    enddo

    ! Creating a new io descriptor

    if (ndims > 0) then 
       num_iodesc = num_iodesc + 1
       if (num_iodesc > max_iodesc) then
          write(iulog,*) trim(subname),' ERROR num_iodesc gt max_iodesc ',max_iodesc
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       endif
       iodnum = num_iodesc
       if (masterproc .and. debug > 1) then
          write(iulog,*) trim(subname),' creating iodesc at iodnum,ndims,dims(1:ndims),xtype',&
               iodnum,ndims,dims(1:ndims),xtype
       endif
    end if

    if (xtype == pio_double ) then
       basetype = PIO_DOUBLE
    else if (xtype == pio_real) then
       basetype  = PIO_DOUBLE
    else if (xtype == pio_int) then
       basetype = PIO_INT
    end if

    call get_elmlevel_gsmap(elmlevel,gsmap)
    gsize = get_elmlevel_gsize(elmlevel)
    gsmap_lsize = mct_gsmap_lsize(gsmap,mpicom)
    gsmap_gsize = mct_gsmap_gsize(gsmap)

    call mct_gsMap_orderedPoints(gsmap,iam,gsmOP)

    fullsize = 1
    do n = 1,ndims
       fullsize = fullsize*dims(n)
    enddo

    vsize = fullsize / gsize
    if (mod(fullsize,gsize) /= 0) then
       write(iulog,*) subname,' ERROR in vsize ',fullsize,gsize,vsize,ndims,dims
       call shr_sys_abort(errMsg(__FILE__, __LINE__))
    endif

    allocate(compDOF(gsmap_lsize*vsize))

    if (present(switchdim)) then
       if (switchdim) then
          cnt = 0
          do m = 1,gsmap_lsize
             do n = 1,vsize
                cnt = cnt + 1
                compDOF(cnt) = (gsmOP(m)-1)*vsize + n
             enddo
          enddo
       else
          write(iulog,*) subname,' ERROR switch dims present must have switchdim true' 
          call shr_sys_abort(errMsg(__FILE__, __LINE__))
       end if
    else         ! currently allow for up to two vertical dimensions
       if (vsize /= 1 .and. vsize /= dims(ndims)) then
          vsize1 = vsize/dims(ndims)
          vsize2 = dims(ndims)
          if (vsize1*vsize2 /= vsize) then
             write(iulog,*)'vsize1= ',vsize1,' vsize2= ',vsize2,' vsize= ',vsize
             call shr_sys_abort('error in vsize1 and vsize2 computation'//errMsg(__FILE__, __LINE__))
          end if
          cnt = 0
          do k = 1,vsize2
             do n = 1,vsize1
                do m = 1,gsmap_lsize
                   cnt = cnt + 1
                   compDOF(cnt) = (k-1)*vsize1*gsmap_gsize + (n-1)*gsmap_gsize +  gsmOP(m) 
                enddo
             enddo
          end do
       else
          cnt = 0
          do n = 1,vsize
             do m = 1,gsmap_lsize
                cnt = cnt + 1
                compDOF(cnt) = (n-1)*gsmap_gsize + gsmOP(m)
             enddo
          enddo
       end if
    end if

    if (debug > 1) then
       do m = 0,npes-1
          if (iam == m) then
             write(iulog,*) trim(subname),' sizes1  = ',iam,gsize,gsmap_gsize,gsmap_lsize
             write(iulog,*) trim(subname),' sizes2  = ',iam,fullsize,npes,vsize
             write(iulog,*) trim(subname),' compDOF = ',iam,size(compDOF),minval(compDOF),maxval(compDOF)
             call shr_sys_flush(iulog)
          endif
          call mpi_barrier(mpicom,status)
       enddo
    endif

    deallocate(gsmOP)

    call pio_initdecomp(pio_subsystem, baseTYPE, dims(1:ndims), compDOF, iodesc_list(iodnum)%iodesc)

    deallocate(compDOF)

    iodesc_list(iodnum)%type  = xtype
    iodesc_list(iodnum)%ndims = ndims
    iodesc_list(iodnum)%dims  = 0
    iodesc_list(iodnum)%dims(1:ndims)   = dims(1:ndims)
    iodesc_list(iodnum)%dimids(1:ndims) = dimids(1:ndims)

  end subroutine ncd_getiodesc

end module ncdio_pio
