!
      module netcdf_module
      use proc,only: nlat,nlon,spval
      use hist,only: history,hdr
      implicit none
      include "netcdf.inc"
!
      integer,private ::
     |  id_time,                   ! dimension id for time (unlimited)
     |  idv_mtime                  ! mtime variable id
!
! For efficient i/o on netcdf history, the f4d array will be 
! dimensioned f4d(h%nlon-2,h%nlat,h%nzp,h%nflds). Each 4-d field 
! is read at the requested mtime by sub nc_rdfld, which is called 
! from a fields loop in nc_rdhist. Then, in the j-loop in getflds, 
! the data is transferred from f4d to flatslice for processing at 
! each latitude. The f4d array is allocated by nc_rdhist, and 
! deallocated after the j-loop in getflds.
!
      real,allocatable :: f4d(:,:,:,:)  ! (nlon,nlat,nlev,nflds)
      real,dimension(nlon-1,nlat) :: tlbc,ulbc,vlbc
!
      contains
!-----------------------------------------------------------------------
      subroutine nc_open(ncid,flnm,status,rdwr)
!
! Open or create a netcdf dataset flnm, returning netcdf file
!   id handle in ncid.
! If status=='NEW' or 'REPLACE', a new netcdf file is created,
!   otherwise (status=='OLD'), open existing dataset).
! If an existing file is opened and rdwr=='READ', the dataset is
!   opened read-only (NF_NOWRITE), otherwise (rdwr=='WRITE' or 
!   'APPEND') the dataset is opened read/write (NF_WRITE).
! If the open fails, print error msg and return ncid == 0.
!
! Args:
      integer,intent(out) :: ncid    ! netcdf file id
      character(len=*),intent(in) :: flnm,status,rdwr
!
! Local:
      integer :: istat
      character(len=120) :: char120
!
      ncid = 0
      if (status=='NEW'.or.status=='REPLACE') then
!
! Create new dataset (clobber old one if it pre-exists):
!
        istat = nf_create(flnm,NF_CLOBBER,ncid)
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_create for netcdf ',
     |      'file ',a,' status=',a)") trim(flnm),status
          call handle_ncerr(istat,char120)
          ncid = 0
        else
          write(6,"('Created netcdf file ',a,' ncid=',i8)") 
     |      trim(flnm),ncid
        endif
!
! Open pre-existing dataset:
! NF_WRITE opens with read-write, where write means append, change, etc
! NF_NOWRITE is read-only.
!
      elseif (status=='OLD') then
        if (rdwr=='WRITE'.or.rdwr=='APPEND') then
          istat = nf_open(flnm,NF_WRITE,ncid)
        elseif (rdwr=='READ') then
          istat = nf_open(flnm,NF_NOWRITE,ncid)
        else
          write(6,"('>>> nc_open: unrecognized rdwr=',a)") trim(rdwr)
          write(6,"('  Will open NF_WRITE (read/write).')")
          istat = nf_open(flnm,NF_WRITE,ncid)
        endif
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_open for netcdf ',
     |      'file ',a,' status=',a)") trim(flnm),status
          call handle_ncerr(istat,char120)
          ncid = 0
        else
          write(6,"('Opened existing netcdf file ',a,' ncid=',i8)") 
     |      trim(flnm),ncid
        endif
!
! Unknown status request:
      else
        write(6,"(/,'>>> nc_open: unrecognized status = ',a)")
     |    status
        ncid = 0
      endif
      end subroutine nc_open
!-------------------------------------------------------------------
      subroutine nc_close(ncid)
!
! Close a netcdf dataset file:
!
! Arg:
      integer,intent(in) :: ncid
!
! Local:
      integer :: istat
!
      istat = nf_close(ncid)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error return from nf_close')
      end subroutine nc_close
!-------------------------------------------------------------------
      subroutine nc_rdhist(ncid,mssfile,dskfile,mtime,itime,h,ier,
     |  iprint)
      use hist,only: print_hist
      use input,only: zkmbot_mgcm
!
! Args:
      integer,intent(in)  :: ncid,iprint
      integer,intent(inout)  :: mtime(3) ! may be changed if < 0 or > ntimes
      integer,intent(out) :: itime,ier
      character(len=*),intent(in) :: mssfile,dskfile
      type(history),intent(inout) :: h
!
! Local:
      integer,parameter :: mxdims=20 ! assume <= mxdims dimensions
      integer,parameter :: mxzp=100  ! assume nzp <= 100
      integer,parameter :: mxfields_read=200  ! max # model fields
      integer,parameter :: nsdtide=10, ndtide=2
      real :: sdtide(nsdtide),dtide(ndtide),dzp
      integer :: i,istat,ntimes,istart2(2),icount2(2),mtime_read(3),
     |  lenvol,len,itype,iddims(mxdims),iyear,iday,nflds_read,igpi,incep
      integer :: ndims,nvars,ngatts,idunlim,natts,ivar1(1),
     |  iscalar,ixf4d(mxfields_read),ixf3d(mxfields_read)
      character(len=80) :: varname,history_type
      real :: var1(1),fscalar
      character(len=8) :: fnames_read(mxfields_read)
      logical :: assoc
