program manyrins *** merge many rinex obs files to array form implicit double precision(a-h,o-z) parameter(MAXOBS=80000000) !*** max obs character*1 letter lin =1 !*** reserved for subroutine use list=2 lout=3 open(lout,file='manyrins.txt',form='formatted') write(*,'(a)') 'program manyrins -- 2012jul18' call readargs(itstrt,itstop,itstep) *** process file list and recon open(list,file='list.txt',form='formatted',status='old') call loadlist(list,lin) call getnfiles(nfiles) write(*,*) 'opened files: ',nfiles *** process files again and store in memory call storeobs(lin,itstrt,itstop,itstep,isodmn,isodmx) call numlng(isize) write(*,'(a,f8.5,a)') * 'program memory use: ',100*dble(isize)/MAXOBS,' precent' *** now write memory to output write(*,*) 'now writing output file ...' call dump1(lout,isodmn,isodmx,itstep) call dump2(lout,isodmn,isodmx,itstep) close(lout) *** go into suspend state until dismissed write(*,*) write(*,'(''Press ENTER key to stop program: '',$)') read(*,'(a)') letter end ************************************************************************ subroutine dumparray(lout,isod,obuf,aobs,lmsksat,lmsksta, * nsta,nsat,noutsta,noutsat) *** dump buffer to output implicit double precision(a-h,o-z) character*2 aobs dimension obuf(nsta,nsat) logical lmsksat(*),lmsksta(*) !*** vectors for condense integer*8 key8,lkey logical getlng parameter(MAXOBS=80000000) !*** max obs dimension odb(MAXOBS) save /dbase/ common/dbase/odb parameter(MXSTA = 1000) !*** max stations dimension dhs(MXSTA),des(MXSTA),dns(MXSTA) character*78 afiles(MXSTA) character*1 asys(MXSTA) character*4 anames(MXSTA) character*20 antnames(MXSTA) save /master/ common/master/dhs,des,dns,nfiles,afiles,asys,anames,antnames,leap character*3 msats !*** master sat table guts common/msttb1/msats(500) if(noutsta.eq.0.or.noutsat.eq.0) return *** 'f' -- obs type record write(lout,'(a,a2,2i5)') 'f ',aobs,noutsta,noutsat *** 't' -- time/station records do ir=1,nsta if(lmsksta(ir)) then key8=lkey(ir,isod,'ZZZ','ZZ') if(.not.getlng(key8,kkey)) then ***** should not enter this clause write(lout,'(a,a,a)') 't ',anames(ir),' no time available' else write(lout,'(a,a,f11.7)') 't ',anames(ir),odb(kkey) endif endif enddo *** 's' -- satellite record write(lout,'(a1$)') 's' do jc=1,nsat if(lmsksat(jc)) write(lout,'(1x,a3$)') msats(jc) enddo write(lout,*) !*** end of line *** dump buffer ('o' records) do ir=1,nsta if(lmsksta(ir)) then write(lout,'(a1$)') 'o' do jc=1,nsat if(lmsksat(jc)) write(lout,'(1x,f14.3$)') obuf(ir,jc) enddo write(lout,*) !*** end of line endif enddo return end ************************************************************************ subroutine makearray(isod,aobs,obuf,lmsksat,lmsksta, * nsta,nsat,noutsta,noutsat) *** dump memory to buffer implicit double precision(a-h,o-z) character*2 aobs dimension obuf(nsta,nsat) integer*8 key8,lkey logical getlng parameter(ONULL=9999999999.999d0) parameter(ONUL2=9999999999.990d0) !*** for checking float.pt. logical lmsksat(*),lmsksta(*) !*** vectors for condense character*3 msats !*** master sat table guts common/msttb1/msats(500) parameter(MAXOBS=80000000) !*** max obs dimension odb(MAXOBS) save /dbase/ common/dbase/odb *** populate buffer do ir=1,nsta do jc=1,nsat key8=lkey(ir,isod,msats(jc),aobs) if(.not.getlng(key8,kkey)) then obuf(ir,jc)=ONULL else obuf(ir,jc)=odb(kkey) endif enddo enddo *** load mask vectors for buffer condensation *** row first noutsta=nsta do ir=1,nsta lmsksta(ir)=.true. !*** set station valid do jc=1,nsat if(obuf(ir,jc).lt.ONUL2) go to 100 !*** live data enddo lmsksta(ir)=.false. !*** fell thru loop - no data noutsta=noutsta-1 100 continue enddo *** now by columns noutsat=nsat do jc=1,nsat lmsksat(jc)=.true. !*** set sat valid do ir=1,nsta if(obuf(ir,jc).lt.ONUL2) go to 200 !*** live data enddo lmsksat(jc)=.false. !*** fell thru loop - no data noutsat=noutsat-1 200 continue enddo return end ************************************************************************ subroutine dump2(lout,isodmn,isodmx,itstep) *** dump memory to output implicit double precision(a-h,o-z) parameter(MXSTA = 1000) !*** max stations dimension dhs(MXSTA),des(MXSTA),dns(MXSTA) character*78 afiles(MXSTA) character*1 asys(MXSTA) character*4 anames(MXSTA) character*20 antnames(MXSTA) save /master/ common/master/dhs,des,dns,nfiles,afiles,asys,anames,antnames,leap save /master2/ common/master2/fmjdm,mjdm,myr,mmn,mdy character*2 mobss !*** master obs table guts common/mobtb1/mobss(50) common/mobtb2/look(47),iptrs(50),nkey,max,ifree !*** do not delete!! parameter(MXSAT = 100) !*** max sats parameter(nbuf = MXSAT*MXSTA) dimension obuf(nbuf) logical lmsksat(MXSAT),lmsksta(MXSTA) !*** vectors for condense call nummst(nmsats) nmstas=nfiles do isod=isodmn,isodmx,itstep !*** outer loop over time steps i =isod ihr =i/3600 i =i-3600*ihr imn =i/60 isec=i-60*imn write(lout,'(a,i4,5i3)') 'e ',myr,mmn,mdy,ihr,imn,isec do iotyp=1,ifree-1 !*** 2nd loop over obs types call makearray(isod,mobss(iotyp),obuf,lmsksat,lmsksta, * nmstas,nmsats,noutsta,noutsat) call dumparray(lout,isod,obuf,mobss(iotyp),lmsksat,lmsksta, * nmstas,nmsats,noutsta,noutsat) enddo !*** loop over all obs types enddo !*** loop over all epochs return end ************************************************************************ subroutine storeobs(lin,itstrt,itstop,itstep,isodmn,isodmx) *** loop over rinex files and store data implicit double precision(a-h,o-z) character*80 card character*78 fname parameter(MXSTA = 1000) !*** max stations dimension dhs(MXSTA),des(MXSTA),dns(MXSTA) character*78 afiles(MXSTA) character*1 asys(MXSTA) character*4 anames(MXSTA) character*20 antnames(MXSTA) save /master/ common/master/dhs,des,dns,nfiles,afiles,asys,anames,antnames,leap save /master2/ common/master2/fmjdm,mjdm,myr,mmn,mdy call newlng !*** initialize long table for data base support call newmst !*** initialize master satellite table write(*,'(a,a)') 'processing files ...' isodmn=+1000000 isodmx=-1000000 do i=1,nfiles open(lin,file=afiles(i),form='formatted',status='old') call rinex1(lin,itstrt,itstop,itstep,i,isodmn,isodmx) close(lin) enddo write(*,*) 'time range: ',isodmn,isodmx *** finalize the reference date if(fmjdm.gt.0.9999d0) then !*** got a piece of yesterday mjdm = mjdm+1 fmjdm= 0.d0 call mjdciv(mjdm,fmjdm,myr,mmn,mdy,mhr,mmn,secm) endif write(*,*) 'ref. date : ',myr,mmn,mdy call nummst(nmsats) write(*,*) 'total sats: ',nmsats return end *----------------------------------------------------------------------- subroutine rinex1(lin,itstrt,itstop,itstep,kname,isodmn,isodmx) *** read a rinex file and store data implicit double precision(a-h,o-z) parameter(ONULL=9999999999.999d0) character*80 card character*1 atim !*** time system code character*2 aotyp(50),typbuf(9) character*3 asats(999),satbuf(12) dimension obss(50),obsbuf(5) logical lallow(50) !*** flag for allowable obs types logical ltick !*** flag for time at interval step logical putlng,getlng,getmob,putmst integer*8 key8,lkey parameter(MXSTA = 1000) !*** max stations dimension dhs(MXSTA),des(MXSTA),dns(MXSTA) character*78 afiles(MXSTA) character*1 asys(MXSTA) character*4 anames(MXSTA) character*20 antnames(MXSTA) save /master/ common/master/dhs,des,dns,nfiles,afiles,asys,anames,antnames,leap parameter(MAXOBS=80000000) !*** max obs dimension odb(MAXOBS) save /dbase/ common/dbase/odb *** process header rewind(lin) 100 read(lin,'(a)') card if(card(61:73).ne.'END OF HEADER') then if(card(61:79).eq.'# / TYPES OF OBSERV') then read(card,'(i6)') nobs if(nobs.gt.50) then write(*,*) 'FATAL -- too many obs in each epoch' stop 45787 endif n=nobs j=0 101 if(n.le.9) then read(card,'(6x,9(4x,a2))') (typbuf(i),i=1,n) do i=1,n j=j+1 aotyp(j)=typbuf(i) !*** obs type enddo else read(card,'(6x,9(4x,a2))') (typbuf(i),i=1,9) do i=1,9 j=j+1 aotyp(j)=typbuf(i) !*** obs type enddo n=n-9 read(lin,'(a)') card !*** obs type continuation rec go to 101 endif if(j.gt.50) stop 12834 !*** should never have to stop here endif go to 100 endif *** set up quick access array for obs type lookup, later do jobs=1,nobs if(.not.getmob(aotyp(jobs),imobs)) then lallow(jobs)=.false. else lallow(jobs)=.true. !*** in master obs type table endif enddo *** process data -- loop over epochs of data *** should be looking at an obs header or an event record 200 read(lin,'(a)',end=777) card read(card,'(28x,i1)') ievntflg if(ievntflg.gt.1) then !*** event read(card,'(29x,i3)') iskip do i=1,iskip read(lin,'(a)') card !*** bypass additional event recs enddo go to 200 !*** reset back to another chunk else !*** epoch/header record read(card,'(5i3,f11.7,3x,i3)') iyr,imo,idy,ihr,imn,sec,nsats if(asys(kname).eq.'R') then !*** GLONASS time system --> GPS sec=sec+leap endif sod = sec + 60.d0*imn + 3600.d0*ihr isod=idnint(sod) !*** integer seconds-of-day if(mod(isod,itstep).eq.0.and. * isod.ge.itstrt.and.isod.le.itstop) then if(isod.lt.isodmn) isodmn=isod if(isod.gt.isodmx) isodmx=isod ltick=.true. else ltick=.false. !*** skip parse and buffer transfer endif if(ltick) then key8=lkey(kname,isod,'ZZZ','ZZ') !*** reserved keys for sod if(.not.putlng(key8,idup)) then if(idup.le.0) then write(*,*) 'FATAL - long table overflow' else write(*,*) 'FATAL - duplicate key elements' endif stop 65600 else if(.not.getlng(key8,kkey)) stop 67800 odb(kkey)=sec !*** store in data base endif endif n=nsats j=0 201 if(n.le.12) then if(ltick) then read(card,'(32x,12(a3))') (satbuf(i),i=1,n) do i=1,n j=j+1 asats(j)=satbuf(i) !*** sats at epoch if(.not.putmst(asats(j),idup)) continue enddo endif !*** end of tick clause else if(ltick) then read(card,'(32x,12(a3))') (satbuf(i),i=1,12) do i=1,12 j=j+1 asats(j)=satbuf(i) !*** sats at epoch if(.not.putmst(asats(j),idup)) continue enddo endif !*** end of tick clause n=n-12 read(lin,'(a)') card !*** obs header continuation rec go to 201 endif ***** data following epoch header record do isat=1,nsats read(lin,'(a)') card !*** a data record n=nobs j=0 202 if(n.le.5) then if(ltick) then read(card,'(5(f14.3,2x))') (obsbuf(i),i=1,n) do i=1,n j=j+1 if(card(1+(i-1)*16:14+(i-1)*16).eq.' ')then obss(j)=ONULL !*** empty obs field else obss(j)=obsbuf(i) !*** obs at epoch endif enddo endif !*** end of tick clause else if(ltick) then read(card,'(5(f14.3,2x))') (obsbuf(i),i=1,5) do i=1,5 j=j+1 if(card(1+(i-1)*16:14+(i-1)*16).eq.' ')then obss(j)=ONULL !*** empty obs field else obss(j)=obsbuf(i) !*** obs at epoch endif enddo endif !*** end of tick clause n=n-5 read(lin,'(a)') card !*** another data record go to 202 endif *** here we store obss(nobs) into memory if(ltick) then do jobs=1,nobs if(lallow(jobs)) then key8=lkey(kname,isod,asats(isat),aotyp(jobs)) if(.not.putlng(key8,idup)) then if(idup.le.0) then write(*,*) 'FATAL - long table overflow' else write(*,*) 'FATAL - duplicate key elements' endif stop 65656 else if(.not.getlng(key8,kkey)) stop 67890 odb(kkey)=obss(jobs) !*** store in data base endif endif !*** end of allowed obs clause enddo !*** end of loop over all obs present endif !*** end of tick clause enddo !*** end of sat loop *** endif !*** end of rinex event/data chunk clause go to 200 *** end of file encountered 777 continue return end ************************************************************************ ************************************************************************ ************************************************************************ subroutine loadlist(list,lin) *** process file name list implicit double precision(a-h,o-z) character*80 card character*78 fname character*1 asystime parameter(tiny = 1.d-10) !*** a tiny number parameter(MXSTA = 1000) !*** max stations dimension dhs(MXSTA),des(MXSTA),dns(MXSTA) character*78 afiles(MXSTA) character*1 asys(MXSTA) character*4 anames(MXSTA) character*20 antnames(MXSTA) save /master/ common/master/dhs,des,dns,nfiles,afiles,asys,anames,antnames,leap save /master2/ common/master2/fmjdm,mjdm,myr,mmn,mdy call newmob !*** initialize master obs code table leap=0 !*** initialize leap sec corrn.. nfiles=0 1 read(list,'(a)',end=777) fname if(fname.eq.' ') go to 1 !*** bypass blank lines if(fname(1:1).eq.'*') go to 1 !*** bypass comments if(access(fname,'r').eq.0) then open(lin,file=fname,form='formatted',status='old') nfiles=nfiles+1 if(nfiles.gt.MXSTA) then write(*,*) 'file name:' write(*,*) fname write(*,*) 'exceeds maximum number of input files: ', MXSTA write(*,*) 'FATAL' stop 98347 endif *** test the version -- first card in rinex obs file is version record read(lin,'(a)') card if(card(61:80).eq.'RINEX VERSION / TYPE') then read(card,'(f9.2,31x,a1)') version,asystime if(asystime.eq.' ') asystime='G' !*** sat time system if(version.gt.(2.12d0 + tiny)) then write(*,*) 'version number:' write(*,*) card(1:9) write(*,*) 'exceeds programming.' write(*,*) 'file name:' write(*,*) fname(1:ntrim(fname)) nfiles=nfiles-1 !*** decrement -- do not store go to 1 endif else write(*,*) 'first record not version type:' write(*,*) fname(1:ntrim(fname)) nfiles=nfiles-1 !*** decrement -- do not store go to 1 endif afiles(nfiles)=fname !*** store file name asys (nfiles)=asystime !*** store time system fmjdm = -1.d0 !*** to get maximum yyyymmdd for ref date mjdm = -1 call scan1(lin) close(lin) endif go to 1 777 call nummob(isize) write(*,*) 'all obs types: ',isize return end subroutine scan1(lin) *** scan current open rinex file header for info implicit double precision(a-h,o-z) character*80 card character*4 aname character*3 atimsys character*20 antname character*2 aobset(9) logical admitobs,putmob parameter(MXSTA = 1000) !*** max stations dimension dhs(MXSTA),des(MXSTA),dns(MXSTA) character*78 afiles(MXSTA) character*1 asys(MXSTA) character*4 anames(MXSTA) character*20 antnames(MXSTA) save /master/ common/master/dhs,des,dns,nfiles,afiles,asys,anames,antnames,leap save /master2/ common/master2/fmjdm,mjdm,myr,mmn,mdy rewind(lin) 100 read(lin,'(a)') card if(card(61:73).ne.'END OF HEADER') then if (card(61:71).eq.'MARKER NAME') then read(card,'(a4)') aname anames(nfiles)=aname elseif(card(61:72).eq.'LEAP SECONDS') then read(card,'(i6)') leap elseif(card(61:72).eq.'ANT # / TYPE') then read(card,'(20x,a20)') antname antnames(nfiles)=antname elseif(card(61:80).eq.'ANTENNA: DELTA H/E/N') then read(card,'(3f14.4)') dh,de,dn dhs(nfiles)=dh des(nfiles)=de dns(nfiles)=dn elseif(card(61:77).eq.'TIME OF FIRST OBS') then read(card,'(5i6,f13.7,5x,a3)') iy,im,id,ih,imn,sec,atimsys if (atimsys.eq.'GPS') then asys(nfiles)='G' elseif(atimsys.eq.'GLO') then asys(nfiles)='R' elseif(atimsys.eq.'GAL') then asys(nfiles)='E' elseif(atimsys.eq.' ') then continue !*** default to version record endif call civmjd(iy,im,id,ih,imn,sec,mjd,fmjd) if(mjd.gt.mjdm) then fmjdm = fmjd mjdm = mjd myr = iy mmn = im mdy = id elseif(mjd.eq.mjdm) then if(fmjd.gt.fmjdm) then fmjdm = fmjd mjdm = mjd myr = iy mmn = im mdy = id endif endif elseif(card(61:79).eq.'# / TYPES OF OBSERV') then read(card,'(i6)') n 1 if(n.le.9) then read(card,'(6x,9(4x,A2))') (aobset(i),i=1,n) do i=1,n if(admitobs(aobset(i))) then if(.not.putmob(aobset(i),idup)) continue endif enddo else read(card,'(6x,9(4x,A2))') (aobset(i),i=1,9) do i=1,9 if(admitobs(aobset(i))) then if(.not.putmob(aobset(i),idup)) continue endif enddo n=n-9 read(lin,'(a)') card if(card(61:79).ne.'# / TYPES OF OBSERV') then write(*,*) 'FATAL -- need # / TYPES OF OBSERV record' stop 93723 endif go to 1 endif endif go to 100 endif return end logical function admitobs(aobs) *** test for an allowable observation code -- no doppler or sig.strength *** obs family as of rinex 2.12: P,C,L / 1,2,5,6,7,8,A,B,C,D character*2 aobs character*1 a1,a2 a1=aobs(1:1) i=ichar(a1) if((i.ge.97).and.(i.le.122)) a1=char(i-32) !*** case fold to upper admitobs=.false. if(a1.ne.'P'.and.a1.ne.'C'.and.a1.ne.'L') return a2=aobs(2:2) admitobs=.true. if(a2.eq.'1') return if(a2.eq.'2') return if(a2.eq.'5') return if(a2.eq.'6') return if(a2.eq.'7') return if(a2.eq.'8') return i=ichar(a2) if((i.ge.97).and.(i.le.122)) a2=char(i-32) !*** case fold to upper if(a2.eq.'A') return if(a2.eq.'B') return if(a2.eq.'C') return if(a2.eq.'D') return admitobs=.false. return end subroutine getnfiles(nfilesx) *** get number of files implicit double precision(a-h,o-z) parameter(MXSTA = 1000) !*** max stations dimension dhs(MXSTA),des(MXSTA),dns(MXSTA) character*78 afiles(MXSTA) character*1 asys(MXSTA) character*4 anames(MXSTA) character*20 antnames(MXSTA) save /master/ common/master/dhs,des,dns,nfiles,afiles,asys,anames,antnames,leap nfilesx=nfiles return end *----------------------------------------------------------------------- subroutine dump1(lout,itstrt,itstop,itstep) *** output beginning of merged file implicit double precision(a-h,o-z) character*78 fname parameter(MXSTA = 1000) !*** max stations dimension dhs(MXSTA),des(MXSTA),dns(MXSTA) character*78 afiles(MXSTA) character*1 asys(MXSTA) character*4 anames(MXSTA) character*20 antnames(MXSTA) save /master/ common/master/dhs,des,dns,nfiles,afiles,asys,anames,antnames,leap *** 'a' recordanames,antnames write(lout,'(a1,1x,4i6)') 'a',itstrt,itstop,itstep,nfiles *** 'b' records do i=1,nfiles write(lout,'(a1,1x,a4,2x,a20,3(1x,f14.4),2x,a1)') * 'b',anames(i),antnames(i),dhs(i),des(i),dns(i),asys(i) if(asys(i).eq.'R'.and.leap.le.0) then fname=afiles(i) write(*,*) 'WARNING - need leap second record for Glonass' write(*,*) fname(1:ntrim(fname)) endif enddo return end ************************************************************************ subroutine readargs(itstrt,itstop,itstep) *** get command line arguments character*20 aindex *** 1st argument -- start time call getarg(1,aindex) !*** return 1st argument if(aindex.eq.' ') then write(*,*) 'FATAL -- missing command line parm 1' stop 87690 endif read(aindex,*) itstrt if(itstrt.lt.0.or.itstrt.gt.86400) then write(*,*) 'FATAL -- illegal command line parm 1' stop 87691 endif *** 2nd argument -- stop time call getarg(2,aindex) !*** return 2nd argument if(aindex.eq.' ') then write(*,*) 'FATAL -- missing command line parm 2' stop 87680 endif read(aindex,*) itstop if(itstop.lt.0.or.itstop.gt.86400) then write(*,*) 'FATAL -- illegal command line parm 2' stop 87681 endif if(itstop.lt.itstrt) then write(*,*) 'FATAL -- stop time less than start time' stop 87683 endif *** 3rd argument -- time step interval call getarg(3,aindex) !*** return 3rd argument if(aindex.eq.' ') then write(*,*) 'FATAL -- missing command line parm 3' stop 87620 endif read(aindex,*) itstep if(itstep.le.0.or.itstep.gt.86400) then write(*,*) 'FATAL -- illegal command line parm 3' stop 87621 endif write(*,'(a,3i6)') 'start, stop, step:',itstrt,itstop,itstep return end ************************************************************************ *** key management routines ******************************************** ************************************************************************ integer*8 function lkey(kname,ktim,aprn,aobs) *** return integer*8 key for rinex elements implicit integer*8(a-z) integer*4 kname,ktim,keyp,keyo character*3 aprn character*2 aobs if(kname.lt.0) stop 84591 if(ktim.lt.0.or.ktim.gt.86400) stop 84592 i8=kname !*** [1..1000+] j8=ktim !*** [0..86400] k8=keyp(aprn) !*** [0..499] l8=keyo(aobs) !*** [0..39] I00000=int8(100000) I0000000000=I00000*I00000 *** composition: iiiijjjjjkkkll lkey= I0000000000*i8 + * I00000*j8 + * 100*k8 + * l8 return end subroutine invkey(key8,kname,ktim,aprn,aobs) *** return rinex elements from integer*8 key implicit integer*8(a-z) integer*4 kname,ktim,kprn,kobs character*3 aprn character*2 aobs I00000=int8(100000) I0000000000=I00000*I00000 *** decomposition: iiiijjjjjkkkll key=key8 i8 =key/I0000000000 !*** [1..1000+] key=key-I0000000000*i8 j8 =key/I00000 !*** [0..86400] key=key-I00000*j8 k8 =key/100 !*** [0..499] l8 =key-100*k8 !*** [0..39] if(i8.lt.0) stop 84593 if(j8.lt.0.or.j8.gt.86400) stop 84594 if(k8.lt.0.or.k8.gt.499) stop 84595 if(l8.lt.0.or.l8.gt.39) stop 84596 kname = i8 ktim = j8 kprn = k8 kobs = l8 call getaprn(kprn,aprn) call getaobs(kobs,aobs) return end *----------------------------------------------------------------------- integer function keyp(aprn) *** return key [000..499] for alpha prn code *** key 999 reserved for ZZZ code character*3 aprn if(aprn.eq.'ZZZ'.or.aprn.eq.'zzz') then keyp=999 return endif *** prn family as of rinex 2.12: G,R,S,E,C if(aprn(1:1).eq.'G'.or.aprn(1:1).eq.'g'.or.aprn(1:1).eq.' ') then i=0 elseif(aprn(1:1).eq.'R'.or.aprn(1:1).eq.'r') then i=100 elseif(aprn(1:1).eq.'S'.or.aprn(1:1).eq.'s') then i=200 elseif(aprn(1:1).eq.'E'.or.aprn(1:1).eq.'e') then i=300 elseif(aprn(1:1).eq.'C'.or.aprn(1:1).eq.'c') then i=400 else stop 84571 endif read(aprn(2:3),'(i2)') j !*** [00..99] if(j.lt.0) stop 84572 keyp=i+j return end subroutine getaprn(kprn,aprn) *** return alpha prn code from key [000..499] *** key 999 reserved for ZZZ code character*2 a character*3 aprn if(kprn.eq.999) then aprn='ZZZ' !*** reserved code return endif i=kprn/100 j=mod(kprn,100) *** prn family as of rinex 2.12: G,R,S,E,C if (i.eq.0) then aprn(1:1)='G' elseif(i.eq.1) then aprn(1:1)='R' elseif(i.eq.2) then aprn(1:1)='S' elseif(i.eq.3) then aprn(1:1)='E' elseif(i.eq.4) then aprn(1:1)='C' else stop 84573 endif write(a,'(i2.2)') j aprn(2:3)=a return end *----------------------------------------------------------------------- integer function keyo(aobs) *** return key [00..39] for alpha obs code *** key 99 reserved for ZZ code character*2 aobs if(aobs.eq.'ZZ'.or.aobs.eq.'zz') then keyo=99 return endif *** obs family as of rinex 2.12: P,C,L,D / 1,2,5,6,7,8,A,B,C,D if (aobs(1:1).eq.'P'.or.aobs(1:1).eq.'p') then i=0 elseif(aobs(1:1).eq.'C'.or.aobs(1:1).eq.'c') then i=10 elseif(aobs(1:1).eq.'L'.or.aobs(1:1).eq.'l') then i=20 elseif(aobs(1:1).eq.'D'.or.aobs(1:1).eq.'d') then i=30 else stop 84581 endif if (aobs(2:2).eq.'1') then j=0 elseif(aobs(2:2).eq.'2') then j=1 elseif(aobs(2:2).eq.'5') then j=2 elseif(aobs(2:2).eq.'6') then j=3 elseif(aobs(2:2).eq.'7') then j=4 elseif(aobs(2:2).eq.'8') then j=5 elseif(aobs(2:2).eq.'A'.or.aobs(2:2).eq.'a') then j=6 elseif(aobs(2:2).eq.'B'.or.aobs(2:2).eq.'b') then j=7 elseif(aobs(2:2).eq.'C'.or.aobs(2:2).eq.'c') then j=8 elseif(aobs(2:2).eq.'D'.or.aobs(2:2).eq.'d') then j=9 else stop 84582 endif keyo=i+j return end subroutine getaobs(kobs,aobs) *** return alpha obs code from key [00..39] *** key 99 reserved for ZZ code character*2 aobs,aobss(0:39) *** obs family as of rinex 2.12: P,C,L,D / 1,2,5,6,7,8,A,B,C,D data aobss/'P1','P2','P5','P6','P7','P8','PA','PB','PC','PD', * 'C1','C2','C5','C6','C7','C8','CA','CB','CC','CD', * 'L1','L2','L5','L6','L7','L8','LA','LB','LC','LD', * 'D1','D2','D5','D6','D7','D8','DA','DB','DC','DD'/ if(kobs.eq.99) then aobs='ZZ' !*** reserved code return endif if(kobs.lt.0.or.kobs.gt.39) stop 84583 aobs=aobss(kobs) return end ************************************************************************ *** long table ********************************************************* ************************************************************************ subroutine newlng *** initialize the table -- 64 bit integer contents parameter(LENTAB=80000000,NUMKEYS=79999987 ) common/lngtb2/look(NUMKEYS),iptrs(LENTAB),nkey,max,ifree nkey=NUMKEYS !*** nearest prime max =LENTAB !*** length of table ifree=1 do 1 i=1,nkey 1 look(i)=0 return end subroutine numlng(isize) *** return current size of table parameter(LENTAB=80000000,NUMKEYS=79999987 ) common/lngtb2/look(NUMKEYS),iptrs(LENTAB),nkey,max,ifree isize=ifree-1 return end logical function getlng(long,ilong) *** table lookup integer*8 longs,long parameter(LENTAB=80000000,NUMKEYS=79999987 ) common/lngtb1/longs(LENTAB) common/lngtb2/look(NUMKEYS),iptrs(LENTAB),nkey,max,ifree logical loklng if(loklng(long,ilook,iptr)) then ilong=iptr getlng=.true. else ilong=ifree getlng=.false. endif return end logical function putlng(long,idup) *** add entry to table *** idup is location of duplicate entry in list *** idup=0 indicates table overflow integer*8 longs,long parameter(LENTAB=80000000,NUMKEYS=79999987 ) common/lngtb1/longs(LENTAB) common/lngtb2/look(NUMKEYS),iptrs(LENTAB),nkey,max,ifree logical loklng if(loklng(long,ilook,iptr)) then ****** duplicate entry idup=iptr putlng=.false. else if(ifree.gt.max) then ****** table overflow idup=0 putlng=.false. else ****** put entry in table if(iptr.eq.0) then look(ilook)=ifree else iptrs(iptr)=ifree endif longs(ifree)=long iptrs(ifree)=0 ifree=ifree+1 putlng=.true. endif endif return end logical function loklng(long,ilook,iptr) *** table lookup & pointer return integer*8 longs,long,i30,nkey8 parameter(LENTAB=80000000,NUMKEYS=79999987 ) common/lngtb1/longs(LENTAB) common/lngtb2/look(NUMKEYS),iptrs(LENTAB),nkey,max,ifree ****** integer*8 type hash function i30 =1073741824 !*** insure 8 byte (square root-ish) nkey8=nkey !*** insure 8 byte i1 =long/i30 i2 =mod(long,i30) key=mod((i1+i2),nkey)+1 ****** integer*8 type hash function ilook=look(key) if(ilook.eq.0) then ****** immediate failure to find entry ilook=key iptr=0 loklng=.false. else iptr=ilook ilook=0 1 if(long.eq.longs(iptr)) then ****** successful location of entry at iptr loklng=.true. else if(iptrs(iptr).eq.0) then ****** failure to find entry loklng=.false. else ****** chain to next location iptr=iptrs(iptr) go to 1 endif endif endif return end ************************************************************************ *** master sat table *************************************************** ************************************************************************ subroutine newmst *** initialize the table common/msttb2/look(499),iptrs(500),nkey,max,ifree nkey=499 max=500 ifree=1 do 1 i=1,nkey 1 look(i)=0 return end subroutine nummst(isize) *** return current size of table common/msttb2/look(499),iptrs(500),nkey,max,ifree isize=ifree-1 return end logical function getmst(msat,imsat) *** table lookup character*3 msats,msat common/msttb1/msats(500) common/msttb2/look(499),iptrs(500),nkey,max,ifree logical lokmst if(lokmst(msat,ilook,iptr)) then imsat=iptr getmst=.true. else imsat=ifree getmst=.false. endif return end logical function putmst(msat,idup) *** add entry to table *** idup is location of duplicate entry in list *** idup=0 indicates table overflow character*3 msats,msat common/msttb1/msats(500) common/msttb2/look(499),iptrs(500),nkey,max,ifree logical lokmst if(lokmst(msat,ilook,iptr)) then ****** duplicate entry idup=iptr putmst=.false. else if(ifree.gt.max) then ****** table overflow idup=0 putmst=.false. else ****** put entry in table if(iptr.eq.0) then look(ilook)=ifree else iptrs(iptr)=ifree endif msats(ifree)=msat iptrs(ifree)=0 ifree=ifree+1 putmst=.true. endif endif return end logical function lokmst(msat,ilook,iptr) *** table lookup & pointer return character*3 msats,msat common/msttb1/msats(500) common/msttb2/look(499),iptrs(500),nkey,max,ifree ****** character type hash function key=mod(ichar(msat(1:1))* * ichar(msat(2:2))* * ichar(msat(3:3)),nkey)+1 ****** character type hash function ilook=look(key) if(ilook.eq.0) then ****** immediate failure to find entry ilook=key iptr=0 lokmst=.false. else iptr=ilook ilook=0 1 if(msat.eq.msats(iptr)) then ****** successful location of entry at iptr lokmst=.true. else if(iptrs(iptr).eq.0) then ****** failure to find entry lokmst=.false. else ****** chain to next location iptr=iptrs(iptr) go to 1 endif endif endif return end ************************************************************************ *** master obs code table ********************************************** ************************************************************************ subroutine newmob *** initialize the table common/mobtb2/look(47),iptrs(50),nkey,max,ifree nkey=47 max=50 ifree=1 do 1 i=1,nkey 1 look(i)=0 return end subroutine nummob(isize) *** return current size of table common/mobtb2/look(47),iptrs(50),nkey,max,ifree isize=ifree-1 return end logical function getmob(mobs,imobs) *** table lookup character*2 mobss,mobs common/mobtb1/mobss(50) common/mobtb2/look(47),iptrs(50),nkey,max,ifree logical lokmob if(lokmob(mobs,ilook,iptr)) then imobs=iptr getmob=.true. else imobs=ifree getmob=.false. endif return end logical function putmob(mobs,idup) *** add entry to table *** idup is location of duplicate entry in list *** idup=0 indicates table overflow character*2 mobss,mobs common/mobtb1/mobss(50) common/mobtb2/look(47),iptrs(50),nkey,max,ifree logical lokmob if(lokmob(mobs,ilook,iptr)) then ****** duplicate entry idup=iptr putmob=.false. else if(ifree.gt.max) then ****** table overflow idup=0 putmob=.false. else ****** put entry in table if(iptr.eq.0) then look(ilook)=ifree else iptrs(iptr)=ifree endif mobss(ifree)=mobs iptrs(ifree)=0 ifree=ifree+1 putmob=.true. endif endif return end logical function lokmob(mobs,ilook,iptr) *** table lookup & pointer return character*2 mobss,mobs common/mobtb1/mobss(50) common/mobtb2/look(47),iptrs(50),nkey,max,ifree ****** character type hash function key=mod(ichar(mobs(1:1))* * ichar(mobs(2:2)),nkey)+1 ****** character type hash function ilook=look(key) if(ilook.eq.0) then ****** immediate failure to find entry ilook=key iptr=0 lokmob=.false. else iptr=ilook ilook=0 1 if(mobs.eq.mobss(iptr)) then ****** successful location of entry at iptr lokmob=.true. else if(iptrs(iptr).eq.0) then ****** failure to find entry lokmob=.false. else ****** chain to next location iptr=iptrs(iptr) go to 1 endif endif endif return end ************************************************************************ *** time conversion **************************************************** ************************************************************************ subroutine setjd0(iyr,imo,idy) *** set the integer part of a modified julian date as epoch, mjd0 *** the modified julian day is derived from civil time as in civmjd() *** allows single number expression of time in seconds w.r.t. mjd0 implicit double precision(a-h,o-z) integer y save /mjdoff/ common/mjdoff/mjd0 if(iyr.lt.1900) stop 34587 if(imo.le.2) then y=iyr-1 m=imo+12 else y=iyr m=imo endif it1=365.25d0*y it2=30.6001d0*(m+1) mjd=it1+it2+idy-679019 *** now set the epoch for future time computations mjd0=mjd return end subroutine civjts(iyr,imo,idy,ihr,imn,sec,tsec) *** convert civil date to time in seconds past mjd epoch, mjd0 *** requires initialization of mjd0 by setjd0() *** imo in range 1-12, idy in range 1-31 *** only valid in range mar-1900 thru feb-2100 (leap year protocols) *** ref: hofmann-wellenhof, 2nd ed., pg 34-35 *** adapted from civmjd() implicit double precision(a-h,o-z) integer y save /mjdoff/ common/mjdoff/mjd0 if(iyr.lt.1900) stop 34589 if(imo.le.2) then y=iyr-1 m=imo+12 else y=iyr m=imo endif it1=365.25d0*y it2=30.6001d0*(m+1) mjd=it1+it2+idy-679019 tsec=(mjd-mjd0)*86400.d0+3600*ihr+60*imn+sec return end subroutine jtsciv(tsec,iyr,imo,idy,ihr,imn,sec) *** convert time in seconds past mjd0 epoch into civil date *** requires initialization of mjd0 by setjd0() *** imo in range 1-12, idy in range 1-31 *** only valid in range mar-1900 thru feb-2100 *** ref: hofmann-wellenhof, 2nd ed., pg 34-35 *** adapted from mjdciv() implicit double precision(a-h,o-z) save /mjdoff/ common/mjdoff/mjd0 mjd=mjd0+tsec/86400.d0 *** the following equation preserves significant digits fmjd=dmod(tsec,86400.d0)/86400.d0 rjd=mjd+fmjd+2400000.5d0 ia=(rjd+0.5d0) ib=ia+1537 ic=(ib-122.1d0)/365.25d0 id=365.25d0*ic ie=(ib-id)/30.6001d0 *** the fractional part of a julian day is (fractional mjd + 0.5) *** therefore, fractional part of julian day + 0.5 is (fractional mjd) it1=ie*30.6001d0 idy=ib-id-it1+fmjd it2=ie/14.d0 imo=ie-1-12*it2 it3=(7+imo)/10.d0 iyr=ic-4715-it3 tmp=fmjd*24.d0 ihr=tmp tmp=(tmp-ihr)*60.d0 imn=tmp sec=(tmp-imn)*60.d0 return end *----------------------------------------------------------------------- subroutine civmjd(iyr,imo,idy,ihr,imn,sec,mjd,fmjd) *** convert civil date to modified julian date *** imo in range 1-12, idy in range 1-31 *** only valid in range mar-1900 thru feb-2100 (leap year protocols) *** ref: hofmann-wellenhof, 2nd ed., pg 34-35 *** operation confirmed against table 3.3 values on pg.34 implicit double precision(a-h,o-z) integer y if(iyr.lt.1900) stop 34588 if(imo.le.2) then y=iyr-1 m=imo+12 else y=iyr m=imo endif it1=365.25d0*y it2=30.6001d0*(m+1) mjd=it1+it2+idy-679019 fmjd=(3600*ihr+60*imn+sec)/86400.d0 return end subroutine mjdciv(mjd,fmjd,iyr,imo,idy,ihr,imn,sec) *** convert modified julian date to civil date *** imo in range 1-12, idy in range 1-31 *** only valid in range mar-1900 thru feb-2100 *** ref: hofmann-wellenhof, 2nd ed., pg 34-35 *** operation confirmed for leap years (incl. year 2000) implicit double precision(a-h,o-z) rjd=mjd+fmjd+2400000.5d0 ia=(rjd+0.5d0) ib=ia+1537 ic=(ib-122.1d0)/365.25d0 id=365.25d0*ic ie=(ib-id)/30.6001d0 *** the fractional part of a julian day is fractional mjd + 0.5 *** therefore, fractional part of julian day + 0.5 is fractional mjd it1=ie*30.6001d0 idy=ib-id-it1+fmjd it2=ie/14.d0 imo=ie-1-12*it2 it3=(7+imo)/10.d0 iyr=ic-4715-it3 tmp=fmjd*24.d0 ihr=tmp tmp=(tmp-ihr)*60.d0 imn=tmp sec=(tmp-imn)*60.d0 return end ************************************************************************ integer function ntrim(card) *** return nonblank length character*(*) card ntrim=0 if(card.eq.' ') return lc=len(card) do i=lc,1,-1 if(card(i:i).ne.' ') then ntrim=i return endif enddo ntrim=0 return end