Right-click or shift-click readgeneral.F to save the file.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C	This is a sample program that will read netCDF data files
C	from NOAA/ESRL PSD.
C
C	In order to run the code you must have the netCDF and UDUNITS
C	libraries installed from UniData.
C
C		netCDF can be found at
C		http://www.unidata.ucar.edu/software/netcdf/index.html
C
C		UDUNITS can be found at
C		http://www.unidata.ucar.edu/software/udunits/index.html
C
C	Don't forget to add them -lnetcdf and -ludunits options to the compile
C	or load command line.
C
C	There are four things that you must change to use this code
C	at your site.  Each place where you must make a change is
C	indicated with a comment of the form like the line below.
C	If you are working on a 64-bit platform, make sure your
C	udunits.inc file defines UD_POINTER to be integer*8.
C
C STEP 1.
C	Change the location of the netcdf.inc and udunits.inc include files
C	to match your installation. (Occurs in 4 places.)
C
C STEP 2.
C	Change the location of the udunits.dat file to match your system
C	installation.
C
C STEP 3.
C	Change PARAMETER values for ilon and jlat to match your file.
C	You can find out the sizes of the lon and lat dimensions by looking
C	at the output of ncdump -h .
C
C STEP 4.
C   Change the name of the input file to match the data file you want to read.
C   This is in the ncopn subroutine call.  Also, change the variable name in
C   the ncvid subroutine call right beneath the ncopn subroutine call.
C
C STEP 5.
C   Choose among the three different ways to read data and remove the comments
C	to get the code you need to do your job.
C
C NOTES:
C	Data is read in one 2D grid at a time in the same lat/lon order
C	it was stored.
C   The gridread subroutine returns the date of the grid read.
C
C     Code written 10/9/96 by Cathy Smith
C
C     NOAA-CIRES ESRL/PSD, Boulder, CO.
C     Based on code writen by Tom Baltzer, Roland Schweitzer and
C     Cathy Smith
C     For help, contact esrl.psd.data@noaa.gov
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C    MAIN CODE
C
c
c This is the latitude, longitude dimension of the grid to be read.
c This corresponds to the lat and lon dimension variables in the netCDF file.
c
C STEP 3.
      parameter (ilon =360)
      parameter (jlat=180)

c
c The data are packed into short integers (INTEGER*2).  The array work
c will be used to hold the packed integers.  The array 'x' will contain
c the unpacked data.
c
C    DATA ARRAY AND WORK ARRAY
      real x(ilon,jlat)
      integer*2 work(ilon,jlat)

c
c Initialize the UDUNITS package.  You may need to change the location
c of the udunits.dat file to match your system.
c
C STEP 2.

       retcode=utopen('/usr/local/etc/udunits.dat')
C STEP 4.
c
c Below in the ncopen call is the file name of the netCDF file.
c You may want to add code to read in the file name and the variable name.
c
c The variable name is the netCDF variable name.  At PSD we name our files
c after the variable name.  If you follow that convention the variable name
c is the same as the first few letters of the file name (the letters up to
c but not including the first '.'.
C
C    OPEN FILE AND GET FILES ID AND VARIABLE ID(S)

      inet=ncopn('sst.mnmean.nc',0,icode)
      ivar=ncvid(inet,'sst',icode)
C STEP 5.
C Pick one of the following scenarios to match your needs.
C
C Comment out this section and remove the comment from one of
C the sections below.

        print*, 'Pick one of the three scenarios for your data.'
        print*, 'See STEP 5 in the comments.'
        stop 'Pick a scenario'
c
C Scenario 1.
C USE THIS LOOP TO READ THE FIRST 'NTIME' TIME STEPS FROM A FILE
C WITH ONLY ONE LEVEL.  IF THE FILE ONLY HAS ONE LEVEL SET ILEV TO 0.
C
C       ntime = 5
C       ilev = 0
C       do it=1,ntime
C          call gridread(inet,ivar,it,x,jlat,ilon,ilev,ny,nm,nd,nh,work)
C          print*,x(2,2),x(72,36),ny,nm,nd,nh
C       enddo
C
C Scenario 2.
C USE THIS LOOP TO READ THE FIRST 'NTIME' TIME STEPS FROM A FILE
C WITH 'NLEV' LEVELS IN IT.  EACH PASS THROUGH THE LOOP WILL RESULT
C IN A GRID AT LEVEL ILEV, AND TIME STEP NTIME.
C

c     ntime = 5
c     nlev = 17
c
c     do it=1,ntime
c       do ilev = 1, nlev
c          call gridread(inet,ivar,it,x,jlat,ilon,ilev,ny,nm,nd,nh,work)
c          print*,x(1,1),ny,nm,nd,nh
c       enddo
c     enddo

c
C Scenario 3.
c THIS CODE BLOCK WILL READ A SPECIFIC TIME AND DATE FROM THE FILE.
c

c
c Set the time of the grid to be read.
c

      ny = 1995
      nm = 3
      nd = 6
      nh = 12

c
c Set this to 0 for a file with one level, loop over it to get all levels,
c or set it to the index of a specific level.
c
c     ilev = 5

c     call timeindex(inet,ny,nm,nd,nh,it)
c     call gridread(inet,ivar,it,x,jlat,ilon,ilev,ny,nm,nd,nh,work)

c
c The value of ny, nm, nd, and nh have been overwritten with the
c values for the grid just read.
c
c     print*,x(1,1),ny,nm,nd,nh

      stop
      end

********************************************************************
********************************************************************
C     MAIN SUBROUTINE TO READ GRID

      subroutine gridread(inet,ivar,nt,x,jlat,ilon,ilev,iyear,imonth,
     &   iday,ihour,idata)

C STEP 1.
      include '/usr/local/include/udunits.inc'
      include '/usr/local/include/netcdf.inc'

********************************************************************
C     UDUNITS STUFF TO UNPACK TIME:
********************************************************************
C     declare udunits function calls

      integer UTMAKE,UTDEC,uttime

C     declare everything else

      integer*4 iyear,imonth,iday,ihour,retcode,itimeid
      integer*4 unitptr
      character*1024 unitstrin
      character*1024 unitstr

      integer*2 idata(ilon,jlat)
      real*4 x(ilon,jlat)
      integer count(3),start(3)
      integer countl(4),startl(4)
      integer icount2(1),istart2(1)
      integer*2 miss
      integer vartype, rtncode, ndims, dims(100), natts
      integer isdfltao, isdfltsf
      character*128 varname
      real*8 xtime


      unitptr=utmake()
      itimeid=NCVID(inet,'time',icode)
      call ncainq(inet,itimeid,'units',iNCCHAR,mlen,i)
      call ncagtc(inet,itimeid,'units',unitstrin,mlen,icode)
      unitstr=unitstrin(1:nblen(unitstrin))
      retcode=utdec(unitstr,unitptr)
      retcode=uttime(unitptr)
      isave=1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C     TEST AND SEE OF FILE HAS LEVELS
C     ALSO SEE IF PACKED

C     GET SCALE FACTOR AND MISSING VALUE IF NECESSARY
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      call NCPOPT(0)	!Turn off netCDF error reaction
      call ncagt(inet,ivar,'scale_factor',xscale,icode1)
      call ncagt(inet,ivar,'add_offset',xoff,icode2)
      call ncagt(inet,ivar,'missing_value',miss,icode3)
C     you need to get level or else use integer level
C     SLAB INFO FOR DATA ARRAY

C     see if packed or not!!!!
      ipack=0
      if((icode1.eq.0).and.(icode2.eq.0))then
        isdfltao = 0
	isdfltsf = 0
	if (abs(xscale - 1.0) .lt. 1.0e-10) then
	    isdfltsf = 1
        endif
	if (abs(xoff) .lt. 1.0e-10) then
	    isdfltao = 1
        endif
	if ((isdfltsf .eq. 1).and.(isdfltao .eq. 1)) then
C        if((xscale.eq.1).and.(xoff.eq.0))then
         call ncvinq(inet, ivar, varname, vartype, ndims, dims, natts,
     +               rtncode)
         if (rtncode .ne. ncnoerr) then
            print *, 'Error calling nf_inq_vartype.'
            stop
         endif
         if (vartype .eq. ncshort) then
            ipack = 1
         else
            ipack = 0
         endif
        else
         ipack=1
        endif
      endif

        start(3)=nt
        if(ilev.eq.0)then
          start(1)=1
          start(2)=1
          start(3)=nt

          count(1)=ilon
          count(2)=jlat
          count(3)=1

C     LOOP THROUGH 1 YEAR OF DATA, UNPACK ARRAY, UNPACK TIME
C     AND READ HEADER

         if(ipack.eq.1)then
           call ncvgt(inet,ivar,start,count,idata,icode)
           call unpack(idata,x,xscale,xoff,miss,ilon,jlat)
         else
           call ncvgt(inet,ivar,start,count,x,icode)
         endif
         call ncvgt(inet,itimeid,it,1,xtime,icode)
         istart2(1)=nt
         icount2(1)=1
         call ncvgt(inet,itimeid,istart2,icount2,xtime,iercode)
         call udparse(unitstrin,xtime,iyear,imonth,iday,ihour)
      else
          startl(1)=1
          startl(2)=1
          startl(3)=ilev
          startl(4)=nt

          countl(1)=ilon
          countl(2)=jlat
          countl(3)=1
          countl(4)=1

C     LOOP THROUGH 1 YEAR OF DATA, UNPACK ARRAY, UNPACK TIME
C     AND READ HEADER

          if(ipack.eq.1)then
            call ncvgt(inet,ivar,startl,countl,idata,icode)
            call unpack(idata,x,xscale,xoff,miss,ilon,jlat)
          else
            call ncvgt(inet,ivar,startl,countl,x,icode)
          endif
          call ncvgt(inet,itimeid,nt,1,xtime,icode)
          istart2(1)=nt
          icount2(1)=1
          call ncvgt(inet,itimeid,istart2,icount2,xtime,iercode)
          call udparse(unitstrin,xtime,iyear,imonth,iday,ihour)
      endif

      return
      end

********************************************************************
      subroutine unpack(x,y,xscale,xadd,miss,ilon,jlat)
      parameter (zmiss=1e30)

      integer*2 x(ilon,jlat),miss
      real y(ilon,jlat)

      do i=1,ilon
         do j=1,jlat
            if(x(i,j).eq.miss)then
               y(i,j)=zmiss
            else
              y(i,j)=(x(i,j)*xscale)+xadd
           endif
         enddo
      enddo

      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      
C
C     This routine (part of the netopen package for PSD netCDF file 
C     handling) takes a netCDF file id, a year, month, day and hour
C     and converts it to the time index for that time in the file
C     and returns that index
C
C     NOTE: WILL NOT HANDLE LESS THAN HOURLY DATA!  Prints error message
C 		and terminates.
C
C     Written by Cathy Smith of PSD on 11/1/94
C     Modified by Tom Baltzer of PSD on 2/14/95
C		- Broke out into own file and added comments
C     Modified by Tom Baltzer of PSD on 2/27/95
CC		  added declarations comments and made to work with 
C		  reanalysis.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine timeindex(netid,nyr,nmn,ndy,nhr,itime)
      include '/usr/local/include/netcdf.inc'
      include '/usr/local/include/udunits.inc'

      integer netid			! NetCDF file id
      integer nyr,nmn,ndy,nhr		! Input year,month,day & hour
      integer itime			! Return time index

      integer inyr2d			! 2 digit in year (used throughout)
      integer inyr4d			! 4 digit input year
      real*8  xtime(10)			! Holds 10 time values from file
      real*8  inxtime			! Input time in netCDF file form
      character*25 cdeltat		! delta_t attribute from file
      character*25 cstat		! Statistic attribute from file
      integer timevid			! Time variable id
      integer timedid			! Time dimension id
      character*100 timeunit 		! unit attribute for time variable
      integer*4 unitptr			! pointer to udunits unit
      integer*4 UTMAKE			! Declare the function.  It's in the
					! uduints.inc at our site, but may 
                                        ! not be at others.
      integer yr1,mn1,dy1,hr1		! First time value in file
      integer yr2,mn2,dy2,hr2		! Second time value in file
      integer*4 xhour1,xhour2,xhourin   ! results of xhour calls 
      integer*4 hourinc			! delta time in hours
      integer monthly			! 1 if data is monthly
      integer seasonal			! 3 if data is monthly
      integer daily			! 1 if data is daily
      integer found			! 0 if not found 1 if found

      integer ltm			! Used to identify a LTM file
      integer equalinc			! Identify non-equal time increments
      
      integer ercode			! To get netCDF error codes
	  integer UTDEC				! Declare the functions.
	  integer UTICALTIME
      integer istart(1),icount(1)	! Used in netCDF calls
      character*15 cname		! Place holder

C     Initialize
c
c I added the condition to find values less that 100.  We assume that such
c a year is a 2 digit year in the 20th century.
c
c The other else if block gets set up for other years, like 1851 the beginning
c of the DOE data set.
c
c     if (nyr.gt.1900.and.nyr.lt.2000) then
c        inyr2d = nyr - 1900
c        inyr4d = nyr
c     else if ( nyr.lt.100 ) then       ! Between 0-99, it's a 2 digit year.
c        inyr2d = nyr
c        inyr4d = nyr + 1900
c     else if ( nyr .ge. 100 .and. nyr .lt. 2000) then  ! Some other year
c        inyr2d = nyr   !Not really, but it's non-zero which is inportant.
c        inyr4d = nyr   !What ever they gave us was the year we want.
c     endif

c
c Do a quick check for ltm, so we can avoid messing with the year value
c if it's a ltm file.
c

      inyr2d = 0
      inyr4d = 0

      timevid=NCVID(netid,'time',ercode)      
      ltm=0
      call NCPOPT(0)
      call NCAGT(netid,timevid,'ltm_range',ltm_range,ercode)
      if (ercode .eq. 0) then
         ltm = 1
      endif
      call NCPOPT(NCVERBOS+NCFATAL)

      if (nyr.ge.1.and.nyr.le.99.and.ltm.eq.0) then
        inyr2d = nyr
        inyr4d = nyr + 1900
        write(0,*) '!*!*!*!*!*!*!*!*!* WARNING *!*!*!*!*!*!*'
        write(0,*) '! Your input of ',nyr,' was converted     !'
        write(0,*) '! to ', inyr4d, '.                           !'
        write(0,*) '! A 4-digit year must be used to open  !'
        write(0,*) '! files for the year 2000 and beyond.  !'
        write(0,*) '! We recommend using 4-digit years     !'
        write(0,*) '! in all cases.                        !'
        write(0,*) '!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*'
      elseif (ltm.eq.0) then
        inyr2d = 99   ! It just needs to be some non-zero value.
        inyr4d = nyr
      endif

      ltm = 0
      equalinc = 0
      monthly = 0
      daily = 0

C     Get all the time variable information including the first 10 or
C     (if < 10) total available values and convert 1st 3 to yr,mn,dy, & hr

      timedid=NCDID(netid,'time',ercode)
      call NCDINQ(netid,timedid,cname,idimt,ercode)
      istart(1)=1
      if (idimt.ge.10) then
         icount(1)=10
      else
         icount(1)=idimt
      endif

      call NCVGT(netid,timevid,istart,icount,xtime,ercode)      
      call NCAGTC (netid,timevid, 'units', timeunit, 100, ercode)

      if (idimt.ge.3) then
         call udparse(timeunit, xtime(1), yr1,mn1,dy1,hr1)
         call udparse(timeunit, xtime(2), yr2,mn2,dy2,hr2)
      elseif (idimt.eq.1) then
         write(0,*) "Only one time increment in the file - returning it"
         itime = 1
         found = 1
         return
      else
         equalinc = 1
         call udparse(timeunit, xtime(1), yr1,mn1,dy1,hr1)
         call udparse(timeunit, xtime(2), yr2,mn2,dy2,hr2)
      endif


C     Turn off error checking and see if this is a PSD data file, and 
C     if so, check time increment < hourly and check for ltm file.
C     Also do first check for equal increment.

      call NCPOPT(0)
      call NCAGTC(netid,timevid,'delta_t',cdeltat,25,ercode)
      if (ercode.eq.0) then
         equalinc = 1
c
c Previously, this was called with variable id hard coded to 1.
c This won't work for multiple variables in the same run.
c
c        call NCAGTC(netid,1,'statistic',cstat,25,ercode)
c
c We will detect long term means from the ltm_range attribute.
c That way we don't need the variable id at all.
c
	call NCAGT(netid,timevid,'ltm_range',ltm_range,ercode)

c
c Ugly kludge to handle different spellings of LTM.  Probably
c should do this as a caseless comparison.  rhs 1/8/96
c Using ltm range gets rid of this ugly kludge.
c
c        if(cstat(1:14).eq.'Long Term Mean' .or.
c    &      cstat(1:14).eq.'Long term mean' .or.
c    &      cstat(1:14).eq.'long term mean') then
c
c If the return code from the above is zero then we have
c and LTM file.
	if (ercode .eq. 0) then
c
c Looking for daily long term means as well as monthly.
c
c           if(inyr2d.ne.0.or.ndy.ne.0.or.nhr.ne.0) then
c              write(0,*) 'timeindex: Warning - non-monthly value ',
c    &                  'requested for monthly LTM file'
c              write(0,*) 'Using month value ',nmn,' to get index'
c           endif
c
c Detect both monthly and daily long term mean.
c
            if(inyr2d.ne.0.or.nhr.ne.0) then
               write(0,*) 'timeindex: Warning - non-zero hourly value ',
     &                  'requested for a LTM file.',
     &                  'Do not know how do deal with hourly LTM.  ',
     &                  'Quiting...'
                stop 'LTM'
            endif
            ltm = 1
         endif
         if (cdeltat(15:16).ne.'00'.or.cdeltat(18:19).ne.'00') then
            write(0,*) 'timeindex: Error: Time increment is < hourly -'
            write(0,*) ' cannot be handled. -Terminating.'
            stop
         endif
         if (cdeltat(7:7).eq.'1') then 
            monthly = 1
         else if (cdeltat(7:7).eq.'3') then 
            seasonal = 1
         else if (cdeltat(10:10).eq.'1') then
            daily = 1
         else if (cdeltat(7:7).ne.'1'.and.cdeltat(7:7).ne.'3'.and.
     $            cdeltat(7:7).ne.'0') then
            write(0,*) 'timeindex: Error:',
     $                 'Cannot handle a delta_t with months increment',
     $                 'other than 0000-01-00 and 0000-03-00.',
     $                 'Your delta_t = ', cdeltat
            stop 'Monthly Delta T'
         endif
      endif
      call NCPOPT(NCVERBOS+NCFATAL)

C     Get the time index.

      if (ltm.eq.1 .and. monthly.eq.1) then
         itime=(nmn-mn1)+1
         found = 1
      elseif ( ltm.eq.1 .and. daily.eq.1) then
         call xhour(yr1,mn1,dy1,hr1,xhour1)
         call xhour(yr2,mn2,dy2,hr2,xhour2)
         hourinc = xhour2 - xhour1
c
c Well, this looks ugly.
c Previously, these xhour values were computed from the year 0.
c xhour turns the year 0 into 1900 and all the calculations are
c done relative to 1900.  1900 is a leap year, which pushes the
c ltm values off by a day.  Instead, since we know it's a ltm
c already from the statistic attribute, we use year 1 for
c the xhour calculations.  rhs 1/8/95
c
         if (daily.eq.1.or.hourinc.eq.24) then
c           call xhour(yr1,mn1,dy1,0,xhour1)
c           call xhour(yr2,mn2,dy2,0,xhour2)
c           call xhour(inyr4d,nmn,ndy,0,xhourin)
            call xhour(1,mn1,dy1,0,xhour1)
            call xhour(1,mn2,dy2,0,xhour2)
            call xhour(1,nmn,ndy,0,xhourin)
         else
            call xhour(inyr4d,nmn,ndy,nhr,xhourin)
         endif
         itime = INT(((xhourin - xhour1)/hourinc) + 0.5)
         itime = itime + 1
         found = 1
      elseif (monthly.eq.1) then
         itime=(((inyr4d*12)+nmn)-((yr1*12)+mn1))
         itime=itime+1
         found = 1
      elseif (seasonal.eq.1) then
         itime=(((inyr4d*12)+nmn)-((yr1*12)+mn1))/3
         itime=itime+1
         found = 1
      elseif (equalinc.eq.1) then
         call xhour(yr1,mn1,dy1,hr1,xhour1)
         call xhour(yr2,mn2,dy2,hr2,xhour2)
         hourinc = xhour2 - xhour1
         if (daily.eq.1.or.hourinc.eq.24) then
            call xhour(yr1,mn1,dy1,0,xhour1)
            call xhour(yr2,mn2,dy2,0,xhour2)
            call xhour(inyr4d,nmn,ndy,0,xhourin)
         else
            call xhour(inyr4d,nmn,ndy,nhr,xhourin)
         endif
         itime = INT(((xhourin - xhour1)/hourinc) + 0.5)
         itime = itime + 1
         found = 1
      else
         write(0,*) 'Cannot assume consistant delta time for ',
     &			'finding index (no delta_t attribute)'
         write(0,*) 'This may take a while ...'
         found = 0

C        Inverse parse the date for comparisons - and move through the 
C        time values searching for the one we want.

         if (timeunit(1:1).eq.'y'.or.timeunit(1:1).eq.'Y') then
            inxtime = inyr4d*10000000000.d0 + nmn*100000000.d0 +
     &                ndy*1000000. + nhr*10000.
         else
            unitptr = UTMAKE()
            ercode = UTDEC(timeunit(1:lnblnk(timeunit)),unitptr)
            if (ercode.ne.0) then
               call uduerr(ercode,'UTDEC','')
               write(0,*) ''
               write(0,*) 'NOTE: You must call netop_init to use 
     &				gridread,'
               write(0,*) '      gridreadx, dayread, and dayreadx'
               stop
            endif
            ercode = UTICALTIME(inyr4d,nmn,ndy,nhr,0,0.,unitptr,inxtime)
            if (ercode.ne.0) then
               call uduerr(ercode,'UTICALTIME','')
               write(0,*) ' '
               write(0,*) 'Something wrong with finding time increment'
               write(0,*) 'Terminating'
               stop
            endif
            call UTFREE(unitptr)
         endif

C        See if time is in first group of time values gotten from file

         do i = 1,icount(1)
            if (xtime(i).eq.inxtime) then
               itime = i
               found = 1
            endif
         enddo
         if (icount(1).eq.idimt) then
            write(0,*) 'Could not find desired time increment'
            write(0,*) ' - Terminating'
            stop
         endif

C        Step through file getting groups of values and seeing if it's in there

         istart(1) = istart(1) + icount(1)
         do while (found.eq.0.and.istart(1)+icount(1).le.idimt)
            call NCVGT(netid,timevid,istart,icount,xtime,ercode)      
            do i = 1,icount(1)
               if (xtime(i).eq.inxtime) then
                  itime = (istart(1) + i) - 1
                  found = 1
               endif
            enddo
         enddo

C        Check the last group of time values in the file 
         if (found.eq.0) then
            icount(1) = (idimt - istart(1)) + 1
            call NCVGT(netid,timevid,istart,icount,xtime,ercode)
            do i = 1,icount(1)
               if (xtime(i).eq.inxtime) then
                  itime = (istart(1) + i) - 1
                  found = 1
               endif
            enddo
         endif

         if (found.eq.0) then
            write(0,*) 'Could not find desired time increment'
            write(0,*) ' - Terminating'
            stop
         endif

      endif

      if(itime.gt.idimt)then
         write(0,*)'You are requesting a time past the end of the file'
         istart(1)=idimt
         icount(1)=1
         call NCVGT(netid,timevid,istart,icount,xtime,ercode)
         call udparse(timeunit,xtime(1),iyear,imonth,iday,ihour)         
         write(0,*)' last day of file is: ',iyear,imonth,iday,ihour
      endif
      if(itime.le.0)then
        print*,'you asked for a time before first day of the dataset'
        print*,iyear,imonth,iday,ihour,' is first day of data'         
        itime=0
      endif

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  This subroutine is part of the netopen/read set of Fortran callable
C  routines for reading PSD netCDF files and returning grids.
C
C  This routine takes in a year, month, day and hour and calculates
C  the number of hours since a date before all our data (year 1 month 1
C  day 1) - this gives a referece value from which differences between
C  dates can be derived.
C
C  Note that this routine will take a 4 or a 2 digit date, if it
C  receives a 2 digit date it assumes that it is in the range 1900 -
C  1999.
C
C  Written by Cathy Smith of PSD on 11/1/94
C  Modified by Tom Baltzer of PSD on Feb 28, 1995
C 	- converted from days to hours so that hourly data can be
C 	  accounted for, and made starting date 1-1-1 from 1-1-1920.
C
C  File: xhour.f (to be included in netopen.subr.f)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine xhour(iy,im,id,ih,xhours)

      integer iy,im,id,ih	! The input year, month, day and hour
      integer*4 xhours		! OUT # of hours between in and 1-1-1

      integer*4 xdays		! Used in calculating hours
      integer   inyear		! Used to work with input year

      integer imonth(12)	! Used in calculating # days in this year
      integer imonthl(12)	! "					"

      data imonth/31,28,31,30,31,30,31,31,30,31,30,31/
      data imonthl/31,29,31,30,31,30,31,31,30,31,30,31/

C     See if date is in 1900s but given as 2 digit.

      if (iy.lt.100) then
         inyear = iy + 1900
      else
         inyear = iy
      endif

C     CALCULATE DAYS FROM JAN 1 Year 1

      xdays=0
      xhours=0

      xdays = INT((inyear-1)*365.25)

      if(im.gt.1)then
        do imm=1,im-1
           if((mod(inyear,4).eq.0).and.(inyear.ne.0))then
             xdays=xdays+imonthl(imm)
           else
             xdays=xdays+imonth(imm)
           endif
        enddo
      endif

      xdays=xdays+id

      xhours = xdays*24

      xhours = xhours + ih

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     This subroutine is part of the netopen Fortran callable PSD
C     group of routines for reading generic netCDF files in the new
C     cooperative standard as presented in:
C	http://ferret.wrc.noaa.gov/noaa_coop/coop_cdf_profile.html
C
C     This routine takes a time value pulled from the netCDF file,
C     determines if it is in the new or old time format (by using the
C     udunits function) and then returns the representative year,
C     month, day and hour represented by the time value.
C
C     Written by Cathy Smith of PSD on 11/1/94
C     Modified by Tom Baltzer of PSD on Feb. 8, 1995
C	- To determine if time is old or new format and parse
C		accordingly
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine udparse(timeunit,intime,oyear,omonth,oday,ohour)

C STEP 1.
      include '/usr/local/include/udunits.inc'

C     Note: POINTER type is integer*4 on the PSD SparcCenter 2000
C     system running SunOS 5.3
      integer UTMAKE,UTDEC,utcaltime

      integer*4 unitptr ! Pointer to udunits "unit" type

      character*100 timeunit
			! The unit attribute for time in the netCDF file
      real*8 intime,xx	! The input time value and var to work with it
      integer oyear,omonth,oday,ohour
			! The output year,month,day and hour
      integer tmin	! Temporary storage for minutes value
      real    tsec	! Temporary storage for seconds value
      integer ercode    ! For determining udunits result

      unitptr = UTMAKE()
      ercode = UTDEC(timeunit(1:nblen(timeunit)),unitptr)
      if (ercode.ne.0) then

C        Assume old PSD standard if time unit is unknown

         if (ercode.eq.UT_EUNKNOWN.and.(timeunit(1:1).eq.'y'.or.
     &	     	timeunit(1:1).eq.'Y')) then
            xx=0.
            oyear=int(intime/10000000000.d0)
            xx=oyear*10000000000.d0
            omonth=int((intime-xx)/100000000.d0)
            xx=xx+omonth*100000000.
            oday=int((intime-xx)/1000000.)
            xx=xx+oday*1000000.
            ohour=int((intime-xx)/10000.)
         else
            call uduerr(ercode,'UTDEC','')
            write(0,*) ''
            write(0,*) 'NOTE: You must call netop_init to use gridread,'
            write(0,*) '      gridreadx, dayread, and dayreadx'
            stop
         endif
      else
         ercode = UTCALTIME(intime,unitptr,oyear,omonth,oday,ohour,
     &                      tmin,tsec)
         if (ercode.ne.0) then
            call uduerr(ercode,'UTCALTIME','')
            stop
         endif
      endif

C     Free up the pointer
      call UTFREE(unitptr)

      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  This function is intended to handle all possible error conditions
C  for the Udunits package  (Version 1.7.1) from Unidata, it takes the
C  error code returned from a udunits function, the function name and
C  the filename the function was accessing (applicable only for utopen
C  function).  It prints an indication to the user of the error and
C  function in which it occurred.
C
C  Written by: Tom Baltzer of PSD  February 10, 1995
C
C  File:  uduerr.f  (To be added to netopen.subr.f)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine uduerr(ercode, funcname, filename)

C STEP 1.
      include '/usr/local/include/udunits.inc'

      integer ercode 		! The Udunits Error number
      character* (*) funcname 	! Function name which returned error
      character* (*) filename	! Name of file on which funct operates
      character*1024 path	! Path to file on which funct really works

      if (filename.eq.'') then
         call getenv('UDUNITS_PATH',path)
         if (path.eq.'') then
            path = '/usr/local/src/udunits'
         endif
      else
         path = filename
      endif

      write(0,*) 'Error in udunits function: ',funcname
      if (ercode.eq.UT_ENOFILE) then
         write(0,*) 'File: ',path(1:nblen(path))
         write(0,*) 'Does not exist!'
      elseif (ercode.eq.UT_ESYNTAX) then
         write(0,*) 'A Udunits syntax error was detected'
         if (path.ne.'') then
            write(0,*) 'Possibly in the file: ',
     &			path(1:nblen(path))
         endif
      elseif (ercode.eq.UT_EUNKNOWN) then
         write(0,*) 'Udunits unknown specification found'
         if (path.ne.'') then
            write(0,*) 'Possibly in the file: ',
     &			path(1:nblen(path))
         endif
      elseif (ercode.eq.UT_EIO) then
         write(0,*) 'I/O Error detected'
         if (path.ne.'') then
            write(0,*) 'Possibly in the file: ',
     &			path(1:nblen(path))
         endif
      elseif (ercode.eq.UT_EINVALID) then
         write(0,*) 'An invalid udunits structure was found'
         write(0,*) ' Note: probably a temporal struct'
      elseif (ercode.eq.UT_ENOINIT) then
         write(0,*) 'Udunits package has not been initialized'
      elseif (ercode.eq.UT_ECONVERT) then
         write(0,*) 'No conversion between the two units is possible'
      elseif (ercode.eq.UT_ENOROOM) then
         write(0,*) 'Not enough room in character string for conversion'
      elseif (ercode.eq.UT_EALLOC) then
         write(0,*) 'Memory allocation falure'
      else
         write(0,*) 'UTOPEN: Unknown error: ',ercode
      endif

      return
      end
      integer function nblen (string)
c
c     given a character string, nblen returns the length of the string
c     to the last non-blank character, presuming the string is left-
c     justified, i.e. if string = '   xs  j   ', nblen = 8.
c
c     called non-library routines: none
c     language: standard fortran 77
c
      integer ls,i
      character*(*) string, blank*1, null*1
      data blank /' '/
c
      null = char(0)
      nblen = 0
      ls = len(string)
      if (ls .eq. 0) return
      do 1 i = ls, 1, -1
         if (string(i:i) .ne. blank .and. string(i:i) .ne. null) go to 2
    1    continue
      return
    2 nblen = i
      return
      end