!
!     write(6,"('nc_rdhist: ncid=',i8,' dskfile=',a)") 
!    |  ncid,trim(dskfile)
      ier = 0
!
! Get id for mtime variable:
      istat = nf_inq_varid(ncid,"mtime",idv_mtime) 
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,'Error getting id of variable mtime')
        ier = 1
        return
      endif
!
! Get number of histories (times) on the file (ntimes):
      istat = nf_inq_unlimdim(ncid,id_time) ! id of unlimited record var
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting unlimited dimension id')
      istat = nf_inq_dimlen(ncid,id_time,ntimes) ! length of time var
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting length of time record variable')
!
! Search for requested model time. When found, retain time index itime:
      istart2(1) = 1
      icount2(1) = 3
      icount2(2) = 1
      itime = 0
      time_loop: do i=1,ntimes
        istart2(2) = i
        istat = nf_get_vara_int(ncid,idv_mtime,istart2,icount2,
     |    mtime_read)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'Error getting mtime values')
!
! If model day is < 0, use first history on the file: 
        if (mtime(1) < 0) then
          itime = 1
          mtime = mtime_read
          write(6,"('Requested day is < 0 -> will read first history ',
     |      'on the file: ',3i4)") mtime
          exit time_loop
        endif
        istat = nf_get_att_text(ncid,NF_GLOBAL,"model_version",
     |    h%version)
        h%ismtgcm = .false. 
        if (trim(h%version)=='MTGCM'.or.h%version(1:5)=='mtgcm')
     |    h%ismtgcm = .true.
!
! If model day > 365, look for day-365 (crossing year-boundary):
! 11/25/09 btf: allow mtgcm model day to be > 365
! 4/15/14 btf: allow non-mtgcm model day to be 366, for leap-year.
!
        if (.not.h%ismtgcm) then
          if (mtime(1) > 366) mtime(1) = mtime(1)-366
        endif
!
        if (iprint > 0) 
     |    write(6,"('Searching for history ',3i4,' Found ',3i4,' i=',i4,
     |      ' ntimes=',i4)") mtime(1:3),mtime_read,i,ntimes
        if (all(mtime_read==mtime)) then
          write(6,"('Found ',3i4,' on file ',a,' n=',i3)") 
     |      mtime,trim(dskfile),i
          itime = i
          exit time_loop
        endif
      enddo time_loop
      if (itime==0) then
        write(6,"(/,'>>> nc_rdhist: mtime ',3i4,' NOT found on ',
     |    'file ',a)") mtime,trim(dskfile)
        ier = 1
        return
      endif
!
! History was found: define history structure from netcdf vars:
! (history h was been initialized by sub gethist)
!
      h%ncid = ncid
      h%lu = 0
      h%mtime(:) = mtime(:)
      h%itime = itime
!     h%histfile = trim(dskfile) ! was set to vols(iv) in gethist
!     istat = nf_get_att_text(ncid,NF_GLOBAL,"mss_path",h%mssvol)
!     if (istat /= NF_NOERR) call handle_ncerr(istat,
!    |  'Error getting global file attribute mss_path')
      istat = nf_get_att_text(ncid,NF_GLOBAL,"model_version",
     |  h%version)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"model_name",h%model_name)
      if (trim(h%model_name)=='time-gcm') then
        h%istimes = .true.
      else
        h%istimes = .false.
      endif
      h%isdyn = .true. ! this is an assumption
      h%ismtgcm = .false. 
      if (trim(h%version)=='MTGCM'.or.h%version(1:5)=='mtgcm') then
        h%istimes = .false.
        h%isdyn = .false.
        h%isjtgcm = .false.
        h%isvtgcm = .false.
        h%ismtgcm = .true.
      endif 
      h%isjtgcm = .false. 
      if (h%version(1:5)=='JTGCM'.or.h%version(1:5)=='jtgcm') then
        h%istimes = .false.
        h%isdyn = .false.
        h%ismtgcm = .false.
        h%isvtgcm = .false.
        h%isjtgcm = .true.
        write(6,"('nc_rdhist: this is a jupiter jtgcm history.')")
      endif 
      h%isvtgcm = .false. 
      if (h%version(1:5)=='VTGCM'.or.h%version(1:5)=='vtgcm') then
        h%istimes = .false.
        h%isdyn = .false.
        h%ismtgcm = .false.
        h%isjtgcm = .false.
        h%isvtgcm = .true.
        h%version = 'VTGCM' 
        write(6,"('nc_rdhist: this is a venus vtgcm history.')")
      endif 
      istat = nf_get_att_text(ncid,NF_GLOBAL,"history_type",
     |  history_type)
      if (history_type(1:3)=='pri') then
        h%issech = .false.
      elseif (history_type(1:3)=='sec') then
        h%issech = .true.
      else
        write(6,"('>>> nc_rdhist: unrecognized history_type=',a)")
     |    trim(history_type)
      endif
!
!  Need to add coupled_ccm to model netcdf output, and use that
!  to define h%iscpl here.
!
! Get number of dims, vars, atts, and id of unlimited dimension:
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
      do i=1,ndims
        istat = nf_inq_dim(ncid,i,varname,len)
        if (trim(varname)=='lon') then
! nlon includes wrap point, history dimension does not,
! (i.e., nlon==73, history lon len==72)
          if (len /= nlon-1) then 
            write(6,"(/,'>>> nc_rdhist: length of lon read from',
     |        ' history does not match parameter nlon.')")
            write(6,"('  len=',i4,' nlon=',i4)") len,nlon
          endif
          h%nlon = len+2 ! 74
        endif
        if (trim(varname)=='lat') then
!         write(6,"('nc_rdhist: h%isvtgcm=',l1,' len=',i3,' nlat=',i3)")
!    |      h%isvtgcm,len,nlat
 
!
! 9/5/06 btf: vtgcm nlat now == full global 36
          if (len /= nlat) then 
            write(6,"(/,'>>> nc_rdhist: length of lat read from',
     |        ' history does not match parameter nlat.')")
            write(6,"('  len=',i4,' nlat=',i4)") len,nlat
          endif

!         if (len /= nlat .and. (h%isvtgcm.and.len/=nlat/2)) then
!           write(6,"(/,'>>> nc_rdhist: length of lat read from',
!    |        ' history does not match parameter nlat.')")
!           write(6,"('  len=',i4,' nlat=',i4)") len,nlat
!         endif

          h%nlat = len
        endif
        if (trim(varname)=='lev') then
          h%nzp = len ! kmaxp1
        endif
      enddo
!
! Allocate history vertical coord arrays:
!
      allocate(h%lev(h%nzp),stat=ier)
      if (ier /= 0) write(6,"('>>> nc_rdhist: error allocating ',
     |  'h%lev: h%nzp=',i3)") h%nzp
      allocate(h%ilev(h%nzp),stat=ier)
      if (ier /= 0) write(6,"('>>> nc_rdhist: error allocating ',
     |  'h%ilev: h%nzp=',i3)") h%nzp
!
      iyear = -1 ; iday = -1
!
! Get number of f4d fields on the history:
      h%nflds = 0
      ixf4d(:) = 0
      ixf3d(:) = 0
      do i=1,nvars
        if (isf4d(ncid,i,h%nzp)) then
          h%nflds = h%nflds+1
          ixf4d(i) = h%nflds
        elseif (isf3d(ncid,i)) then
          h%nflds = h%nflds+1
          ixf3d(i) = h%nflds
        endif
      enddo

!     write(6,"('nc_rdhist: found ',i3,' total (3d+4d) fields on ',
!    |  'input file')") h%nflds
!
! Allocate main 4-d fields array:
      if (allocated(f4d)) deallocate(f4d)
      allocate(f4d(h%nlon-2,h%nlat,h%nzp,h%nflds),stat=ier)
      if (ier /= 0) then
        call allocerr(ier,
     +    'error allocating space for f4d in nc_rdhist')
      else
!       write(6,"('Allocated f4d: h%nlon-2=',i3,' h%nlat=',
!    |    i3,' h%nzp=',i3,' h%nflds=',i3,' total = ',i8)")
!    |    h%nlon-2,h%nlat,h%nzp,h%nflds, (h%nlon-2)*h%nlat*h%nzp*h%nflds
      endif

      h%isnew = .false.
      h%lbc = spval
!
! Define history vars:
      varloop: do i=1,nvars
        istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts)
        select case(trim(varname))
          case('lon')
!           istat = nf_get_var_real(ncid,i,glon) ! do not need these now
          case('lat') 
!           istat = nf_get_var_real(ncid,i,glat) ! do not need these now
          case('lev')
            istat = nf_get_var_real(ncid,i,h%lev)
            if (istat /= NF_NOERR) 
     |        call handle_ncerr(istat,"reading lev coord var")
            h%zpt = h%lev(h%nzp)
            h%zpb = h%lev(1)
!           write(6,"(/,'nc_rdhist: h%zpb,zpt=',2f9.3,' h%lev=',/,
!    |        (8f9.3))") h%zpb,h%zpt,h%lev
!
! Interface levels are on new histories only (after about Nov, 2005 for tiegcm)
          case('ilev')
            istat = nf_get_var_real(ncid,i,h%ilev)
            if (istat /= NF_NOERR) 
     |        call handle_ncerr(istat,"reading ilev coord var")
            h%zpit = h%ilev(h%nzp)
            h%zpib = h%ilev(1)
!           write(6,"('nc_rdhist: h%zpib,zpit=',2f9.3,' h%ilev=',/,
!    |        (8f9.3))") h%zpib,h%zpit,h%ilev

!         case('calday')  ; istat = nf_get_var1_int(ncid,i,itime,iday)
 
          case('year')
            istat = nf_get_var1_int(ncid,i,itime,iyear)
            h%year = iyear
!           write(6,"('nc_rdhist: h%year=',i3)") h%year
          case('day')
            istat = nf_get_var1_int(ncid,i,itime,iday)
            h%day = iday
!           write(6,"('nc_rdhist: h%day=',i3)") h%day
          case('ut')   
            istat = nf_get_var1_real(ncid,i,itime,fscalar)
            h%ut = fscalar
