      subroutine readspec(datanm,what,itime_index,nyear)
c
c There is an example main program at the end of the file.
c
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c !                                                              !
c ! The names of the data files are hard-wired into the code.    !
c !                                                              !
c ! You'll need, the netCDF and UDUNITS libraries from UniData.  !
c !                                                              !
c ! netCDF can be found at:                                      !
c ! http://www.unidata.ucar.edu/packages/netcdf/index.html       !
c !                                                              !
c ! UDUNITS can be found at:                                     !
c ! http://www.unidata.ucar.edu/packages/udunits/index.html      !
c !                                                              !
c ! If you want to convert the data to a gaussian grid see:      !
c ! http://www.scd.ucar.edu/css/software/spherepack/             !
c ! and                                                          !
c ! http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/sphere/    !
c !                                                              !
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
c    Reads the CDC spectral coefficient netCDF files from 
c    the NCEP/NCAR Reanalysis. 
c    Spectral coefficients returned in complex array datanm
c    (dimension 2016x28).
c
c    input variables:
c    what is a character*4 field identifier.
c    allowable values are: 'vort','div ','shum','vair',
c    'lnps' and 'orog' (vorticity, divergence, specific
c    humidity, virtual temperature, natural log of surface
c    pressure (in centibars) and orographic height (meters)).
c
c    NOTE:  if what = 'lnps' or 'orog', values are returned
c    in (datanm(nm,1),nm=1,2016). Other 27 levels of
c    this array is unchanged.
c
c    itime_index is an integer time index denoting the 
c    record to read from the netCDF file.     
c
c    nyear is an integer denoting the year of data desired
c    (1984, 1989 etc)
c
      parameter (nlon = 192)
      parameter (nlat = 94)
      parameter (mwaves = 63)
      parameter (nmdim = (mwaves+1)*mwaves/2)
      parameter (levels = 28)
c
      complex datanm(nmdim,levels)
      character*4 what
c
      if (what .eq. 'vort' .or. what .eq. 'div ' .or.
     * what .eq. 'shum' .or. what .eq. 'vair') then
	call specread(datanm,what,itime_index,nyear)
      else if (what .eq. 'lnps') then
	call specreadps(datanm(1,1),itime_index,nyear)
      else if (what .eq. 'orog') then 
	call specreadzs(datanm(1,1))
      else
        write(6,*) 'invalid field identifier!' 
        write(6,*) 'allowed are:'
        write(6,*) 'vort,div ,shum,vair,lnps,orog'
        stop
      end if
c
      return
      end
c
	subroutine specread(datanm,what,time_index,nyear)
c
c==> reads spectral sigma files for what = 'vort', 'div',
c    'shum' and 'vair'.
c    Input: what (an character*4 denoting the name of the variable),
c    time_index (an integer denoting the output time
c    to read from the file) and nyear (an integer denoting
c    the year).
c    datanm is the output complex 2016x28 array of spectral coefficients.
c 
        INCLUDE '/usr/local/include/netcdf.inc'
c     
	parameter (mwaves = 63) 
	parameter (nmdim = (mwaves+1)*mwaves/2)
	parameter (nmdim2 = 2*nmdim)
	parameter (levels = 28)
	parameter (levels2 = 56)
c
	complex datanm(nmdim, levels)
c
        character*4 what
        character*85 ncfilename
c
	real mean(levels2), add_offset(levels), scale_factor(levels)
	integer*2 cscaled(2*nmdim)
c
	integer time_index,rcode
C
	integer ncid, specid, meanid, add_offsetid, scale_factorid
c
	integer cstart(3), ccount(3)

	data cstart /1,1,1/
	data ccount /nmdim2, 1, 1/

	integer mstart(2), mcount(2)
	data mstart /1,1/
	data mcount /levels2,1/
c      
	integer astart(2), acount(2)
	data astart /1,1/
	data acount /levels,1/


	   mstart(2) = time_index
	   astart(2) = time_index
	   cstart(3) = time_index

        len_what = index(what,' ') - 1
        if (len_what .eq. -1) len_what = 4
        if ( len_what .eq. 3 ) then
          write(ncfilename, 10) what(1:len_what), nyear
 10   format(
     *'/Datasets.large/archive0/Datasets.public/',
     *'ncep.reanalysis/spectral/',
     *   a3,'.spec.',i4,'.nc')
        else if ( len_what .eq. 4 ) then
          write(ncfilename, 20) what(1:len_what), nyear
 20   format('/Datasets.large/archive0/Datasets.public/',
     *'ncep.reanalysis/spectral/',
     *   a4,'.spec.',i4,'.nc')
        end if

        if (ncfilename(1:15) .eq. '/Datasets.large' ) then

           print*, "Change the path to the data files!"
           stop 99

        endif


	   ncid = ncopn(ncfilename, NCREAD, RCODE)
	   specid = ncvid(NCID,what,RCODE)
           meanid = ncvid(ncid, 'mean', RCODE)
           add_offsetid = ncvid(ncid, 'add_offset', RCODE)
           scale_factorid = ncvid(ncid, 'scale_factor', RCODE)
	   call ncvgt(ncid, meanid, mstart, 
     $                mcount, mean, RCODE)
	   call ncvgt(ncid, add_offsetid, astart, 
     $                acount, add_offset, RCODE)
	   call ncvgt(ncid, scale_factorid, astart, 
     $                acount, scale_factor, RCODE)

	   l = 1
	     do ii = 1, levels
	        cstart(2) = ii
		call ncvgt(ncid, specid, cstart, 
     $                     ccount, cscaled, RCODE)
		k = 3
		do i = 1, nmdim
			if( i .eq. 1 ) then
			   datanm(i,ii) = cmplx(mean(l), mean(l+1))
			   l = l + 2
			else
			   unpackr = scale_factor(ii)*
     $			             float(cscaled(k)) + add_offset(ii)
			   unpackc = scale_factor(ii)*
     $			             float(cscaled(k+1)) + add_offset(ii)
			   datanm(i,ii) = cmplx(unpackr, 
     $				             unpackc )
			   k = k + 2
			endif
		enddo
	     enddo


	call ncclos(ncid, RCODE)
        return
	end
