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 psl.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