!           write(6,"('nc_rdhist: h%ut=',f8.2)") h%ut
          case('timestep')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%step = iscalar
!           write(6,"('nc_rdhist: h%step=',i5)") h%step
          case('nflds') 
            istat = nf_get_var_int(ncid,i,ivar1)
!           h%nflds = ivar1(1) ! was already calculated using isf4d above
          case('dtide')
            istart2(1) = 1
            istart2(2) = itime
            icount2(1) = ndtide
            icount2(2) = 1
            istat = nf_get_vara_real(ncid,i,istart2,icount2,dtide)
            h%dtide(:) = dtide(:)
            hdr%dtide(:) = dtide(:)
!           write(6,"('nc_rdhist: h%dtide=',2e12.4)") h%dtide
          case('sdtide')
            istart2(1) = 1
            istart2(2) = itime
            icount2(1) = nsdtide
            icount2(2) = 1
            istat = nf_get_vara_real(ncid,i,istart2,icount2,sdtide)
            h%sdtide(:) = sdtide(:)
            hdr%sdtide(:) = sdtide(:)
!           write(6,"('nc_rdhist: h%sdtide=',/,(5e12.4))") h%sdtide
          case('f107d')
            istat = nf_get_var1_real(ncid,i,itime,var1)
            h%f107d = var1(1)
            hdr%f107d = var1(1)
!           write(6,"('nc_rdhist: h%f107d=',e12.4)") h%f107d
          case('f107a')
            istat = nf_get_var1_real(ncid,i,itime,var1)
            h%f107a = var1(1)
            hdr%f107a = var1(1)
!           write(6,"('nc_rdhist: h%f107a=',e12.4)") h%f107a
          case('hpower')
            istat = nf_get_var1_real(ncid,i,itime,var1)
            h%hpower = var1(1)
            hdr%hp = var1(1)
!           write(6,"('nc_rdhist: h%hpower=',e12.4)") h%hpower
          case('ctpoten')
            istat = nf_get_var1_real(ncid,i,itime,var1)
            h%ctpoten = var1(1)
            hdr%cp = var1(1)
!           write(6,"('nc_rdhist: h%ctpoten=',e12.4)") h%ctpoten
          case('byimf')
            istat = nf_get_var1_real(ncid,i,itime,var1)
            h%byimf = var1(1)
            hdr%byimf = var1(1)
!           write(6,"('nc_rdhist: h%byimf=',e12.4)") h%byimf
          case('gpi')
            istat = nf_get_var1_int(ncid,i,itime,igpi)
            h%gpi = igpi
!           write(6,"('nc_rdhist: h%gpi=',i3)") h%gpi
          case('ncep')
            istat = nf_get_var1_int(ncid,i,itime,incep)
            h%ncep = incep
!           write(6,"('nc_rdhist: h%ncep=',i3)") h%ncep
          case('gswmdi')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%gswmdi = iscalar
!           write(6,"('nc_rdhist: h%gswmdi=',i3)") h%gswmdi
          case('gswmsdi')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%gswmsdi = iscalar
!           write(6,"('nc_rdhist: h%gswmsdi=',i3)") h%gswmsdi
          case('gswmnmidi')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%gswmnmidi = iscalar
!           write(6,"('nc_rdhist: h%gswmnmidi=',i3)") h%gswmnmidi
          case('gswmnmisdi')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%gswmnmisdi = iscalar
!           write(6,"('nc_rdhist: h%gswmnmisdi=',i3)") h%gswmnmisdi
! tiegcm1.8 and later:
          case('gswmnmdi')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%gswmnmidi = iscalar
!           write(6,"('nc_rdhist: h%gswmnmidi=',i3,
!    |        ' (read as gswmnmdi)')") h%gswmnmidi
          case('gswmnmsdi')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%gswmnmisdi = iscalar
!           write(6,"('nc_rdhist: h%gswmnmisdi=',i3,
!    |        ' (read as gswmnmsdi)')") h%gswmnmisdi
!
! Lower boundary heights for mtgcm, from mgcm:
          case('zkmbot_mgcm')
            istart2(1) = 1
            istart2(2) = itime
            icount2(1) = h%nlat
            icount2(2) = 1
            istat = nf_get_vara_real(ncid,i,istart2,icount2,zkmbot_mgcm)
            if (istat /= NF_NOERR) 
     |        call handle_ncerr(istat,"reading zkmbot_mgcm")
            if (iprint > 0)
     |        write(6,"('nc_rdhist: zkmbot_mgcm=',/,(6e12.4))") 
     |          zkmbot_mgcm
          case('TLBC')
            call nc_rdlbc(ncid,i,itime,h,tlbc)
          case('ULBC')
            call nc_rdlbc(ncid,i,itime,h,ulbc)
          case('VLBC')
            call nc_rdlbc(ncid,i,itime,h,vlbc)
          case('LBC')
            istat = nf_get_var1_real(ncid,i,itime,var1)
            h%lbc = var1(1)
            h%isnew = .true.
          case default
