program readers c date - 7/31/96 pgmr - Anita C. Brenner c 7/96 Brian D. Beckley c Function -Sample program to read ers1 pathfinder data base and related files c Shows how to read the main sea surface height residual file c and the auxiliary files c Before Running This Program you must do the following c 1) set the names of the files in routine ers1init c 2) see comments in ers1init to change the record lengths c of the files based on the number of repeat cycles c 3) dimension in this main program, issh and iflag by the number c of repeat cycles present in the data base - see README file c This program will print out all available information for the sea surface c height near the johnston isle tide gauge. Note that we already know c the number of the closest track and the indice within that track c for this tide gauge - a georeference directory will be available c shortly to select data by geographic location c The output for this is in readers1.out at the same location c where this source was; please check your output c against this file to make certain you have loaded the correct files implicit none logical*1 lbit16 INTEGER*2 ITRK,IDX,ISSH,iflag,mss,iundf2 integer*4 nreclephm,nreclmss,nreclsshfl,nrecltime,ncycles integer*4 numrevs,nrecldir,nperiod,nperiodp,ncyclesp,numrevsp integer*4 nrectrk,nflag integer*4 ilat,ilon,nrec,i,index,icycle,mjd,nytrk integer*4 ntrack,in5,iout6,inephm,inssh,intime,inflag,inmss,indir real amss,rssh,ssh,stdd,amssh,acc character*80 errflag,ofile dimension lbit16(16) c ers1 specific parameters parameter (nperiodp=6040) parameter (numrevsp=501) c the ncycles parameter need to be defined based on the number of cycles c in the version of the pathfinder data base you are reading parameter (ncyclesp=86) dimension ilat(nperiodp),ilon(nperiodp) dimension nrectrk(nperiodp),mss(nperiodp) double precision fday,xlat,xlatold,xlon,xday,xprev integer*4 iytr(90),iytr2(90),iynb,iynb2,jsign,inyc real yssh(90),yssh2(90) double precision yday(90),yday2(90),ss,ss2,ss3,sss,sss2 c following arrays must be dimensioned by the number of repeats in the c data base - this is a function of data base version dimension issh(ncyclesp),iflag(ncyclesp) c following arrays must be dimensioned by the number of indices in a c reference rev - this is a function of mission c processing area c --------------- real xlatmin,xlatmax,xlonmin,xlonmax,spacing c Namelist/Narea/xlatmin,xlatmax,xlonmin,xlonmax,spacing common/cflag/errflag(16) common/cmission/nreclephm,nreclmss,nreclsshfl,nrecltime, . nrecldir,ncycles,numrevs,nperiod common/constsi2/iundf2 common/cunits/in5,iout6,inephm,inssh,intime,inflag,inmss,indir SAVE call ersinit call openfiles write(*,*) 'Enter min and max latitude, min and max longitude' read(*,*) xlatmin,xlatmax,xlonmin,xlonmax write(*,*) 'Enter spacing between altimetric points (0.1-1.0 deg)' read(*,*) spacing write(*,*) 'Enter output file name ' read(*,*) ofile open(89,file=ofile) c open(90,file='ers2.pos') c write(*,*) xlatmin,xlatmax c c This example initializes the track number, and and indice to be close c to the Johnston isle tide gauge do ntrack = 1,numrevsp c do ntrack = 1,20 write(*,'(A,i6)') ' Now processing track number ',ntrack c Load orbital ephemeris for this track, returns in common block crforb call rtloaders(ntrack,nperiod,ilat,ilon) c load the mean sea surface for this track read(inmss,rec=ntrack)mss c load directory for this track read(indir,rec=ntrack)nrectrk c find first defined index number, if there is no data for index,i c then nrectrk(i) is set to 0 xlatold = 100.0 xprev = 999.9 do 200 i=1,nperiodp if (nrectrk(i).eq.0) go to 200 index=i nrec=nrectrk(i) c load the residual sea surface heights read(inssh,rec=nrec,err=1997)idx,itrk,issh c double check that idx = index and itrk is ntrack if (itrk.ne.ntrack) go to 1998 if (idx.ne.index) go to 1999 c load the flag words read(inflag,rec=nrec,err=2997)idx,itrk,iflag c double check that idx = index and itrk is ntrack if (itrk.ne.ntrack) go to 2998 if (idx.ne.index) go to 2999 c write out location specific quantities; mean sea surface, latitude c and longitude xlat=dfloat(ilat(index))*1.d-6 xlon=dfloat(ilon(index))*1.d-6 amss=dfloat(mss(index))*1.d-2 stdd= 0.01 c write(40,551) itrk,xlat,xlon,amss,stdd 551 format (i6,2f11.5,2f9.3) xprev = xlat c if (xlon.gt.180) xlon = xlon - 360.0 if (xlon.lt.xlonmin.or.xlon.gt.xlonmax) go to 200 if (abs(xlat-xlatold).lt.spacing) go to 200 xlatold = xlat if (xlonmin.lt.0.or.xlonmax.lt.0) then if (xlon.gt.180.0) xlon = xlon-360.0 end if if (xlat.lt.xlatmin.or.xlat.gt.xlatmax) go to 200 c if (abs(xlat).lt.xlatmin) go to 200 if (xlon.lt.xlonmin.or.xlon.gt.xlonmax) go to 200 c write(iout6,10000)ntrack,index c decompose flag word iynb = 0 ss = 0.0 ss2 = 0.0 do 100 icycle=1,ncycles C calculate time call OBSTIME(mjd,fday,iCYCLE,ntrack,index,*80) if (issh(icycle).ne.iundf2) then iynb = iynb + 1 yssh(iynb)=float(issh(icycle))*1.e-3 ss = ss + yssh(iynb) ss2 = ss2 + yssh(iynb)**2 iytr(iynb) = icycle yday(iynb) = (mjd-46066+fday)*86400.0 endif go to 100 80 continue c write(iout6,10300)icycle 100 continue iynb2 = 0 if (iynb.gt.0) then ss=ss/iynb if (iynb.gt.1) ss2=sqrt((ss2/(iynb-1)-ss*ss)) if (iynb.eq.1) then sss = ss sss2 = ss2 iynb2 = 1 yssh2(1) = yssh(1) iytr2(1) = iytr(1) yday2(1) = yday(1) end if if (iynb.ge.2) then iynb2 = 0 sss = 0.0 sss2 = 0.0 ss3 = 3.0*ss2 c if (ss3.gt.0.3) ss3 = 0.3 do icycle = 1,iynb if (abs(yssh(icycle)-ss).lt.ss3) then iynb2 = iynb2 + 1 yssh2(iynb2)=yssh(icycle) sss = sss + yssh(icycle) sss2 = sss2 + yssh(icycle)**2 iytr2(iynb2) = iytr(icycle) yday2(iynb2) = yday(icycle) end if end do if (iynb2.gt.0) sss=sss/iynb2 if (iynb2.gt.1) sss2=sqrt((sss2/(iynb2-1)-sss*sss)) end if if (abs(sss).lt.3.5.and.iynb2.gt.1) then if (sss2.gt.5.0) sss2 = 5.0 acc = 0.20 - 0.15*iynb2/60.0 if (acc.lt.0.05) acc = 0.05 nytrk = ntrack*2 write(89,9991)nytrk,xlat,xlon,amss,sss,acc 9991 format (i9,6f13.5) end if end if 200 continue end do 300 stop 1997 write(iout6,12000)inssh,nrec,ntrack,index stop 16 1998 write(iout6,12100)itrk,ntrack stop 16 1999 write(iout6,12200)idx,index stop 16 2997 write(iout6,12300)inflag,nrec,ntrack,index stop 16 2998 write(iout6,12400)itrk,ntrack stop 16 2999 write(iout6,12500)idx,index stop 16 10000 format('Ocean Pathfinder Output for Track no.',i4,' at indice ', . i5) 10100 format(' latitude (deg N) = ',f8.3,' longitude (deg E) = ',f8.3/ . ' mean sea surface (m) = ',f8.2) 10200 format(' For cycle number ',i4,':' , . 'Mjd = ',i5,' Fday = ',f15.14/5x, . ' sea surface height relative to mean sea surface (m) = ',f8.3/ . 5x,' sea surface height relative to ellipsoid (m) = ',f8.3) 10210 format(' for cycle number ',i4,' the sea surface height is', . ' undefined') 10300 format(' no data available for cycle ',i4,' at this location') 10400 format(' following flags have been set for this sea surface height .:') 10500 format(a) 12000 format('EXECUION TERMINATING'/ . ' i/o error reading residual sea surface heights from', . ' unit;',i3/' trying to read record ',i10,' for track ',i4, . ' at indice ',i5) 12100 format('EXECUION TERMINATING'/ . ' track number on residual sea surface height file is ',i5 . /' this disagrees with the track number expected from the ', . 'directory file',i5/' check that both files are for the same ', . 'data base version') 12200 format('EXECUION TERMINATING'/ . ' indice number on residual sea surface height file is ',i5 . /' this disagrees with the indice number expected from the ', . 'directory file',i5/' check that both files are for the same ', . 'data base version') 12300 format('EXECUION TERMINATING'/ . ' i/o error reading flag word from', . ' unit;',i3/' trying to read record ',i10,' for track ',i4, . ' at indice ',i5) 12400 format('EXECUION TERMINATING'/ . ' track number on flag word file is ',i5 . /' this disagrees with the track number expected from the ', . 'directory file',i5/' check that both files are for the same ', . 'data base version') 12500 format('EXECUION TERMINATING'/ . ' indice number on flag word file is ',i5 . /' this disagrees with the indice number expected from the ', . 'directory file',i5/' check that both files are for the same ', . 'data base version') end