c
	subroutine specreadps(datanm,time_index,nyear)
c
c==> reads spectral sigma files for natural log of sfce pressure.
c    Input: time_index (an integer denoting the output time
c    to read from the file) and nyear (an integer denoting
c    the year).
c    datanm is the output complex 2016 element array of spectral coefficients.
c 

        INCLUDE '/usr/local/include/netcdf.inc'
c     
	parameter (mwaves = 63) 
	parameter (nmdim = (mwaves+1)*mwaves/2)
	parameter (nmdim2 = 2*nmdim)
c
	complex datanm(nmdim)
c
        character*4 what
        character*95 ncfilename
c
        real mean(2), add_offset, scale_factor
	integer*2 cscaled(2*nmdim)
c
	integer time_index,rcode
C
	integer ncid, specid, meanid, add_offsetid, scale_factorid
c
      integer cstart(2), ccount(2)

      data cstart /1,1/
      data ccount /nmdim2, 1/

      integer mstart(2), mcount(2)
      data mstart /1,1/
      data mcount /2,1/

      integer astart
      data astart /1/


	   mstart(2) = time_index
	   astart    = time_index
	   cstart(2) = time_index

          what = 'pres'
          write(ncfilename, 10) nyear
 10   format('/Datasets.large/archive0/Datasets.public',
     * '/ncep.reanalysis/spectral/',
     * 'pres.nlog.sfc.spec.',
     * i4,'.nc')

        if (ncfilename(1:15) .eq. '/Datasets.large' ) then

           print*, "Change the path to the data files!"
           stop 99

        endif

	   ncid = ncopn(ncfilename, NCREAD, RCODE)
	   specid = ncvid(NCID,what,RCODE)
           meanid = ncvid(ncid, 'mean', RCODE)
           add_offsetid = ncvid(ncid, 'add_offset', RCODE)
           scale_factorid = ncvid(ncid, 'scale_factor', RCODE)
	   call ncvgt(ncid, meanid, mstart, 
     $                mcount, mean, RCODE)
	   call ncvgt1(ncid, add_offsetid, astart, 
     $                 add_offset, RCODE)
	   call ncvgt1(ncid, scale_factorid, astart, 
     $                 scale_factor, RCODE)

		call ncvgt(ncid, specid, cstart, 
     $                     ccount, cscaled, RCODE)
		k = 3
		do i = 1, nmdim
			if( i .eq. 1 ) then
			   datanm(i) = cmplx(mean(1), mean(2))
			else
			   unpackr = scale_factor*
     $			             float(cscaled(k)) + add_offset
			   unpackc = scale_factor*
     $			             float(cscaled(k+1)) + add_offset
			   datanm(i) = cmplx(unpackr, 
     $				             unpackc )
			   k = k + 2
			endif
		enddo


	call ncclos(ncid, RCODE)

        return
	end
c
	subroutine specreadzs(datanm)
c
c==> reads the spectral orography file.
c    datanm is the output complex 2016 element array of spectral coefficients.
c
	parameter (mwaves = 63) 
	parameter (nmdim = (mwaves+1)*mwaves/2)
c
        character*80 input_file
        complex datanm(nmdim)
        character*4 what
        integer istart(1),icount(1)
c
        INCLUDE '/usr/local/include/netcdf.inc'
c
        input_file='/Datasets.large/archive0/Datasets.public/'//
     *  'ncep.reanalysis/spectral/orog.spec.nc'

        if (input_file(1:15) .eq. '/Datasets.large' ) then

           print*, "Change the path to the data files!"
           stop 99

        endif

        ncid = ncopn(input_file, NCREAD, IERCODE)
        what = 'orog'
        ispecid = ncvid(NCID,what,IERCODE)
        istart(1) = 1
        icount(1) = 2*nmdim
        call ncvgt(ncid, ispecid, istart, 
     $             icount, datanm, IERCODE)
	call ncclos(ncid,IERCODE)
c
        return
        end        
      program driver
c
c Read a bunch of spectral data.
c

      parameter (nlon = 192)
      parameter (nlat = 94)
      parameter (mwaves = 63)
      parameter (nmdim = (mwaves+1)*mwaves/2)
      parameter (levels = 28)
c
      complex datanm(nmdim,levels)
      character*4 what

      what = 'orog'

      itime_index = 1

      nyear = 1948

      call readspec(datanm,what,itime_index,nyear)

      do i = 1, 10
         write(6,10) what, datanm(i,1)
 10     format(a4,' ---> ',e16.10,'  ',e16.10)
      enddo

      what = 'div'

      call readspec(datanm,what,itime_index,nyear)

      do i = 1, 10
         write(6,10) what, datanm(i,1)
      enddo

      what = 'lnps'

      call readspec(datanm,what,itime_index,nyear)

      do i = 1, 10
         write(6,10) what, datanm(i,1)
      enddo


      what = 'vort'

      call readspec(datanm,what,itime_index,nyear)

      do i = 1, 10
         write(6,10) what, datanm(i,1)
      enddo


      what = 'shum'

      call readspec(datanm,what,itime_index,nyear)

      do i = 1, 10
         write(6,10) what, datanm(i,1)
      enddo


      what = 'vair'

      call readspec(datanm,what,itime_index,nyear)

      do i = 1, 10
         write(6,10) what, datanm(i,1)
      enddo

      stop
      end