!             write(6,"('Note nc_rdhist: unused variable: ',a)") 
!    |          trim(varname)
        end select
      enddo varloop ! i=1,nvars
      if (iyear > -1 .and. iday > -1) then
        h%iyd = (iyear-(iyear/100*100))*1000 + iday
!       write(6,"('nc_rdhist: iyear=',i5,' iday=',i3,' h%iyd=',i8)")
!    |    iyear,iday,h%iyd
      endif
!
! Old history has only lev array, which has interface values,
! therefore h%zpb,zpt were read as interface values. Also,
! all fields from old history were converted to interfaces.
! 
      if (h%zpib==0..and.h%zpit==0.) then
        dzp = (h%zpt-h%zpb)/float(h%nzp-1)
        h%zpib = h%zpb
        h%zpit = h%zpt
        h%ilev = h%lev
      endif
!
! Read 3d and 4d fields at current time:
!
      nflds_read = 0 ! number of 3d fields read from the history
      do i=1,nvars
        istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts)
!
! Read 4d field (3d+time):
        if (ixf4d(i) > 0) then
          call nc_rdfld(ncid,varname,itime,h,h%nzp,h%nlat,nflds_read,
     |      ier)
          if (ier > 0) then
            write(6,"('>>> WARNING nc_rdhist: error reading f4d ',
     |        'field ',a)") trim(varname)
          else
            fnames_read(nflds_read) = varname(1:8)
          endif
!
! Read 3d field (2d lat/lon + time):
! These will be defined in f4d at nflds_read, but levels data will be
! redundant (data will be identical at all levels).
!
        elseif (ixf3d(i) > 0) then
          call nc_rdfld3d(ncid,varname,itime,h,nflds_read,ier)
          if (ier > 0) then
            write(6,"('>>> WARNING nc_rdhist: error reading f3d ',
     |        'field ',a)") trim(varname)
          else
            fnames_read(nflds_read) = varname(1:8)
          endif
        endif
      enddo
!
! Allocate space for and define field names: 
!
      allocate(h%fnames(h%nflds),stat=ier)
      if (ier /= 0) call allocerr(ier,
     +  'error allocating space for field names in nc_rdhist')
      do i=1,h%nflds
        h%fnames(i) = fnames_read(i)
      enddo
!
      if (iprint > 0) then
        write(6,"(' ')")
        call print_hist(h)
      endif

      end subroutine nc_rdhist
!-------------------------------------------------------------------
      subroutine nc_rdfld(ncid,name,itime,h,nzp,nlatrd,ixfld,ier)
      use fields,only: flds ! field structures
!
! Read a 3-d field from netcdf history on ncid at time itime.
! These are 4-d fields on the history (because they include time).
! If there are not 4 dimensions, or the 1st 3 dimensions do not
! correspond to nlon,nlat,nlev, then return ier > 0.
!
! Args:
      integer,intent(in) :: ncid,itime,nzp,nlatrd
      integer,intent(inout) :: ixfld
      integer,intent(out) :: ier
      character(len=*),intent(in) :: name
      type(history),intent(in) :: h
!
! Local:
      integer :: istat,idvar,i,k,j,start_4d(4),count_4d(4)
      character(len=120) :: char120
      real :: f3d(nlon-1,nlatrd,nzp),fmin,fmax
      integer :: ixf,nfields,nnans
      character(len=56) :: long_name,vdimname
      character(len=16) :: units
      character(len=80) :: char80
      integer,external :: ixfindc
      integer :: iddims(4)
!
! If check_nan is set, check for presence of NaN's in data when reading.
! Sub check_nan currently works only on IBM:
! 11/29/05 btf: This should be checked (check_nan may also work on 
!               other platforms)
!
!#ifdef AIX
!      logical :: check_nan = .true.
!#else
!      logical :: check_nan = .false.
!#endif
       logical :: check_nan = .true.
!
!     write(6,"('Enter nc_rdfld: name=',a,' itime=',i3,' nzp=',i3,
!    |  ' ixfld=',i3)") trim(name),itime,nzp,ixfld
!
      ier = 0
!
! Get var id from name:
      istat = nf_inq_varid(ncid,name,idvar)
      if (istat /= NF_NOERR) then
        write(char120,"('nc_rdfld: Error return from nf_inq_varid',
     +    ' for field ',a,' itime=',i2)") trim(name),itime
        call handle_ncerr(istat,char120)
        ier = 1
        return
      else
!       write(6,"('nc_rdfld: field ',a,': idvar=',i3)") trim(name),
!    |    idvar
      endif
!
! Determine name of vertical dimension for this field ("lev" or "ilev")
! (at this point we can safely assume that lev or ilev is the 3rd dim)
!
      istat = nf_inq_vardimid(ncid,idvar,iddims)
      istat = nf_inq_dimname(ncid,iddims(3),vdimname)
      if (vdimname /= 'lev'.and.vdimname /= 'ilev') then
        write(6,"('>>> nc_rdfld: unknown lev dimname for field ',
     |    a,': ',a)") name,vdimname
      endif
!
! Read the data into f3d:
      start_4d(1:3) = 1
      start_4d(4) = itime
      count_4d(1) = nlon-1
      count_4d(2) = nlatrd
      count_4d(3) = nzp
      count_4d(4) = 1 
