c ****************************************************************** c DESCRIPTION c This program uses the netCDF subroutines to read the SuomiNet GPS data c from the sgp30suomigps*.cdf and gec30suomigps*.cdf files c and output data from a specified station. c c SYNOPSIS c To compile: f77 -g -o suomi suomi.f -lnetcdf c c To execute: suomi < path/filename.cdf > < ARM SGP EF > c or: suomi < path/filename.cdf > 0 < site > c where c < path/filename.cdf > is the path and name of the netCDF file c < ARM SGP EF > is the number of the ARM/SGP Extended Facility where c the desired Suominet station is located. c < site > is the station ID or 'barrow', 'manus', 'nauru', or 'darwin' c c EXAMPLES c To read data from SGP EF 13: c suomi sgp30suomigpsX1.c1.20041203.000000.cdf 13 c To read data from TWP Nauru: c suomi gec30suomigpsX1.c1.20041203.000000.cdf 0 nauru c c LIBRARIES c libnetcdf.a must have been previously installed. This c library may be acquired via anonymous ftp to c unidata.ucar.edu or via the fink package manager c if you had the foresight to purchase a Macintosh. c c AUTHOR c Jim Liljegren c Argonne National Laboratory c Decision and Information Sciences Division, Bldg. 900 c 9700 South Cass Avenue c Argonne, IL 60439 c c ***************************************************************** program suomi include '/sw/include/netcdf.inc' character*280 filename character*9 string character*31 name integer cdfid, rcode integer start(2)/1,1/, count(2) integer itime, tarray(9) parameter (igps0=315964800) ! seconds from 1 Jan 1970 to 6 Jan 1980 parameter (nrec=48,maxstation=130,nchar_station_name=4) character station(maxstation)*4 real*8 lat(maxstation), lon(maxstation), height(maxstation) real*8 offset(nrec) real pwv(nrec,maxstation), pwv_err(nrec,maxstation) real wet_delay(nrec,maxstation),dry_delay(nrec,maxstation) real pifact(nrec,maxstation) real pres(nrec,maxstation),temp(nrec,maxstation),rh(nrec,maxstation) logical missing parameter (nsgpef=27) character sgpef(nsgpef)*4, site*4, twpcf(3)*4,nsacf*4 data sgpef/'SG12','None','SG15','None','SG13','SG14','SG16','SG11','SG04', & 'None','SG10','SG08','SG01','None','SG09','None','None','None', & 'SG20','None','None','SG19','None','SG18','SG17','None','SG34'/ data twpcf/'SA42','SA40','SA39'/ data nsacf/'SG27'/ c Get the file name call getarg(1,filename) c Get the SGP Extended Facility number of the Suominet station c or '0' to indicate the Suominet station ID or the TWP/NSA site name will be input instead. call getarg(2,string) read(string,*) ief c Check that there is a Suominet station at this EF or assign the station name if ( ief .gt. nsgpef .or. sgpef(ief) .eq. 'None' ) then write(*,'("No SuomiNet station at EF-",i2.2)') ief stop else if ( ief .eq. 0 ) then call getarg(3,site) if ( site(1:1).eq.'B' .or. site(1:1).eq.'b' ) then site = nsacf else if ( site(1:1).eq.'M' .or. site(1:1).eq.'m' ) then site = twpcf(1) else if ( site(1:1).eq.'N' .or. site(1:1).eq.'n' ) then site = twpcf(2) else if ( site(1:1).eq.'D' .or. site(1:1).eq.'d' ) then site = twpcf(3) endif else site = sgpef(ief) endif c Make errors non-fatal to allow us to do error handling. c The default is all errors are fatal = call ncpopt(NCVERBOS+NCFATAL) c call NCPOPT(0) c call NCPOPT(NCVERBOS) c Open the netCDF file as read-only NCNOWRIT; get its "ID" CDFID: c cdfid = NCOPN(filename(1:lnblnk(filename)), NCNOWRIT, rcode) if ( rcode.NE.0 ) stop 'Error opening file' c Get the ID of the Station dimension and the number of stations idsta = NCDID(cdfid, 'station', rcode) call NCDINQ(cdfid, idsta, name, nstation, rcode) if ( nstation .gt. maxstation ) stop c Link variable names and variable IDs: c idsta = NCVID(cdfid, 'station', rcode) idlat = NCVID(cdfid, 'lat', rcode) idlon = NCVID(cdfid, 'lon', rcode) idhgt = NCVID(cdfid, 'height', rcode) idoffset = NCVID(cdfid, 'time_offset', rcode) idpwv = NCVID(cdfid, 'pwv', rcode) idpwve = NCVID(cdfid, 'pwv_err', rcode) idwdel = NCVID(cdfid, 'wet_delay', rcode) idddel = NCVID(cdfid, 'final_dry_delay', rcode) idpif = NCVID(cdfid, 'pifact', rcode) idpres = NCVID(cdfid, 'pres', rcode) idtemp = NCVID(cdfid, 'temperature', rcode) idrh = NCVID(cdfid, 'rh', rcode) c Get station info count(1) = nchar_station_name count(2) = nstation nchar = nchar_station_name * nstation call NCVGTC(cdfid, idsta, start, count, station, nchar, rcode) count(1) = nstation call NCVGT(cdfid, idlat, start, count, lat, rcode) call NCVGT(cdfid, idlon, start, count, lon, rcode) call NCVGT(cdfid, idhgt, start, count, height, rcode) c write(*,'(a,2x,3f10.4)') (station(i),lat(i),lon(i),height(i),i=1,nstation) c Determine which station's data to output do i = 1,nstation if ( station(i) .eq. site ) then ista = i c write(*,110) ief,station(ista),lat(ista),lon(ista),height(ista) 110 format(i2,1x,a,2x,3f10.4) missing = .false. goto 1 endif enddo c write(*,111) site, ief 111 format("Station ",a," at EF-",i2.2, " is missing in data file") missing = .true. c Get data 1 count(1) = nrec count(2) = nstation call NCVGT(cdfid, idoffset, start, count, offset, rcode) call NCVGT(cdfid, idpwv, start, count, pwv, rcode) call NCVGT(cdfid, idpwve, start, count, pwv_err, rcode) call NCVGT(cdfid, idwdel, start, count, wet_delay, rcode) call NCVGT(cdfid, idddel, start, count, dry_delay, rcode) call NCVGT(cdfid, idpif, start, count, pifact, rcode) call NCVGT(cdfid, idpres, start, count, pres, rcode) call NCVGT(cdfid, idtemp, start, count, temp, rcode) call NCVGT(cdfid, idrh, start, count, rh, rcode) c Get the start_time (in seconds since 6 January 1980) call NCAGTC (cdfid, NCGLOBAL, 'start_time', string, 9, rcode) if ( rcode .eq. 0 ) then read(string,'(i9)') istart c write(*,'(a,2x,i9)') string, istart else write(*,*) "Error reading start_time" stop endif c Close the netCDF file call NCCLOS(cdfid, rcode) c Convert GPS time to something useful and output data ibase = istart + igps0 do i = 1,nrec itime = ibase + offset(i) call gmtime(itime,tarray) isec = tarray(1) imin = tarray(2) ihr = tarray(3) iday = tarray(4) imth = tarray(5)+1 ! January = 0 and December = 11 iyr = tarray(6)+1900 ! 4-digit year idayoftheweek = tarray(7) ! Sunday = 0 idoy = tarray(8)+1 ! Day of year (Jan 1 = day 0) doy = float(iday) + float( isec + 60*( imin + 60*ihr ) ) / 86400. if ( .not. missing ) then write(*,113) itime, doy, & pres(i,ista), temp(i,ista), rh(i,ista), & 0.1*pwv(i,ista), 0.1*pwv_err(i,ista), & wet_delay(i,ista), pifact(i,ista), dry_delay(i,ista) & else write(*,113) itime, doy, & -9999., -9999., -9999., & -9999., -9999., & -9999., -9999., & -9999. endif enddo stop 112 format(i4.4,2i2.2,1x,i3,1x,i2.2,':',i2.2,':',i2.2,f10.5) 113 format(i10,9('\t',f11.5)) end