Hello,
I'm a student at the Philipps-University Marburg (Germany) and want to work
with the NCEP/NCAR Reanalysis Data but I've got problems compiling
'readgeneral.f' from CDC, which is for reading netcdf format files from CDC.
I have installed netcdf-3.5.1 and udunits-1.11.7 from source on this system:
Linux linux 2.4.21-99-athlon #1 Wed Sep 24 13:34:32 UTC 2003 i686 athlon i386
GNU/Linux
I set following environments:
setenv CPPFLAGS "-DNDEBUG -Df2cFortran"
setenv FC g77
setenv LD_MATH "-lm"
I tried to compile readgeneral.f with the following command:
g77 readgeneral.F -lnetcdf -ludunits
and get this output:
/usr/local/udunits/include/udunits.inc: In subroutine `gridread':
In file included from readgeneral.F:174:
/usr/local/udunits/include/udunits.inc:27: error: undefined or invalid #
directive
/usr/local/udunits/include/udunits.inc:36:
PTR utmake
1 2
Unrecognized statement name at (1) and invalid form for assignment or
statement-function definition at (2)
readgeneral.F:217: warning:
call ncagt(inet,ivar,'scale_factor',xscale,icode1)
1
readgeneral.F:219: (continued):
....
/usr/local/udunits/include/udunits.inc:36:
PTR utmake
1 2
Unrecognized statement name at (1) and invalid form for assignment or
statement-function definition at (2)
Attached is 'udunits.inc' and 'readgeneral.F'
Thank you very much for your help!
Regards
Katja Trachte
__________________________________________________________
Mit WEB.DE FreePhone mit hoechster Qualitaet ab 0 Ct./Min.
weltweit telefonieren! http://freephone.web.de/?mc=021201 integer UT_EOF, UT_ENOFILE, UT_ESYNTAX, UT_EUNKNOWN
integer UT_EIO, UT_EINVALID, UT_ENOINIT, UT_ECONVERT
integer UT_EALLOC, UT_ENOROOM, UT_ENOTTIME
parameter (UT_EOF = 1)
parameter (UT_ENOFILE = -1)
parameter (UT_ESYNTAX = -2)
parameter (UT_EUNKNOWN = -3)
parameter (UT_EIO = -4)
parameter (UT_EINVALID = -5)
parameter (UT_ENOINIT = -6)
parameter (UT_ECONVERT = -7)
parameter (UT_EALLOC = -8)
parameter (UT_ENOROOM = -9)
parameter (UT_ENOTTIME = -10)
integer UT_MAXNUM_BASE_QUANTITIES
parameter (UT_MAXNUM_BASE_QUANTITIES = 10)
C The FORTRAN API:
C
C NB: `PTR in the following stands for whatever FORTRAN type is
C appropriate for storing a pointer to a structure. Because this is
C necessarily platform-dependent, IT IS UP TO THE USER TO DECLARE AND USE
C THE CORRECT TYPE.
#ifndef PTR
# define PTR integer
#endif
C
C Initialize the units package:
integer utopen
C (character*(*) fpath)
C Create a new unit:
PTR utmake
C ()
C Is a unit a temporal one?
integer uttime
C (PTR unit)
C Indicate whether or not a unit has a non-zero origin (0=>no, 1=>yes).
integer utorigin
C (PTR unit)
C Clear a unit:
C utclr
C (PTR unit)
C Copy a unit:
C utcpy
C (PTR source,
C PTR dest)
C Shift the origin of a unit:
C utorig
C (PTR source,
C doubleprecision amount,
C PTR result)
C Scale a unit:
C utscal
C (PTR source,
C doubleprecision factor,
C PTR result)
C Multiply two units:
C utmult
C (PTR term1,
C PTR term2,
C PTR result)
C Invert a unit:
C utinv
C (PTR source,
C PTR result)
C Divide one unit by another:
C utdiv
C (PTR numer,
C PTR denom,
C PTR result)
C Raise a unit to a power:
C utexp
C (PTR source,
C PTR power,
C PTR result)
C Decode a formatted specification into a unit:
integer utdec
C (character*(*) fspec,
C PTR unit)
C Convert a temporal value into a UTC Gregorian date and time:
integer utcaltime
C (real value,
C PTR unit,
C integer year
C integer month,
C integer day,
C integer hour,
C integer minute,
C real second)
C Convert a Gregorian/Julian date and time into a temporal value:
integer uticaltime
C (integer year,
C integer month,
C integer day,
C integer hour,
C integer minute,
C real second,
C PTR unit,
C doubleprecision value)
C Encode a unit into a formatted specification:
integer utenc
C (PTR unit,
C char fspec)
C Convert from one unit to another:
integer utcvt
C (PTR from,
C PTR to,
C doubleprecision slope,
C doubleprecision intercept)
C Free a unit thingy created by utmake():
C utfree
C (PTR unit)
C Close the units package:
C utcls
C ()
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C This is a sample program that will read netCDF data files
C from the NOAA-CIRES Climate Diagnostics Center.
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/packages/netcdf/index.html
C
C UDUNITS can be found at
C http://www.unidata.ucar.edu/packages/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
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 <your_file>.
C
C STEP 4.
C Change the name of the input file to match the data file you want to read.
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 Climate Diagnostics Center, Boulder, CO.
C based on code writen by Tom Baltzer, Roland Schweitzer and
C Cathy Smtih
C For help, contact Roland Schweitzer or Cathy Smith
C Roland.Schweitzer@xxxxxxxx or
cathy.smith@xxxxxxxx
C :)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C MAIN CODE
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/udunits/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 CDC 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
ntime = 5
ilev = 0
do it=1,ntime
call gridread(inet,ivar,it,x,jlat,ilon,ilev,ny,nm,nd,nh,work)
print*,x(2,2),x(72,36),ny,nm,nd,nh
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/udunits/include/udunits.inc'
include '/usr/local/netcdf/include/netcdf.inc'
********************************************************************
C UDUNITS STUFF TO UNPACK TIME:
********************************************************************
C declare udunits function calls
integer*4 UTMAKE
c integer 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
c unitptr=utmake()
retcode=utdec(unitstr,unitptr)
retcode=uttime(unitptr)
isave=1
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))
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
if((xscale.eq.1).and.(xoff.eq.0))then
ipack=0
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,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)
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 CDC 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 CDC on ???
C Modified by Tom Baltzer of CDC on 2/14/95
C - Broke out into own file and added comments
C Modified by Tom Baltzer of CDC on 2/27/95
C added declarations comments and made to work with
C reanalysis.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine timeindex(netid,nyr,nmn,ndy,nhr,itime)
C STEP 1.
include '/usr/local/netcdf/include/netcdf.inc'
include '/usr/local/udunits/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
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
! udunits.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 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
c integer UTDEC ! Declare the functions.
c 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
if (nyr.gt.1900.and.nyr.lt.3000) then
inyr2d = nyr - 1900
inyr4d = nyr
else if ( nyr.lt.100 ) then ! Between 0-99, it's a 2 digit year.
inyr2d = nyr
inyr4d = nyr + 1900
else if ( nyr .ge. 100 .and. nyr .lt. 2000) then ! Some other year
inyr2d = nyr !Not really, but it's non-zero which is inportant.
inyr4d = nyr !What ever they gave us was the year we want.
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
timevid=NCVID(netid,'time',ercode)
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 CDC 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') monthly = 1
if (cdeltat(10:10).eq.'1') daily = 1
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 (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:nblen(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 CDC 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 CDC on ???
C Modified by Tom Baltzer of CDC 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 CDC
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 CDC on ???
C Modified by Tom Baltzer of CDC 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/udunits/include/udunits.inc'
C Note: POINTER type is integer*4 on the CDC SparcCenter 2000
c system running SunOS 5.3
c 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 CDC 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 CDC 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/udunits/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/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
