NOTICE: This version of the NSF Unidata web site (archive.unidata.ucar.edu) is no longer being updated.
Current content can be found at unidata.ucar.edu.
To learn about what's going on, see About the Archive Site.
NOTICE: This version of the NSF Unidata web site (archive.unidata.ucar.edu) is no longer being updated.
Current content can be found at unidata.ucar.edu.
To learn about what's going on, see About the Archive Site.
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
netcdfgroup archives: