! netcdf4dnew.f    last update 4/7/98 (ccm)
! Updated so that can read in vars that don't have level as 3rd dim.
! The netcdf V3 version of netcdf4d.f
!
!     This contains three subroutines: ini4dnc initializes a netcdf file
! which will contain one idim x jdim x kdim x unlimited array and writes
! the first idim x jdim x kdim "record".  wr4dnc appends another 
! idim x jdim x kdim record onto an extant file.  rd4dnc reads a single
! idim x jdim array from the netcdf file, along with title and units.
! Upon writing an array, ndim >= idim*jdim*kdim
! Upon reading an array, ndim >= idim*jdim
!

      subroutine ini4dnc ( array, ndim, idim, jdim, kdim, &
                uarray, uidim, ujdim, ukdim, urdim, &
                validim, valjdim, valkdim, valrdim,  &
                aname, title, filen)

      include 'netcdf.inc'

      integer   idcdf, ierr, iomode, idlat, idlon, idlev
      integer   idtim, la, ntime
      integer   ivdims(4), istart(4), icnt(4)
      integer   idim, jdim, kdim
      integer   ndim, nn, j, i, k, ntry
      real      array(ndim), validim(idim), amiss
      real      valjdim(jdim), valkdim(kdim), valrdim
      real      addoff, scalefac
      character*(*) uarray, uidim, ujdim, ukdim
      character*(*) urdim, aname, title, filen

      data addoff, scalefac / 0., 1. /

      amiss = 9.e+20

!     Dimension error checking
!     ------------------------
      if ( idim .eq. 0 .or. jdim .eq. 0 .or. kdim .eq. 0 ) then
         print *, 'ERROR in ini4dnc: one of the dimensions is zero'
         print *, 'idim = ', idim
         print *, 'jdim = ', jdim
         print *, 'kdim = ', kdim
         stop 1
      else if ( idim .gt. 999999 .or. jdim .gt. 999999 &
           .or. kdim .gt. 999999 ) then
         print *, 'ERROR in ini4dnc: one of the dimensions is'
         print *, 'greater than 999999'
         print *, 'idim = ', idim
         print *, 'jdim = ', jdim
         print *, 'kdim = ', kdim
         stop 1
      end if

!     Data array error checking
!     -------------------------
      nn = 0
      do 5 j = 1, jdim
         do 6 i = 1, idim
            do 7 k = 1, kdim
               nn = nn + 1
               if ( abs(array(nn)) .gt. 1.e+35 ) then
                  !print *, 'WARNING in ini4dnc: ', aname
                  !print *, 'contains an absolute value greater than'
                  !print *, '1.e+29'
                  !print *, 'i, j, k, value = ', i, j, k, array(nn)
               end if
 7          continue
 6       continue
 5    continue

!     Initial starting point and length for each dimension
!     ----------------------------------------------------
      data istart / 1, 1, 1, 1 /

      icnt(1) = idim
      icnt(2) = jdim
      icnt(3) = kdim
      icnt(4) = 1


!     Create and open a new netcdf file (w/ clobber) in define mode
!     -------------------------------------------------------------
      ntry = 0
 10   ierr = NF_CREATE (filen, NF_CLOBBER, idcdf)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         ntry = ntry + 1
         if ( ntry .gt. 3 ) stop
         print *, filen
         print *, 'directory incorrect or file already exists'
         print *, 'Enter another complete file name to try again: '
         read (*,'(a80)') filen
         go to 10
      end if

!     Create a global title
!     ---------------------
      ierr = NF_PUT_ATT_TEXT (idcdf, NF_GLOBAL, 'title', 80, title)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 2
      end if

