program line c c This program demonstrates the use of a set of subroutines c for reading altimeter data from register containing Geosat c ERM data adjusted at the Ohio State University, computing c cross-over differences, estimating bias/tilt parameters, and c correction of the altimeter data. C PERFORMS AN EDITING OF THE TRACKS REMOVING THE ONES WITH c RMS > 2 STD GENERAL STD. c c Programmed by Per Knudsen, KMS - Denmark, Oct 5, 1990. c implicit real*8 (a-h,o-z),integer*4(I-N) parameter(maxdat=150000) character*60 indfil,udfil,indfil2 C C COMMON SPACE OF P.T. 563250*4 BYTES (USED BY SUBROUTINE ABIAS) C common/field/ krev(maxdat),lat(maxdat),lon(maxdat), . kalt(maxdat),kerr(maxdat) LOGICAL LSAVE,LTILT,LQUAD,LNOERR,LNPASS,lplane c read (*,'(A)') indfil read (*,'(A)') indfil2 read (*,'(A)') udfil read(*,*) alat1,alat2,alon1,alon2 open(15,form='formatted',status='old',file=indfil) open(30,form='unformatted',file=udfil) c c C******************************************************** C READ DATA AND STORE IN COMMON/FIELD c IUIN=15 ndat=0 call rdata(iuin,ndat,alat1,alat2,alon1,alon2) write(*,*) 'ndat=',ndat if(ndat.eq.0) goto 80 ndat2=ndat C C************************************************************ C RUNS OALINE TO GET TRACK INFO IUOUT=0 LSAVE=.TRUE. CALL OALINE(IUOUT,NDAT,LSAVE) C************************************************************ C store longest lines first C IUIN=15 ndat=0 call rdata(iuin,ndat,alat1,alat2,alon1,alon2) write(*,*) 'ndat=',ndat if(ndat.eq.0) goto 80 ndat2=ndat c************************************************************ IUOUT=0 LSAVE=.TRUE. CALL ALINE(IUOUT,NDAT,LSAVE,alat1,alat2,alon1,alon2) 80 continue c------------------------------------------- c same for geosat (skipped due to new processing) goto 79 IUIN=16 ndat=0 call rdata(iuin,ndat,alat1,alat2,alon1,alon2) write(*,*) 'ndat=',ndat if(ndat.eq.0) goto 79 ndat2=ndat C C************************************************************ C RUNS OALINE TO GET TRACK INFO IUOUT=0 LSAVE=.TRUE. CALL OALINE(IUOUT,NDAT,LSAVE) C************************************************************ C store longest lines first C IUIN=16 ndat=0 call rdata(iuin,ndat,alat1,alat2,alon1,alon2) write(*,*) 'ndat=',ndat if(ndat.eq.0) goto 79 ndat2=ndat c************************************************************ IUOUT=0 LSAVE=.TRUE. CALL ALINE(IUOUT,NDAT,LSAVE,alat1,alat2,alon1,alon2) 79 continue end C C ******************************************************************* C ***** ***** C ***** SUBROUTINE RDATA C ***** ***** C ******************************************************************* C SUBROUTINE RDATA(IUIN,NDAT,alat1,alat2,alon1,alon2) C ***************************************************************** C C PROGRAMS ARE COPYRIGHT BY THE AUTHOR AND KORT- OG MATRIKELSTYRELSEN. C C SUBROUTINE RDATA READS DATA FROM FILES UNIT NUMBER "IUIN" AND C STORES THE DATA IN THE COMMON ARRAYS. C C PROGRAMMED BY PER KNUDSEN C KORT- OG MATRIKELSTYRELSEN C RENTEMESTERVEJ 8 C DK-2400 KOEBENHAVN NV 15.10.90 C C THIS VERSION: OCT 15, 1990. PER KNUDSEN. C C ***************************************************************** C C CONNECTED FILES: IUIN UNIT NUMBER OF DATA FILE C C CONTENT OF ARRAYS: KREV(.) REVOLUTION NUMBER C LAT(.) LATITUDE (*1000000.) C LON(.) LONGITUDE (*1000000.) C KALT(.) SEA HEIGTH (* 1000.) C KERR(.) STD. DEVIATION (* 1000.) C C ***************************************************************** C IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) parameter(maxdat=150000) common/field/ krev(maxdat),lat(maxdat),lon(maxdat), . kalt(maxdat),kerr(maxdat) c rewind(IUIN) 100 read(iuin,*,end=101) irev,alat,alon,dum,dssh,std irev = irev +1 if(abs(dssh).gt.5.0) goto 100 c if(nnn.lt.5) go to 100 if(alat.lt.alat1.or.alat.gt.alat2) go to 100 if(alon.lt.alon1.or.alon.gt.alon2) go to 100 if((alat.lt.15.15).or.(alat.gt.24.85)) go to 100 if((alon.lt.186.15).or.(alon.gt.192.85)) go to 100 if(std .gt. 0.15d0) go to 100 ndat=ndat+1 krev(ndat)=irev lat(ndat)=alat*1.0d6 lon(ndat)=alon*1.0d6 kalt(ndat)=dssh*1.0d3 kerr(ndat)=std*1.0D3 go to 100 101 continue return end c C ******************************************************************* C ***** ***** C ***** SUBROUTINE ALINE ***** C ***** ***** C ******************************************************************* C SUBROUTINE ALINE(IUOUT,NDAT,LSAVE,alat1,alat2,alon1,alon2) C ***************************************************************** C C PROGRAMMED BY PER KNUDSEN C C ***************************************************************** C C CONNECTED FILES: IUOUT EQUATION SYSTEM FOR EACH TRACK C (TO SUBROUTINE ABIAS) C C CONTENT OF ARRAYS: KREV(.) REVOLUTION NUMBER C LAT(.) LATITUDE (*1000000.) C LON(.) LONGITUDE (*1000000.) C KALT(.) GEOID HEIGTH (* 1000.) C KERR(.) STD. DEVIATION (* 1000.) C LNS(.) = .TRUE. SOUTH-GOING TRACK C = .FALSE. NORTH-GOING TRACK. C C MAX NUMBER OF DATA 100000 C C ***************************************************************** C IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) parameter(maxdat=150000,MAXREV=500000) common/field/ krev(maxdat),lat(maxdat),lon(maxdat), . kalt(maxdat),kerr(maxdat) integer*4 ifirst(maxrev),ilast(maxrev) real*8 rr(maxrev) LOGICAL LSAVE C C ******************************************************************** c open(unit=34,file='fort.34',status='old') open(unit=35,file='fort.35',status='old') rewind(34) rewind(35) c reads info from oaline on tracks to be able to edit tracks. read(35,*) i1,rms if (rms.lt.0.04) rms = 0.04 rrms = 2.0*rms write(*,*) 'from fort 35 ',i1,rrms do i= 1,i1 read(34,*) i3,i2,i3,x1,x1,x1,x1,rms rr(i2) = rms end do c do i=1,ndat c write(36,*) krev(i),lat(i),lon(i) c end do C*** EVALUATION OF NUMBER OF TRACKS AND BEGINNING AND END OF EACH C TRACK. C IREV=1 KREVCC=KREV(1) IFIRST(IREV)=1 DO 40 I=2,NDAT IF(KREV(I).EQ.KREVCC) GO TO 40 ILAST(IREV)=I-1 IREV=IREV+1 IFIRST(IREV)=I KREVCC=KREV(I) 40 CONTINUE ILAST(IREV)=NDAT write(*,*)' ALINE: ndat,ntrack= ',ndat,irev C lat1=alat1*1.0d6 lat2=alat2*1.0d6 lon1=alon1*1.0d6 lon2=alon2*1.0d6 c jrev=0 dum=0.0 do iii=1,irev kkk=iii do jjj=1,irev if((ilast(jjj)-ifirst(jjj)).gt.(ilast(kkk)-ifirst(kkk))) kkk=jjj enddo if((ilast(kkk)-ifirst(kkk)+1).gt.4) then c write(*,*) kkk,krev(ifirst(kkk)),(ilast(kkk)-ifirst(kkk)+1) do 45 i = ifirst(kkk),ilast(kkk) write(27,*) i,rrms,krev(i),lat(i),lon(i),kalt(i),kerr(i) if (lat(i).gt.lat1.and.lat(i).lt.lat2.and. . lon(i).gt.lon1.and.lon(i).lt.lon2) goto 46 45 continue goto 55 46 continue jrev=jrev+1 if (jrev.gt.1598) goto 51 do 50 i=ifirst(kkk),ilast(kkk) alat=lat(i)/1.d6 alon=lon(i)/1.d6 ssh=kalt(i)/1.d3 err=kerr(i)/1.d3 c write(28,500) krev(i),alat,alon,dum,ssh,err c write(28,500) i,alat,alon,dum,ssh,err c if (rr(krev(i)).lt.rrms) write(30) krev(i),lat(i),lon(i),kalt(i),kerr(i) write(29,500) krev(i),alat,alon,dum,ssh,err c write(28,*) rr(krev(i)),rrms,krev(i),lat(i),lon(i),kalt(i),kerr(i) 50 continue 51 continue endif 55 ilast(kkk)=ifirst(kkk)-1 500 format(i6,2f12.5,f4.1,2f9.3) enddo write(*,*)' ALINE: inside region: ntrack= ',jrev return end C ******************************************************************* C ***** ***** C ***** SUBROUTINE ORIGALINE ***** C ***** ***** C ******************************************************************* C SUBROUTINE OALINE(IUOUT,NDAT,LSAVE) C ***************************************************************** C C PROGRAMS ARE COPYRIGHT BY THE AUTHOR AND KORT- OG MATRIKELSTYRELSEN. C C SUBROUTINE ALINE COMPUTES FOR EACH TRACK BIAS/TILT PARAMETERS FROM C RESIDUAL SEA SUEFACE HEIGHTS. THE ELEMENTS IN THE NORMAL MATRIX C (2X2) WHICH IS DIAGONAL USING THE MEAN LONGITUDE AS REFERENCE, AND C THE RIGHT HAND SIDE ARE STORED (BINARY) ON UNIT=IUOUT FOR LATER USE. C C PROGRAMMED BY PER KNUDSEN C KORT- OG MATRIKELSTYRELSEN C RENTEMESTERVEJ 8 C DK-2400 KOEBENHAVN NV 15.10.90 C C THIS VERSION: OCT 15, 1990. PER KNUDSEN. C C ***************************************************************** C C CONNECTED FILES: IUOUT EQUATION SYSTEM FOR EACH TRACK C (TO SUBROUTINE ABIAS) C C CONTENT OF ARRAYS: KREV(.) REVOLUTION NUMBER C LAT(.) LATITUDE (*1000000.) C LON(.) LONGITUDE (*1000000.) C KALT(.) GEOID HEIGTH (* 1000.) C KERR(.) STD. DEVIATION (* 1000.) C LNS(.) = .TRUE. SOUTH-GOING TRACK C = .FALSE. NORTH-GOING TRACK. C C MAX NUMBER OF DATA 100000 C c performs aline but writes parameters to 34 and 35 for editing c of tracks in subsequent process c oa+hf 07-05 C ***************************************************************** C IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) parameter(maxdat=150000,MAXREV=500000) COMMON/field/ krev(maxdat),lat(maxdat),lon(maxdat), . kalt(maxdat),kerr(maxdat) INTEGER*4 IFIRST(MAXREV),ILAST(MAXREV) LOGICAL LSAVE C C ******************************************************************** C open(34) open(35) WRITE(6,490) 490 FORMAT( '1SUBROUTINE ALINE, VERSION 15.10.90,',/) C SUM=0.0D0 SUM2=0.0D0 DO 10 I=1,NDAT DH=KALT(I)/1.D3 SUM = SUM + DH*1.0D0 SUM2= SUM2+(DH*1.0D0)**2 10 CONTINUE IF(NDAT.GT.1) SUM2=DSQRT((SUM2-(SUM**2)/NDAT)/(NDAT-1)) SUM=SUM/NDAT C C*** EVALUATION OF NUMBER OF TRACKS AND BEGINNING AND END OF EACH C TRACK. C IREV=1 KREVCC=KREV(1) IFIRST(IREV)=1 DO 40 I=2,NDAT IF(KREV(I).EQ.KREVCC) GO TO 40 ILAST(IREV)=I-1 IREV=IREV+1 IFIRST(IREV)=I KREVCC=KREV(I) 40 CONTINUE ILAST(IREV)=NDAT write(6,*) SUM,SUM2 C WRITE(6,510) NDAT,IREV,SUM,SUM2 510 FORMAT(' ENTERED WITH ',I6,' DATA FROM ',I4,' TRACKS.', +//, ' MEAN VALUE OF (OBS-REF):',F8.3,' M', +/, ' STD.DEV. - - :',F8.3,' M.') C WRITE(*,520) 520 FORMAT(///, .' I REV NUMBER MEAN.LONG BIAS TILT RMS ', .' RMS',/) C SUM=0.0D0 SUM2=0.0D0 C C** FOR EACH TRACH DO: C jrev = irev DO 200 I=1,IREV NDATTR=ILAST(I)-IFIRST(I)+1 IF(NDATTR.LT.4) jrev = jrev -1 IF(NDATTR.LT.4) GO TO 200 C C FIND MEAN LONGITUDE C AMLON=0.0D0 ERR2=0.0 DO 100 K=IFIRST(I),ILAST(I) ALON=LON(K)/1.0D6 ERR=KERR(K)/1.0D3+0.0001 AMLON=AMLON+ALON ERR2=ERR2+ERR**2 100 CONTINUE AMLON=AMLON/NDATTR ERR2=ERR2/NDATTR C C SET UP NORMAL EQUATION SYSTEM FOR BIAS/TILT ESTIMATION C (THE NORMAL MATRIX IS DIAGONAL WHEN THE MEAN LONGITUDE IS C USED AS REFERENCE) C DIAG1=NDATTR*1.0D0 DIAG2=0.0D0 RIGHT1=0.0D0 RIGHT2=0.0D0 DO 110 K=IFIRST(I),ILAST(I) DLON=LON(K)/1.0D6-AMLON SSH=KALT(K)/1.0D3 DIAG2=DIAG2+DLON**2 RIGHT1=RIGHT1+SSH RIGHT2=RIGHT2+DLON*SSH 110 CONTINUE BIAS=RIGHT1/DIAG1 TILT=RIGHT2/DIAG2 C S=0.0D0 S2=0.0D0 DO 120 K=IFIRST(I),ILAST(I) DLON=LON(K)/1.0D6-AMLON SSH=KALT(K)/1.0D3 HBT=BIAS+DLON*TILT S=S+SSH**2 SSH=SSH-HBT S2=S2+SSH**2 IF(LSAVE) KALT(K)=SSH*1000.0+0.5 SUM=SUM+SSH SUM2=SUM2+SSH**2 120 CONTINUE S=DSQRT(S/NDATTR) S2=DSQRT(S2/NDATTR) DIAG1=DIAG1/ERR2 DIAG2=DIAG2/ERR2 RIGHT1=RIGHT1/ERR2 RIGHT2=RIGHT2/ERR2 if (BIAS.lt.-99) bias = -99 if (BIAS.gt.99) bias = 99 if (TILT.gt.99) TILT = 99 if (TILT.lt.-99) TILT = -99 if (S .gt. 10) S = 10 if (S2 .gt. 10) S2 = 10 if(iuout.gt.6) .WRITE(IUOUT) KREV(IFIRST(I)),AMLON,DIAG1,DIAG2,RIGHT1,RIGHT2 c WRITE(*,530) I,KREV(IFIRST(I)),NDATTR,AMLON,BIAS,TILT,S,S2 if (bias .eq. 99 .or. bias .eq. -99) write(6,*)I,S,S2 WRITE(34,530) I,KREV(IFIRST(I)),NDATTR,AMLON,BIAS,TILT,S,S2 c .ifirst(i),ilast(i) 530 FORMAT(I6,2I8,F10.4,4F8.3,2F12.3) C 200 CONTINUE IF(NDAT.GT.1) SUM2=DSQRT((SUM2-(SUM**2)/NDAT)/(NDAT-1)) SUM=SUM/NDAT WRITE(35,*) JREV,SUM2,NDAT WRITE(6,*) I c if(iuout.gt.6) REWIND(IUOUT) C close(unit=34) close(unit=35) RETURN END