!
      istat = nf_get_vara_real(ncid,idvar,start_4d,count_4d,f3d)
      if (istat /= NF_NOERR) then

        write(6,"('Error return from nf_get_vara_real',
     |    ' for field var ',a,' itime=',i2)") trim(name),itime
        write(6,"('itime=',i4,' nlon-1=',i4,' nlatrd=',i4,' nzp=',
     |    i4,' ixfld=',i4)") itime,nlon-1,nlatrd,nzp,ixfld
        write(6,"('f3d min,max=',2e12.4)") minval(f3d),maxval(f3d)
        if (any(f3d==NF_FILL_DOUBLE)) then
          write(6,"('Found NF_FILL_DOUBLE=',e12.4,' in f3d.')")
        endif
        call check_nans(f4d(:,:,:,ixfld),nlon-1,nlatrd,nzp,
     |    name(1:8),nnans,0,0.,1,0)

!     real :: f3d(nlon-1,nlatrd,nzp),fmin,fmax
!       do k=1,nzp 
!         do j=1,nlatrd 
!           write(6,"('nc_rdfld: itime=',i4,' k=',i4,' j=',i4,
!    |        ' f3d(1:nlon-1,j,k)=',/,(6e12.4))") itime,k,j,
!    |        f3d(1:nlon-1,j,k)
!         enddo
!       enddo

        write(char120,"('Error return from nf_get_vara_real',
     |    ' for field var ',a,' itime=',i2)") trim(name),itime
        call handle_ncerr(istat,char120)
        ier = 1
!       return
      endif ! error reading
      ixfld = ixfld+1

! real :: f3d(nlon-1,nlatrd,nzp)
!     write(6,"('nc_rdfld: Read 3d field ',a,' ixfld=',i3,' ixf=',i3,
!    |  ' nlon-1=',i4,' nlatrd=',i4,' nzp=',i4,' min,max=',2e12.4)") 
!    |  name(1:8),ixfld,ixf,nlon-1,nlatrd,nzp,minval(f3d),maxval(f3d)

!
! NF_FILL_DOUBLE is netcdf's pre-fill value, i.e., if the variable
! was not written to the netcdf file by the model (e.g., the first
! history on secondary files), then netcdf pads the var with this 
! prefill value.  The prefill value is printed as "_" by ncdump. 
! Here, replace any prefill values read with spval, so it will not 
! be plotted: (the default NF_FILL_DOUBLE value is 0.9969E+37)
!
      where(f3d==NF_FILL_DOUBLE) f3d = spval
!
! Transfer to 4d global array in module data:
! allocate(f4d(h%nlon-2,h%nlat,h%nzp,h%nflds),stat=ier)
! real f3d(nlon-1,nlat,nzp),fmin,fmax
!
      f4d(:,:,:,ixfld) = f3d(:,:,:)
!
! Check for NaNs (sub check_nans is in util.F):
! subroutine check_nans(f,id1,id2,id3,name,ispval,spval,iprint,ifatal)
!   If ispval > 0 -> replace any INF or NaNs with spval
!   If iprint==1  -> print warnings only if INF or NaNs are found
!   If iprint==2  -> always print number of INF and NaNs found
!   If ifatal > 0 -> stop program when first INF or NaNs are found
!
! Do not replace NaNs here with spval, just print warnings -- plot 
! routines will print -NaNQ for min,max on the plots.
!
      if (check_nan) 
     |  call check_nans(f4d(:,:,:,ixfld),nlon-1,nlatrd,nzp,
     |    name(1:8),nnans,0,0.,1,0)
!
!     call fminmax(f4d(1,1,1,ixfld),(nlon-1)*nlatrd*nzp,fmin,fmax)
!     write(6,"('Read field ',a,' ixfld=',i3,' min,max=',
!    |  2e12.4)") name(1:8),ixfld,fmin,fmax
!
!     if (trim(name)=='Z') then
!     if (trim(name)=='TN') then
!       do j=1,nlatrd 
!         do k=1,nzp 
!           write(6,"('Read field ',a,' itime=',i4,' ixfld=',i3,
!    |        ' k=',i4,' j=',i4,' f4d(1:nlon-1,j,k,ixfld)=',/,
!    |        (6e12.4))") name(1:8),itime,ixfld,k,j,
!    |        f4d(1:nlon-1,j,k,ixfld)
!         enddo
!       enddo
!     endif
!
! If field is unknown, read long_name and units, if present:
!   (added by btf, 6/23/04)
!
      nfields = size(flds)
      ixf = ixfindc(flds%fname8,nfields,name)
      if (ixf > 0) then
        if (.not.flds(ixf)%known) then ! field is unknown
!
! Get long_name (istat==-43 means attribute was not found):
          long_name = ' '
          istat = nf_get_att_text(ncid,idvar,'long_name',long_name)
          if (istat /= NF_NOERR.and.istat /= -43) then
            write(char120,"('Error return from nf_get_att_text',
     |        ' for long_name of field ',a)") trim(name)
            call handle_ncerr(istat,char120)
          elseif (istat /= -43) then ! got long_name
            flds(ixf)%fname56 = long_name