!     Create lat,lon, level, and time dimensions
!     ------------------------------------------
      ierr = NF_DEF_DIM (idcdf, 'lon', idim, idlon)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 3
      end if
      ierr = NF_DEF_DIM (idcdf, 'lat', jdim, idlat)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 3
      end if
      ierr = NF_DEF_DIM (idcdf, 'level', kdim, idlev)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 3
      end if
      ierr = NF_DEF_DIM (idcdf, 'time', NF_UNLIMITED, idtim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 3
      end if

!     Create floating point variables for dimension variables
!     -------------------------------------------------------
      ivdims(1) = idlon
      ierr = NF_DEF_VAR (idcdf, 'lon', NF_FLOAT, 1, ivdims, idvlon)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 4
      end if
      ivdims(1) = idlat
      ierr = NF_DEF_VAR (idcdf, 'lat', NF_FLOAT, 1, ivdims, idvlat)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 4
      end if
      ivdims(1) = idlev
      ierr = NF_DEF_VAR (idcdf, 'level', NF_FLOAT, 1, ivdims, idvlev)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 4
      end if
      ivdims(1) = idtim
      ierr = NF_DEF_VAR (idcdf, 'time', NF_FLOAT, 1, ivdims, idvtim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 4
      end if

!     Create a float pt variable of dimension (idim,jdim,kdim,ldim)
!     -------------------------------------------------------------
      la = index(aname,'$') - 1
      if (la .le. 1) la = len(aname)
      ivdims(1) = idlon
      ivdims(2) = idlat
      ivdims(3) = idlev
      ivdims(4) = idtim
      ierr = NF_DEF_VAR (idcdf, aname(1:la), NF_FLOAT, 4, ivdims, idvar)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 4
      end if

!     Create units attributes for all variables
!     -----------------------------------------
      ierr = NF_PUT_ATT_TEXT (idcdf, idvlon, 'units', 80, uidim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 5
      end if
      ierr = NF_PUT_ATT_TEXT (idcdf, idvlat, 'units', 80, ujdim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 5
      end if
      ierr = NF_PUT_ATT_TEXT (idcdf, idvlev, 'units', 80, ukdim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 5
      end if
      ierr = NF_PUT_ATT_TEXT (idcdf, idvtim, 'units', 80, urdim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 5
      end if
      ierr = NF_PUT_ATT_TEXT (idcdf, idvar, 'units', 80, uarray)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 5
      end if

!     Create attributes add_offset and scale_factor for each var
!     ----------------------------------------------------------
      ierr = NF_PUT_ATT_REAL (idcdf, idvlon, 'add_offset',NF_FLOAT, 1, addoff)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvlon, 'scale_factor',NF_FLOAT, 1, scalefac)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvlat, 'add_offset',NF_FLOAT, 1, addoff)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvlat, 'scale_factor',NF_FLOAT, 1, scalefac)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvlev, 'add_offset',NF_FLOAT, 1, addoff)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvlev, 'scale_factor',NF_FLOAT, 1, scalefac)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvtim, 'add_offset',NF_FLOAT, 1, addoff)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvtim, 'scale_factor',NF_FLOAT, 1, scalefac)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvar, 'add_offset',NF_FLOAT, 1, addoff)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if
      ierr = NF_PUT_ATT_REAL (idcdf, idvar, 'scale_factor',NF_FLOAT, 1, scalefac)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if

!     Create attribute missing_value for 4d var
!     -----------------------------------------
      ierr = NF_PUT_ATT_REAL (idcdf, idvar, 'missing_value',NF_FLOAT, 1, amiss)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 7
      end if

!     Leave define mode and enter data mode
!     -------------------------------------
      ierr = NF_ENDDEF (idcdf)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 8
      end if

!     Set nofill mode
!     ---------------
      ierr = NF_SET_FILL (idcdf, NF_NOFILL, iomode)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 8
      end if