!           write(6,"('nc_rdfld: long_name for unknown field ',a,
!    |        ': ',a)") trim(name),trim(long_name)
          endif
!
! Get units (istat==-43 means attribute was not found):
          units = ' '
          istat = nf_get_att_text(ncid,idvar,'units',units)
          if (istat /= NF_NOERR.and.istat /= -43) then
            write(char120,"('Error return from nf_get_att_text',
     |        ' for units of field ',a)") trim(name)
            call handle_ncerr(istat,char120)
          elseif (istat /= -43) then ! got units
            flds(ixf)%units = units
!           write(6,"('nc_rdfld: units for unknown field ',a,': ',a)")
!    |        trim(name),trim(units)
          endif
        endif   ! unknown field
!
! 3/13/06: Set vertical coordinate components according to lev vs ilev
! (vertical dimension name was determined with nf_inq_dimname above)
! This is for fields on the history only (known and unknown). Vertical 
!   components of derived fields will be set in sub setflev (getflds.F)
!
! However, default processing of "new" histories for T,U,V is conversion 
!   from midpoints to interfaces with T,U,VLBC as bottom interface level.
!   (see sub proclat in proclat.F)
! Note any field zptype can be overridden by the user with namelist read
!   grid_levels (see sub setflev in getflds.F)
!
!
        if (associated(flds(ixf)%lev)) deallocate(flds(ixf)%lev)
        allocate(flds(ixf)%lev(nzp),stat=istat)

        if (istat /= 0) then
          write(char80,"('nc_rdfld: error allocating lev for field ',
     |      a,': nzp=',i3)") name,nzp
          call allocerr(ier,char80)
        endif

        if (trim(name)=='TN'.or.trim(name)=='UN'.or.trim(name)=='VN') 
     |    then
          flds(ixf)%lev = h%ilev
          flds(ixf)%zptype = "INTERFACES"
        elseif (vdimname=='lev') then      ! midpoints
          flds(ixf)%lev = h%lev
          flds(ixf)%zptype = "MIDPOINTS"
        else                               ! interfaces
          flds(ixf)%lev = h%ilev
          flds(ixf)%zptype = "INTERFACES"
        endif
!
        flds(ixf)%nlev = nzp
        flds(ixf)%dlev = (flds(ixf)%lev(nzp)-flds(ixf)%lev(1))/
     |    float(nzp-1)
!
! For old histories, all fields are plotted at interfaces:
        if (.not.h%isnew) flds(ixf)%zptype = "INTERFACES"
!
!       write(6,"('nc_rdfld: ixf=',i3,' field ',a,' zptype=',a,
!    |    ' nlev=',i3,' levs=',f7.2,' to ',f7.2,' by ',
!    |    f7.2)") ixf,flds(ixf)%fname8,flds(ixf)%zptype,flds(ixf)%nlev,
!    |    flds(ixf)%lev(1),flds(ixf)%lev(flds(ixf)%nlev),
!    |    flds(ixf)%dlev
!
      else
!       write(6,"('nc_rdfld: did not find name ',a,
!    |    ' in flds%fname8 ixf=',i4)") trim(name),ixf
      endif
      end subroutine nc_rdfld
!-----------------------------------------------------------------------
      subroutine nc_rdfld3d(ncid,name,itime,h,ixfld,ier)
!
! Read a 3d field (lon,lat,time), defining f4d(lon,lat,lev,ixfld), where
! the lev dimension is redundant.
!
! Args:
      integer,intent(in) :: ncid,itime
      integer,intent(inout) :: ixfld
      character(len=*),intent(in) :: name
      type(history),intent(in) :: h
      integer,intent(out) :: ier
!
! Local:
      integer :: k,istat,idvar,start_3d(3),count_3d(3)
      character(len=120) :: char120
      real :: f2d(h%nlon-2,h%nlat)

      ier = 0
!
! Get var id from name:
      istat = nf_inq_varid(ncid,name,idvar)
      if (istat /= NF_NOERR) then
        ier = 1
        write(char120,"('>>> nc_rdfld: Error return from nf_inq_varid',
     +    ' for field ',a,' itime=',i2)") trim(name),itime
        call handle_ncerr(istat,char120)
        return
      else
!       write(6,"('nc_rdfld3d: field ',a,': idvar=',i3)") trim(name),
!    |    idvar
      endif
!
! Read the data into f2d:
      start_3d(1:2) = 1
      start_3d(3) = itime
      count_3d(1) = h%nlon-2
      count_3d(2) = h%nlat
      count_3d(3) = 1 
!
      istat = nf_get_vara_real(ncid,idvar,start_3d,count_3d,f2d)
      if (istat /= NF_NOERR) then
        write(char120,"('>>> nc_rdfld3d: Error reading f3d field',a)")
     |    trim(name)
        ier = 1
        return
      endif    
      ixfld = ixfld+1
!
! For this 3d lat-lon field, levels will be redundant:
      do k=1,h%nzp
        f4d(:,:,k,ixfld) = f2d(:,:)
      enddo
!
! Print min,max of first level (data is the same for all levels):
!     write(6,"('nc_rdfld3d: read field ',a,' ixfld=',i3,' min,max=',
!    |  2e12.4)") trim(name),ixfld,minval(f4d(:,:,1,ixfld)),
!    |                             maxval(f4d(:,:,1,ixfld))

      end subroutine nc_rdfld3d
!-----------------------------------------------------------------------
      logical function isf4d(ncid,idvar,nzp)
      implicit none
!
! Args:
      integer,intent(in) :: ncid,idvar,nzp
!
! Local:
      integer :: istat,itype,ndims,iddims(4),natts,idimsizes(4),i
      character(len=80) :: rdname,dimname
!
      isf4d = .false.
      istat = nf_inq_var(ncid,idvar,rdname,itype,ndims,iddims,natts) 
      if (ndims /= 4) return
      do i=1,ndims
        istat = nf_inq_dim(ncid,iddims(i),dimname,idimsizes(i))
      enddo
!     write(6,"('isf4d: idvar=',i3,' idimsizes=',4i4,' name=',a)") 
!    |  idvar,idimsizes,trim(rdname)
      if (idimsizes(1) /= nlon-1 .or. idimsizes(2)/=nlat .or.
     |    idimsizes(3) /= nzp) then
!       write(6,"('isf4d returning because dims for field ',a,
!    |    ' do not match: idimsizes=',4i4,' nlon-1=',i4,' nlat=',i4,
!    |    ' nzp=',i4)") trim(rdname),idimsizes,nlon-1,nlat,nzp
        return
      endif
      isf4d = .true.

!     write(6,"('isf4d found f4d field: idvar=',i3,' name=',a)")
!    |  idvar,trim(rdname)

      end function isf4d
!-----------------------------------------------------------------------
      logical function isf3d(ncid,idvar)
      implicit none
!
! Args:
      integer,intent(in) :: ncid,idvar
!
! Local:
      integer :: istat,itype,ndims,iddims(4),natts,idimsizes(4),i
      character(len=80) :: rdname,dimname
!
      isf3d = .false.
      istat = nf_inq_var(ncid,idvar,rdname,itype,ndims,iddims,natts) 
      if (ndims /= 3) return
      do i=1,ndims
        istat = nf_inq_dim(ncid,iddims(i),dimname,idimsizes(i))
      enddo

!     write(6,"('isf3d: idvar=',i3,' idimsizes=',4i4,' name=',a)") 
!    |  idvar,idimsizes,trim(rdname)

!     if (idimsizes(1) /= nlon-1 .or. idimsizes(2)/=nlat) then
!       write(6,"('isf3d returning because dims for field ',a,
!    |    ' do not match: idimsizes=',4i4,' nlon-1=',i4,' nlat=',i4
!    |    )") trim(rdname),idimsizes,nlon-1,nlat
!       return
!     endif
      isf3d = .true.

!     write(6,"('isf3d found f3d field: idvar=',i3,' name=',a)")
!    |  idvar,trim(rdname)

      end function isf3d
!-------------------------------------------------------------------
      subroutine nc_rdlbc(ncid,idvar,itime,h,flbc)
!
! Args:
      integer,intent(in) :: ncid,idvar,itime
      real,intent(out) :: flbc(:,:)
      type(history),intent(in) :: h
!
! Local:
      integer :: istat,itype,ndims,iddims(3),natts,idimsizes(3),i
      character(len=80) :: rdname,dimname
      character(len=120) :: char120
      integer :: start_3d(3),count_3d(3)
      real :: f2d(nlon-1,h%nlat)
!
      istat = nf_inq_var(ncid,idvar,rdname,itype,ndims,iddims,natts) 
      if (ndims /= 3) then
        write(6,"('>>> WARNING nc_rdlbc: field ',a,' is not 3-d:',
     |    ' ndims=',i3)") trim(rdname),ndims
        return
      endif
      do i=1,ndims
        istat = nf_inq_dim(ncid,iddims(i),dimname,idimsizes(i))
      enddo
      if (idimsizes(1) /= nlon-1 .or. idimsizes(2) /= h%nlat) 
     |  then
        write(6,"('>>> WARNING nc_rdlbc: field ',a,' dims do not',
     |    ' match: idimsizes=',3i4,' nlon-1=',i3,' h%nlat=',i3)")
     |    trim(rdname),idimsizes,nlon-1,h%nlat 
        return
      endif
!
      start_3d(1:2) = 1
      start_3d(3) = itime
      count_3d(1) = nlon-1
      count_3d(2) = h%nlat
      count_3d(3) = 1
      istat = nf_get_vara_real(ncid,idvar,start_3d,count_3d,f2d)
      flbc = f2d
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_get_vara_real',
     +    ' for lbc var ',a,' itime=',i2)") trim(rdname),itime
        call handle_ncerr(istat,char120)
      else
!       write(6,"('nc_rdlbc: read ',a,': min,max=',2e12.4)")
!    |    trim(rdname),minval(flbc),maxval(flbc)
      endif
      end subroutine nc_rdlbc
!-------------------------------------------------------------------
      end module netcdf_module