!     Write values into dimension variables and data variable
!     -------------------------------------------------------
      ierr = NF_PUT_VAR_REAL (idcdf, idvlon, validim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 9
      end if
      ierr = NF_PUT_VAR_REAL (idcdf, idvlat, valjdim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 9
      end if
      ierr = NF_PUT_VAR_REAL (idcdf, idvlev, valkdim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 9
      end if
      ierr = NF_PUT_VAR1_REAL (idcdf, idvtim, istart, valrdim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 9
      end if

      ierr = NF_PUT_VARA_REAL (idcdf, idvar, istart, icnt, array)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 9
      end if

!     Close cdf file
!     --------------
      ierr = NF_CLOSE (idcdf)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 10
      end if

      return
      end

!...........................................................................


      subroutine wr4dnc ( array, ndim, idim, jdim, kdim, ntime, &
                         valrdim, aname, filen)

      include 'netcdf.inc'

      integer   idcdf, ierr, idvar
      integer   istart(4), icnt(4)
      integer   idim, jdim, kdim, i, j, k, nn, ndim
      real      array(ndim), valrdim
      character*(*) aname, filen


!     Time index error checking
!     -------------------------
      if (ntime .le. 0) then
         print *, 'ERROR in wr4dnc: time index ntime is'
         print *, 'less than or equal to zero'
         print *, 'ntime = ', ntime
         stop 1
      else if (ntime .gt. 999999) then
         print *, 'ERROR in wr4dnc: time index ntime is'
         print *, 'greater than 999999'
         print *, 'ntime = ', ntime
         stop 1
      end if

!     Data array error checking
!     -------------------------
      nn = 0
      do 5 j = 1, jdim
         do 6 i = 1, idim
            do 7 k = 1, kdim
               nn = nn + 1
               if ( abs(array(nn)) .gt. 1.e+35 ) then
                  print *, 'WARNING in wr4dnc: ', aname
                  print *, 'contains an absolute value greater than'
                  print *, '1.e+29'
                  print *, 'i, j, k, value = ', i, j, k, array(nn)
               end if
 7          continue
 6       continue
 5    continue

!     Initial starting point and length for each dimension
!     ----------------------------------------------------
      data istart / 1, 1, 1, 1 /

      icnt(1) = idim
      icnt(2) = jdim
      icnt(3) = kdim
      icnt(4) = 1

      istart(4) = ntime


!     Open file in write mode
!     -----------------------
      ierr = NF_OPEN (filen, NF_WRITE, idcdf)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 2
      end if

!     Get variable id's
!     -----------------
      ierr = NF_INQ_VARID (idcdf, 'time', idvtim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 3
      end if
      la = index(aname,'$') - 1
      if (la .le. 1) la = len(aname)
      ierr = NF_INQ_VARID (idcdf, aname(1:la), idvar)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 3
      end if

!     Write values into dimension variables and data variable
!     -------------------------------------------------------
      istart(1) = 1
      ierr = NF_PUT_VARA_REAL (idcdf, idvar, istart, icnt, array)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 4
      end if

      istart(1) = ntime
      ierr = NF_PUT_VAR1_REAL (idcdf, idvtim, istart, valrdim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 5
      end if

!     Close cdf file
!     --------------
      ierr = NF_CLOSE (idcdf)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         stop 6
      end if

      return
      end


!....................................................................

      subroutine rd4dnc ( array, ndim, idim, jdim, nlevel, ntime, &
                         aname, units, title, filen, ierr )

      implicit none
      include 'netcdf.inc'
      integer       :: ndim
      integer       :: idim
      integer       :: jdim
      integer       :: nlevel
      integer       :: ntime
      character*(*) :: aname
      character*(*) :: units
      character*(*) :: title
      character*(*) :: filen
      integer       :: ierr
      real          :: array(ndim)
      integer       :: istart(4)
      integer       :: icount(4)
      character*80  :: junk
      integer       :: idcdf
      integer       :: idtim
      integer       :: nrecords
      integer       :: idlev
      integer       :: nlevs
      integer       :: la
      integer       :: idvar
      integer       :: ndims
!     Open netcdf file in no write mode
!     ---------------------------------
      idcdf=10
      ierr = NF_OPEN (filen, NF_NOWRITE, idcdf)
!      WRITE(*,*)ierr,filen, NF_NOWRITE, idcdf
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         return
      end if


!     Time, level index error checking
!     --------------------------------
      idtim=20
      ierr = NF_INQ_DIMID (idcdf, 'time', idtim)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         print *, 'error in obtaining time id in rd4dnc'
         stop 1
      end if
      ierr = NF_INQ_DIMLEN (idcdf, idtim, nrecords)

      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         print *, 'error in obtaining time size in rd4dnc'
         stop 1
      end if
      ierr = NF_INQ_DIMID (idcdf, 'level', idlev)
      if ( ierr .ne. NF_NOERR ) then
!*         print *, NF_STRERROR(ierr)
!*         print *, 'error in obtaining level id in rd4dnc'
         nlevs = -999
      else
         ierr = NF_INQ_DIMLEN (idcdf, idlev, nlevs)
         if ( ierr .ne. NF_NOERR ) then
            print *, NF_STRERROR(ierr)
            print *, 'error in obtaining level size in rd4dnc'
            stop 1
         end if
         if (nlevel .gt. nlevs) then
            print *, 'ERROR in wr4dnc: nlevel is'
            print *, 'greater than maximum value in file'
            print *, 'nlevel = ', nlevel
            stop 1
         end if
      end if
      if (ntime .le. 0) then
         print *, 'ERROR in wr4dnc: ntime is'
         print *, 'less than or equal to zero'
         print *, 'ntime = ', ntime
         stop 1
      else if (ntime .gt. nrecords) then
         print *, 'ERROR in wr4dnc: ntime'
         print *, 'greater than maximum value in file'
         print *, 'ntime = ', ntime
         stop 1
      end if

!     Get id for the data array
!     -------------------------
      la = index(aname,'$') - 1
      if (la .le. 1) la = len(aname)
      ierr = NF_INQ_VARID (idcdf, aname(1:la), idvar)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         print *, 'error in obtaining data id in rd4dnc'
         stop 2
      end if

!     Check to see if it has 3 or 4 dims (has level?)
!     -----------------------------------------------
      ierr = NF_INQ_VARNDIMS(idcdf, idvar, ndims)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         print *, 'error in obtaining data ndims in rd4dnc'
         stop 2
      end if
      if (ndims .lt. 4) nlevs = -999

!     Get some attributes
!     -------------------
      ierr = NF_GET_ATT_TEXT (idcdf, idvar, 'units', units)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         print *, 'error in obtaining units in rd4dnc'
         stop 2
      end if
      ierr = NF_GET_ATT_TEXT (idcdf, NF_GLOBAL, 'title', title)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         print *, 'error in obtaining title in rd4dnc'
!        stop 2
      end if

!     Get the data array for level nlevel and time ntime
!     --------------------------------------------------
      if (nlevs .ne. -999) then
         istart(1) = 1
         istart(2) = 1
         istart(3) = nlevel
         istart(4) = ntime
         icount(1) = idim
         icount(2) = jdim
         icount(3) = 1
         icount(4) = 1
      else
         istart(1) = 1
         istart(2) = 1
         istart(3) = ntime
         icount(1) = idim
         icount(2) = jdim
         icount(3) = 1
      end if
      ierr = NF_GET_VARA_REAL (idcdf, idvar, istart, icount, array)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         print *, 'error in obtaining data in rd4dnc'
         stop 3
      end if

!     Close netcdf file
!     -----------------
      ierr = NF_CLOSE (idcdf)
      if ( ierr .ne. NF_NOERR ) then
         print *, NF_STRERROR(ierr)
         print *, 'error in closing file in rd4dnc'
         stop 4
      end if


      return
      end
