C *************** P R O G R A M A L T G E O S ************************** C 00001**2 C PROGRAM ALTGEOS 00002**2 C 00003**2 C PURPOSE READ GEOS-3 ALTIMETRY SUBSETS (BIO TAPES) , EDIT 00004**2 C THE ALTIMETRY RECORDS , AND (ON OPTION) MERGE THE 00005**2 C DATA WITH PRECISE SATELLITE EPHEMERIS 00006**2 C 00007**2 C I/O DEVICES LUREAD - CARD READER 00008**2 C LULIST - LINE PRINTER 00009**2 C LUBUG - DEBUGGING LOGICAL UNIT (LINE PRINTER) 00010**2 C MTIN - INPUT MAGNETIC TAPE LOGICAL UNIT 00011**2 C MTOUT - OUTPUT MAGNETIC TAPE (EDITED ALTIMETRY) 00012**2 C LOGICAL UNIT 00013**2 C MTAPE1 , MTAPE2 - OUTPUT (TAPE OR DISK) FILE UNITS 00014**2 C 00015**2 C CALLING ROUTINES CSI2 , EQX , DATE , EDIT , PCUBIC , 00016**2 C SPL2 , TIME , CHEBY , ZERO , CLEAR , 00017**2 C XYZPLH , MERGEF 00018**2 C 00019**2 C 00020**2 DOUBLE PRECISION X(289) , XS(32) , YS(32) , YY(32) , F(32) , 00021**2 # G(32) , C(4,32) , CHEB(10) , DCHEB(10) , 00022**2 # CXYZ(10,3) , SC(6) 00023**2 DOUBLE PRECISION FODAY , FDAY , TNODE , XLNODE , STLAT , ENLAT , 00024**2 # STLONG , ENLONG , TB , TS , FT , TK , T , Y , 00025**2 # P , XL , H , YI , PI , RD 00026**2 C 00027**2 DIMENSION IHEAD(18) , IC(6) , IAS(64) , IFORM(650) , IOP(650) 00028**2 DIMENSION DUMMY(2) , IDUMMY(2) 00029**2 C 00030**2 COMMON /EPHEM/ IORD , LOKR , IMIN , CXYZ , IWD1 , IWD2 , ISWT 00031**2 COMMON IEC(5) , IR(5) , JR(15) 00032**2 C 00033**2 DATA DUMMY /0.,-8888./ , IDUMMY /0,-8888/ 00034**2 DATA TB , TS , NS / -53186.56D3 , 10240512.D2 , 32 / 00035**2 DATA LUREAD /5/ , LULIST /6/ , MTIN /10/ , MTOUT /14/ , 00036**2 # MTAPE1 /11/ , MTAPE2 / 12/ , LUBUG /13/ 00037**2 C 00038**2 PI = DARCOS(-1.D0) 00039**2 RD = PI / 180.D0 00040**2 CALL ERRSET (218,256,-1,1,1,1) 00041**2 DO 11 I = 1,5 00042 IR(I) = 0 00043 DO 11 J = I,5 00044 L = (2*5-I)*(I-1)/2 + J 00045 11 JR(L) = 0 00046 C 00047 C READ NUMBER OF FILES TO PROCESS AND HEADER 00048 C 00049 READ(LUREAD,1001) LASTF,NFILE,IHEAD 00050 READ(LUREAD,1003) MERGE , NEFIL , IWRITP 00051**3 1001 FORMAT(2I4,18A4) 00052 WRITE(LULIST,1002) IHEAD 00053 1002 FORMAT('1',30X,18A4,//,1X,'FILE REV YYMMDD START STOP EQXING' 00054 # ,' EQXING STATUS LONG PULSE SHORT PULSE RECORDS DUR START',00055 # ' START STOP STOP',/,17X,3(7H HHMMSS),' LONG STA END'00056 #,2(12H 1 2 3),' STA END',8X,'LAT LON ', 00057 # ' LAT LON',/) 00058 C 00059 C LOOP OVER ALL FILES 00060 C 00061 DO 500 IFILE = 1,NFILE 00062 IFERR = 0 00063 NEDR = 0 00064 C 00065 C LOOP OVER ALL RECORDS 00066 C 00067 DO 400 IREC = 1,650 00068 READ(MTIN,1003,END=410,ERR=900) JREC,JRATE,LREC,(X(I),I=1,LREC)00069 1003 FORMAT(3I6,165F12.0,124F12.0) 00070 GO TO 950 00071 C 00072 900 WRITE(LUBUG,1009) IREC,IFILE 00073 1009 FORMAT (' *** READING ERROR ON RECORD',I6,' OF FILE',I6,' ***') 00074 IFERR = IFERR + 1 00075 IF(IFERR .GT. 10) STOP 00076 JAS = 7 00077 WRITE(LUBUG,1014) IREC,IDUMMY(1),DUMMY(2),(IDUMMY(2),I=1,12), 00078 # (IDUMMY(1),I=1,JAS) 00079 GO TO 400 00080 950 CONTINUE 00081 C OUTPUT RECORD 00082 C 00083 IF(IREC .GT. 1) GO TO 10 00084 C 00085 C PROCESS FOR FIRST RECORD IN FILE 00086 C 00087 DO 5 I = 1,6 00088 5 IC(I) = 0 00089 C 00090 IST = 0 00091 IRC = 0 00092 IFORM (1) = 0 00093 IOP(1) = 0 00094 C 00095 WRITE(LUBUG,1010) IFILE 00096 1010 FORMAT(//,' FILE',I3,/,' REC FM ',7X,'FODAY SMOALT SIG ', 00097 # 3X,'LAT LONG TROP GEOID TIDE SEAH RES AVAGC SIGMA', 00098 # ' H1/3 STATUS') 00099 10 CONTINUE 00100 C 00101 C COMPUTE TELEMETRY FORMAT 00102 C 1, 2, 3 = FORMAT 1,2,3 WITH LONG PULSE MODE 00103 C 11, 12, 13 = FORMAT 1,2,3 WITH SHORT PULSE MODE 00104 C 00105 IFORM(IREC) = X(6) / 10000 + 10 * (X(3) - 40) 00106 J = IFORM(IREC) 00107 IF(J .GT. 3) J = J - 7 00108 IF(J .LT. 0 .OR. J .GT. 6) GO TO 150 00109 IC(J) = IC(J) + 1 00110 C 00111 C STORE CHANGES IN THE ALTIMETER STATUS 00112 C 00113 I1 = 62 00114 I2 = 81 00115 IF(JRATE .EQ. 2) I1 = 86 00116 IF(JRATE .EQ. 2) I2 = 149 00117 N = 1 00118 IAS(1) = X(I1) 00119 I1 = I1 + 1 00120 DO 20 I = I1,I2 00121 IF(X(I) .EQ. X(I-1)) GO TO 20 00122 N = N + 1 00123 IAS(N) = X(I) 00124 20 CONTINUE 00125 C 00126 C FLAG WHETHER ALTIMETER IS OPERATIONAL 00127 IOP(IREC) = 0 00128 IF(N .EQ. 1 .AND. IAS(1) .EQ. 204) IOP(IREC) = 204 00129 IF(N .EQ. 1 .AND. IAS(1) .EQ. 76) IOP(IREC) = 76 00130 IF(N .EQ. 1 .AND. IAS(1) .EQ. 78) IOP(IREC) = 78 00131 IF(N .EQ. 1 .AND. IAS(1) .EQ. 79) IOP(IREC) = 79 00132 C 00133 C SAVE DATA FROM FIRST OPERATIONAL RECORD 00134 C 00135 IF(IOP(IREC) .EQ. 0 .OR. IST .NE. 0) GO TO 25 00136 IST = IOP(IREC) 00137 ISTF = IFORM(IREC) 00138 ISTR = IREC 00139 C 00140 C PROCESS FOR SUMMARY 00141 C 00142 C GET TIME FROM RECORD 00143 IYR = X(20) 00144 NDAY = X(19) 00145 CALL DATE (IYR,NDAY,IMON,IDAY) 00146 FODAY = X(8) 00147 ISTD = IYR * 10000 + IMON * 100 + IDAY 00148 CALL TIME (FODAY,IHR,MIN,ISEC) 00149 ISTT = IHR * 10000 + MIN * 100 + ISEC 00150 C 00151 C IDENTIFY REVOLUTION NUMBER OF PASS 00152 CALL EQX(IYR,NDAY,FODAY,IREV,TNODE,XLNODE) 00153 ID = IDINT(TNODE) 00154 SMUSEC = (TNODE - DFLOAT(ID)) * 86400E6 00155 CALL TIME(SMUSEC,JEHR,JEM,JES) 00156 IEQXT = JEHR * 10000 + JEM * 100 + JES 00157 C 00158 C GET LATITUDE AND LONGITUDE FROM RECORD 00159 STLAT = X(10) / 1.D4 00160 STLONG = X(11) / 1.D4 00161 C 00162**2 C 00163**2 C PREPARE FOR PRECISE EPHEMERIS MERGING 00164**2 C 00165**2 LOKG = NDAY*1440 + IHR*60 + MIN 00166**2 IF(MERGE.EQ.0) GO TO 15 00167**2 C 00168**2 CALL MERGEF (LOKG,NEFIL,IMCH) 00169**2 GO TO (16,15) , IMCH 00170**2 16 TK = 2.D0 / (IMIN-1) 00171**2 WRITE(LUBUG,603) IREV 00172**2 603 FORMAT(//,10X,'*** PRECISE EPHEMERIS UPDATE FOUND FOR REV',I5, 00173**2 # ' ***',//) 00174**2 GO TO 25 00175**2 15 WRITE(LUBUG,602) IREV 00176**2 602 FORMAT(//,'*** NO UPDATED EPHEMERIS FOR REV ',I5,' ***',/ ) 00177**2 25 CONTINUE 00178 C 00179 C SAVE DATA FROM LAST OPERATIONAL RECORD 00180 C 00181 IF(IOP(IREC) .EQ. 0 .AND. IST .NE. 0) GO TO 30 00182 IEN = IOP(IREC) 00183 IENF = IFORM(IREC) 00184 IENR = IREC 00185 CALL TIME(X(8),IEHR,IEMIN,IESEC) 00186 IENT = IEHR * 10000 + IEMIN * 100 + IESEC 00187 ENLAT = X(10) / 1.D4 00188 ENLONG = X(11) / 1.D4 00189 NUREC = IENR - ISTR + 1 00190 IDUR = IDINT(X(8) - FODAY) / 1E6 00191 30 CONTINUE 00192 C 00193 C LIST ONE LINE PER RECORD SUMMARY 00194 C 00195 35 IRC = IREC 00196 ISM = 165 00197 IF(JRATE .EQ. 2) ISM = 289 00198 FDAY = X(8) 00199 LSMALT = X(9) 00200 LSIG = X(12) 00201 LAT = X(10) 00202 LONG = X(11) 00203 LTROP = X(14) 00204 LGEOID = X(16) 00205 LTIDE = X(17) 00206 MSEAH = X(ISM) 00207 C 00208 CALL CLEAR (XS,NS) 00209 CALL CLEAR (YS,NS) 00210 CALL CLEAR (F ,NS) 00211 CALL CLEAR (G ,NS) 00212 CALL ZERO (C,4,NS) 00213 C 00214**2 IF(IOP(IREC).EQ.0) GO TO 41 00215**2 C 00216**2 C 00217**2 C IF MERGE.NE.0 MERGE DATA WITH PRECISE SATELLITE EPHEMERIS 00218**2 C 00219**2 IF(MERGE.EQ.0 .OR. IMCH.EQ.2) GO TO 21 00220**2 C 00221**2 C INTERPOLATE SATELLITE X , Y , Z COORDINATES 00222**2 C AT CURRENT EPOCH WITHIN FRAME 00223**2 C 00224**2 CALL TIME (FDAY,JHR,JMIN,JSEC) 00225**2 T = NDAY*1440.D0 + FDAY/6.D7 00226**2 Y = (T-DFLOAT(LOKR))*TK - 1.D0 00227**2 C 00228**2 CALL CHEBY (CHEB,DCHEB,2,IORD,Y,TK) 00229**2 C 00230**2 DO 18 J = 1,3 00231**2 SC(J) = CXYZ(1,J)/1.D3 00232**2 DO 18 K = 2,IORD 00233**2 18 SC(J) = SC(J) + CHEB(K)*CXYZ(K,J)/1.D3 00234**2 C 00235**2 DO 19 J = 1,3 00236**2 19 SC(J) = SC(J)*1.D3 00237**2 C 00238**2 C 00239**2 C CONVERT X , Y , Z SATELLITE COORDINATES INTO PHI, LAMDA , H 00240**2 C AND STORE INTO WORDS #13 , 15 , 5 (CURRENTLY ZERO) OF RECORD 00241**2 C 00242**2 C 00243**2 C COMPUTED SATELLITE HEIGHT IS BASED ON REFERENCE ELLIPSOID DEFINED 00244**2 C BY SEMI-MAJOR AXIS A=6378145. AND FLATTENING F=1/298.255 00245**2 C 00246**2 CALL XYZPLH (SC(1),SC(2),SC(3),P,XL,H) 00247**2 C 00248**2 P = (P/RD)*1.D4 00249**2 IF(XL.LT.0.D0) XL = (360.D0+(XL/RD)) * 1.D4 00250**2 C 00251**2 X(5) = H * 1.D2 00252**2 X(13) = P 00253**2 X(15) = XL 00254**2 C 00255**2 LSMALT = X(5) 00256**2 LAT = X(13) 00257**2 LONG = X(15) 00258**2 C 00259**2 21 CONTINUE 00260**2 C 00261**2 C 00262 C 00263 C COMPUTE A SMOOTH SEA SURFACE HEIGHT PER MAJOR 00264 C FRAME USING A CUBIC SPLINE FIT. 00265 C 00266 GO TO (37,38) , JRATE 00267 37 ISSH = 142 00268 JSSH = ISM - 4 00269 IF(IMCH.EQ.1) ISSH = 21 00270**2 IF(IMCH.EQ.1) JSSH =40 00271**2 GO TO 39 00272 C 00273 38 ISSH = 254 00274 JSSH = ISM - 4 00275 IF(IMCH.EQ.1) ISSH = 21 00276**2 IF(IMCH.EQ.1) JSSH = 52 00277**2 39 L = 0 00278 FT = X(18) * 1.D4 00279 DO 40 K1 = ISSH,JSSH 00280 C 00281 C DISREGARD DATA FLAGGED FOR REJECTION 00282 C 00283**2 IF(IMCH.EQ.1 .AND. X(K1).LT.0.D0) GO TO 40 00284**2 IF(IMCH.EQ.2 .AND. DABS(X(K1)).GT.88.D2) GO TO 40 00285**2 C 00286 L = L + 1 00287 KM1 = L - 1 00288 YS(L) = X(K1) 00289 C 00290 C TIME TAG SAMPLE DATA USING THE ALGORITHM GIVEN IN NASA'S 00291**2 C PREPROCESSING REPORT (LEITAO ET AL , 1975 , P.54) 00292**2 C 00293**2 C TB - TIME WITHIN FRAME 00294**2 C FT - LOW DATA RATE TIME FROM FRAME TIME 00295**2 C 00296**2 C 00297 XS(L) = (FT+TB+DFLOAT(KM1)*TS) / 1.D4 00298 C 00299 40 CONTINUE 00300 C 00301 IF(L.LT.3) GO TO 41 00302 YY(1) = 0.0D0 00303 YY(L) = 0.0D0 00304 C 00305 CALL SPL2 (XS,YS,YY,L,F,G) 00306 CALL CSI2 ( XS,YS,YY,L,C) 00307 YI = PCUBIC (FDAY,XS,L,C) 00308 LSEAH = YI 00309 IF(IMCH.EQ.1) LSEAH = H*1.D2 - YI + LTROP 00310**2 GO TO 42 00311 C 00312 41 LSEAH = -8888 00313 42 CONTINUE 00314 C 00315 IRES = X(16) + X(17) - X(ISM) 00316 LAVAGC = X(ISM-2) 00317 LH = X(ISM-1) 00318 C 00319 I1 = 102 00320 I2 = 121 00321 IF(JRATE .EQ. 2) I1 = 182 00322 IF(JRATE .EQ. 2) I2 = 213 00323 NM1 = I2 - I1 00324 VV = 0. 00325 DO 36 I = I1,I2 00326 36 VV = VV + ( X(ISM-2) - X(I) )**2 00327 LSAGC = SQRT(VV/NM1) 00328 1014 FORMAT(1X,I4,2X,I2,1X,F12.0,I9,I6,2I8,8I6,7I4,9(/,101X,7I4)) 00329 C 00330 C EDIT RECORD 00331 C 00332 CALL EDIT(LSIG,LSEAH,LAVAGC,LSAGC,LH,LREJ) 00333 C 00334 IF(LREJ .EQ. 1) GO TO 400 00335 NEDR = NEDR + 1 00336 C 00337 C SAVE DATA FROM FIRST EDITED RECORD 00338 C 00339 IF(NEDR .GT. 1) GO TO 411 00340 SLAT = LAT / 1.E4 00341 SLONG = LONG / 1.E4 00342 CALL TIME(FDAY,LHOUR,LMIN,LSEC) 00343 ISTIME = LHOUR*10000 + LMIN*100 + LSEC 00344 GO TO 413 00345 C 00346 C SAVE DATA FROM LAST EDITED RECORD 00347 C 00348 411 ELAT = LAT / 1.E4 00349 ELONG = LONG / 1.E4 00350 CALL TIME(FDAY,NHOUR,NMIN,NSEC) 00351 IETIME = NHOUR*10000 + NMIN*100 + NSEC 00352 C 00353 C OUTPUT EDITED RECORD 00354 C 00355 413 IF(N .GT. 7) N = 7 00356 412 WRITE(MTAPE1,1024) NEDR,IREC,IFORM(IREC),FDAY,LSMALT,LSIG,LAT, 00357 # LONG,LTROP,LGEOID,LTIDE,LSEAH,IRES,LAVAGC,LSAGC,LH,MSEAH,N, 00358 # (IAS(I),I=1,N) 00359 1024 FORMAT(2I5,I4,F13.0,I9,I6,2I8,10I6,7I4) 00360 C 00361 WRITE(LUBUG,1014) IREC,IFORM(IREC),FDAY,LSMALT,LSIG,LAT,LONG, 00362**2 # LTROP,LGEOID,LTIDE,LSEAH,IRES,LAVAGC,LSAGC,LH,(IAS(I),I=1,N) 00363**3 C 00364**2 C 00365**3 C REPLACE THE LAST WORD OF RECORD (CURRENTLY NOT USED) WITH 00366**3 C THE REVOLUTION NUMBER OF THE CURRENT PASS 00367**3 C 00368**3 IF(IWRITP.EQ.0) GO TO 414 00369**3 IF(JRATE.EQ.1) X(165) = DFLOAT(IREV) 00370**3 IF(JRATE.EQ.2) X(289) = DFLOAT(IREV) 00371**3 C 00372**3 WRITE(MTOUT,1003) JREC,JRATE,LREC,(X(I),I=1,LREC) 00373**3 C 00374**3 414 CONTINUE 00375**3 400 CONTINUE 00376 C 00377 C 00378 C LIST ONE LINE PER FILE SUMMARY 00379 C 00380 410 WRITE(LULIST,1013) IFILE,IREV,ISTD,ISTT,IENT,IEQXT,XLNODE, 00381 # IST,IEN,IC,ISTR,IENR,IDUR,STLAT,STLONG,ENLAT,ENLONG 00382 1013 FORMAT(2I5,4I7,F7.2,11I4,4F9.4) 00383 C 00384 END FILE MTAPE1 00385 IF(IWRITP.EQ.1) END FILE MTOUT 00386**3 C 00387 LASTF = LASTF + 1 00388 WRITE(MTAPE2,1023) LASTF,IREV,ISTD,ISTIME,IETIME,IEQXT,XLNODE, 00389 # SLAT,SLONG,ELAT,ELONG 00390 1023 FORMAT(6I7,5F9.4) 00391 500 CONTINUE 00392 C 00393 END FILE MTAPE2 00394 C 00395 WRITE(LULIST,1029) 00396 1029 FORMAT('1',30X,'*** DATA REJECTION ***',//) 00397 DO 12 I = 1,5 00398 12 WRITE(LULIST,1030) IR(I),( JR((2*5-J)*(J-1)/2 + I),J=1,I) 00399 1030 FORMAT(22X,I6,4X,5I6,/) 00400 STOP 00401 150 WRITE(LULIST,1015) J 00402 1015 FORMAT(5X,'INCONSISTENT INPUT TAPE FORMAT',4X,'FORMAT=',I4) 00403 STOP 00404 END 00405 SUBROUTINE DATE (IYR,NDAY,IMON,IDAY) C 00407**2 C CONVERT DAY NUMBER IN YEAR INTO MONTH AND DAY OF THE MONTH 00408**2 C 00409**2 C PARAMETERS IYR - THE LAST TWO DIGITS OF THE YEAR 00410**2 C NDAY - DAY NUMBER IN YEAR IYR 00411**2 C IMON - MONTH IN IYR 00412**2 C IDAY - DAY IN IMON 00413**2 C 00414**2 DIMENSION MON(12) 00415 DATA MON /31,28,31,30,31,30,31,31,30,31,30,31/ 00416 MON(2) = 28 00417 IF(MOD(IYR,4) .EQ. 0) MON(2) = 29 00418 IDAY = NDAY 00419 DO 10 I = 1,12 00420 IMON = I 00421 IDAY = IDAY - MON(I) 00422 IF(IDAY .LE. 0) GO TO 20 00423 10 CONTINUE 00424 20 IDAY = IDAY + MON(IMON) 00425 RETURN 00426 END 00427 SUBROUTINE EDIT(LSIG,LSEAH,LAVAGC,LSAGC,LH,LREJ) C 00429**2 C SUBROUTINE EDIT 00430**2 C 00431**2 C PURPOSE EDIT FRAME MEASUREMENTS AND DELETE NOISY DATA 00432**2 C 00433**2 C PARAMETERS LSIG - ALTITUDE MEASUREMENT STANDARD DEVIATION (CM)00434**2 C LSEAH - SMOOTH SEA SURFACE HEIGHT (CM) 00435**2 C LAVAGS - RADAR ALTIMETER AUTOMATIC GAIN (AGCV) 00436**2 C CONTROL VOLTAGE (DBM * 1.E2) 00437**2 C LH - THE H 1/3 VALUE (SEA STATE MEASURE) 00438**2 C 00439**2 C REFERENCES RAPP , R.H. (1978) - "MEAN GRAVITY ANOMALIES AND 00440**2 C SEA SURFACE HEIGHTS DERIVED FROM GEOS-3 ALTIMETRY 00441**2 C DATA" , NASA CR 156843. 00442**2 C 00443**2 C 00444**2 C COMMENTS 00445**2 C IR , JR - ARE VECTORS RELATED TO THE NUMBER OF 00446**2 C MEASUREMENTS DELETED BY ONE OR MORE OF THE 00447**2 C THE APPLIED EDIT CRITERIA 00448**2 C 00449**2 C IR(I) INDICATES THE NUMBER OF MEASUREMENTS DELETED BECAUSE 00450**2 C THE I-TH CRITERION ALONE WAS MET 00451**2 C 00452**2 C JR IS THE VECTOR TRANSCRIPTION OF AN ARRAY WHOSE I-TH DIAGONAL 00453**2 C ELEMENTS INDICATE THE NUMBER OF MEASUREMENTS DELETED BECAUSE00454**3 C THE I-TH CRITERION AND POSSIBLY ONE OR MORE OTHER CRITERIA 00455**2 C WERE MET, AND THE OFF DIAGONAL ELEMENTS (I,J) INDICATE 00456**3 C THE NUMBER OF MEASUREMENTS DELETED BECAUSE THE I-TH AND 00457**2 C THE J-TH CRITERIA WERE MET TOGETHER 00458**2 C 00459**2 COMMON IEC(5),IR(5),JR(15) 00460 C 00461 DO 100 I = 1,5 00462 100 IEC(I) = 0 00463 C 00464 IS = 0 00465 LREJ = 0 00466 C 00467 C 00468**2 C A FRAME MEASUREMENT IS ACCEPTED IF THE FOLLOWING 5 CRITERIA ARE MET 00469**2 C 00470**2 IF(LSIG .GT. 25 .AND. LSIG .LT.8E2) GO TO 10 00471 IEC(1) = 1 00472 10 IF(IABS(LSEAH) .LT. 17E3) GO TO 20 00473 IEC(2) = 2 00474 20 IF(LAVAGC.GT.(-78E2) .AND. LAVAGC.LT.(-62E2)) GO TO 30 00475 IEC(3) = 3 00476 30 IF(LSAGC .LT. 2E2) GO TO 40 00477 IEC(4) = 4 00478 40 IF( IABS(LH) .LT. 88E2 ) GO TO 50 00479 IEC(5) = 5 00480 C 00481 50 DO 200 I = 1,5 00482 200 IS = IS + IEC(I) 00483 IF(IS .EQ. 0) GO TO 90 00484 GO TO (1,2,3,4,5) , IS 00485 60 GO TO 70 00486 C 00487 1 IR(1) = IR(1) + 1 00488 GO TO 80 00489 2 IR(2) = IR(2) + 1 00490 GO TO 80 00491 3 IF(IEC(3) .EQ. 0) GO TO 70 00492 IR(3) = IR(3) + 1 00493 GO TO 80 00494 4 IF(IEC(4) .EQ. 0) GO TO 70 00495 IR(4) = IR(4) + 1 00496 GO TO 80 00497 5 IF(IEC(5) .EQ. 0) GO TO 70 00498 IR(5) = IR(5) + 1 00499 GO TO 80 00500 C 00501 70 CONTINUE 00502 DO 300 I = 1,5 00503 IF( IEC(I) .EQ. 0) GO TO 300 00504 J = I 00505 L = (10 - I)*(I-1)/2 + J 00506 JR(L) = JR(L) + 1 00507 IF(I .EQ. 5) GO TO 300 00508 J = I + 1 00509 DO 301 K = J,5 00510 IF( IEC(K) .EQ. 0 ) GO TO 301 00511 L = (10-I)*(I-1)/2 + K 00512 JR(L) = JR(L) + 1 00513 301 CONTINUE 00514 300 CONTINUE 00515 C 00516 80 LREJ = 1 00517 90 RETURN 00518 END 00519 SUBROUTINE EQX(IYR,NDAY,FODAY,IREV,TNODE,XLNODE) C 00521**2 C SUBROUTINE EQX 00522**2 C 00523**2 C PURPOSE GIVEN TIME TO COMPUTE REVOLUTION NUMBER AND 00524**2 C NODAL CROSSING FOR THE GEOS-3 SATELLITE 00525**2 C 00526**2 C PARAMETERS IYR - LAST TWO DIGITS OF YEAR 00527**2 C NDAY - DAY NUMBER 00528**2 C FODAY - TIME OF DAY IN MICROSECONDS 00529**2 C IREV - REVOLUTION NUMBER 00530**2 C TNODE - TIME OF NODAL CROSSING 00531**2 C XLNODE - LONGITUDE OF NODAL CROSSING 00532**2 C 00533**2 DOUBLE PRECISION FODAY,TNODE,XLNODE,P,PD,W, 00534 # WD,TN,XN,PER,DAY,DT 00535 DATA P /101.78528D0/, 00536 # PD /-26.D-8/, 00537 # W /25.32335D0/, 00538 # WD /-6.D-8/, 00539 # KO /3200/, 00540 # TN /326.160439815D0/, 00541 # XN /17.11D0/ 00542 PER = P + PD * KO 00543 DAY = NDAY + FODAY / 86400.D6 00544 DT = DAY - TN 00545 IF(IYR .EQ. 76) DT = DT + 365. 00546 IF(IYR .EQ. 77) DT = DT + 731. 00547 NREV = (DT / PER) * 1440.D0 00548 N = IABS(NREV) 00549 IF(NREV .LT. 0) N = N + 1 00550 NSIGN = ISIGN(1,NREV) 00551 IREV = KO 00552 TNODE = TN 00553 XLNODE = XN 00554 DO 10 I = 1,N 00555 IREV = IREV + NSIGN 00556 PER = P + PD * IREV 00557 WM = W + WD * IREV 00558 TNODE = TNODE + PER * NSIGN / 1440.D0 00559 XLNODE = XLNODE - WM * NSIGN 00560 IF(XLNODE .LT. -180.D0) XLNODE = XLNODE + 360.D0 00561 IF(XLNODE .GT. 180 .D0) XLNODE = XLNODE - 360.D0 00562 10 CONTINUE 00563 RETURN 00564 END 00565 SUBROUTINE TIME(FODAY,IHR,MIN,ISEC) C 00567**2 C CONVERT TIME GIVEN IN MICROSECONDS INTO 00568**2 C HOURS , MINUTES AND SECONDS 00569**2 C 00570**2 DOUBLE PRECISION FODAY 00571 SEC = FODAY / 1E6 00572 IHR = SEC / 3600. 00573 MIN = AMOD(SEC,3600.) / 60. 00574 ISEC = AMOD(SEC,60.) 00575 RETURN 00576 END 00577 SUBROUTINE SPL2 ( X,Y,Y2,N,F,G ) C 00579 C GIVEN A SET OF FUNCTIONAL VALUES (Y) CORRESPONDING TO A SET C OF POINTS (X) TOGETHER WITH THE SECOND DERIVATIVES AT THE C END POINTS Y2(1) , Y2(N) SUBROUTINE SPL2 COMPUTES THE C SECOND DERIVATIVES AT ALL OTHER POINTS C 00584 C DESCRIPTION OF PARAMETERS 00585 C X - VECTOR OF ARGUMENTS 00586 C Y - ASSOCIATED VECTOR OF FUNCTIONAL VALUES 00587 C Y2 - VECTOR OF SECOND DERIVATIVES 00588 C N - NUMBER OF MESH POINTS 00589 C F,G - AUXILIARY VECTORS 00590 C 00591 IMPLICIT REAL*8(A-H,O-Z) 00592 DIMENSION X(N),Y(N),Y2(N),F(N),G(N) 00593 C 00594 F(1) = .0 00595 G(1) = .0 00596 NM1 = N - 1 00597 DO 2 K = 1,NM1 00598 H2 = X(K+1) - X(K) 00599 R2 = (Y(K+1) - Y(K))/H2 00600 IF(K.EQ.1) GO TO 1 00601 Z = 1./(2.*(H1+H2) - H1*G(J)) 00602 G(K) = Z*H2 00603 H = 6.*(R2-R1) 00604 IF(K.EQ.2) H = H - H1*Y2(1) 00605 IF(K.EQ.NM1) H = H - H2*Y2(N) 00606 F(K) = Z*(H-H1*F(J)) 00607 1 J = K 00608 H1 = H2 00609 R1 = R2 00610 2 CONTINUE 00611 Y2(NM1) = F(NM1) 00612 IF(NM1.LE.2) RETURN 00613 NJ = NM1 - 1 00614 DO 3 J1 = 2,NJ 00615 K = N-J1 00616 Y2(K) = F(K) - G(K)*Y2(K+1) 00617 3 CONTINUE 00618 RETURN 00619 END 00620 SUBROUTINE CSI2 (X,Y,Y2,N,C ) C 00622 C SUBROUTINE CSI2 COMPUTES THE POLYNOMIAL COEFFICIENTS OF EACH 00623 C DEGREE POLYNOMIAL PIECE OF THE CUBIC SPLINE AT EACH INTERVAL C BETWEEN POINTS I AND (I+1) UTILIZING THE COMPUTED MOMENTS AT C THESE POINTS C 00626 C DESCRIPTION OF PARAMETERS 00627 C X - VECTOR OF ARGUMENTS 00628 C Y - ASSOCIATED VECTOR OF FUNCTIONAL VALUES 00629 C Y2 - VECTOR OF SECOND DERIVATIVES 00630 C N - NUMBER OF MESH POINTS 00631 C C - MATRIX OF COMPUTED COEFFICIENTS 00632 C 00633 IMPLICIT REAL*8(A-H,O-Z) 00634 DIMENSION X(N),Y(N),Y2(N),C(4,N) 00635 C 00636 NM1 = N - 1 00637 DO 1 I = 1,NM1 00638 DX = X(I+1)- X(I) 00639 DY = Y(I+1)- Y(I) 00640 C(1,I) = Y(I) 00641 C(2,I) = (DY/DX) - (DX*(Y2(I+1) + 2.*Y2(I))/6. ) 00642 C(3,I) = Y2(I)/2. 00643 1 C(4,I) = (Y2(I+1) - Y2(I))/(6.*DX) 00644 RETURN 00645 END 00646 FUNCTION PCUBIC( XI,X,N,C ) C 00648 C SUBROUTINE FUNCTION PCUBIC EVALUATES THE FINAL INTERPOLATION C VALUE OF THE COMPUTED PIECEWISE CUBIC FUNCTION FOR C A GIVEN ARGUMENT 'XI' C 00651 IMPLICIT REAL*8(A-H,O-Z) 00652 DIMENSION X(N) , C(4,N) 00653 I = 1 00654 DX = XI - X(I) 00655 IF(DX) 1,3,2 00656 1 IF(I.EQ.1) GO TO 3 00657 I = I - 1 00658 DX = XI - X(I) 00659 IF(DX) 1,3,3 00660 4 I = I + 1 00661 DX = DDX 00662 2 IF(I.EQ.N) GO TO 3 00663 DDX = XI - X(I+1) 00664 IF(DDX) 3,4,4 00665 3 PCUBIC = C(1,I) + DX*( C(2,I) + DX*( C(3,I) + DX*C(4,I))) 00666 RETURN 00667 END 00668 SUBROUTINE CLEAR (AA,N) C 00670**2 C 00671**2 C SUBROUTINE CLEAR 00672**2 C 00673**2 C PURPOSE TO CLEAR A VECTOR OF LENGTH 'N' 00674**2 C 00675**2 DIMENSION AA(1) 00676 DO 1 M = 1,N 00677 1 AA(M) = 0.0 00678 RETURN 00679 END 00680 SUBROUTINE ZERO (A,N,M) C 00682**2 C SUBROUTINE ZERO 00683**2 C 00684**2 C PURPOSE TO CLEAR AN ARRAY OF SIZE N BY M 00685**2 C 00686**2 DIMENSION A(1) 00687 I = N * M 00688 DO 1 J = 1,I 00689 1 A(J) = 0.0 00690 RETURN 00691 END 00692 SUBROUTINE MERGEF (LOKG,NEFIL,IMCH) C 00694**2 C SUBROUTINE MERGEF 00695**2 C 00696**2 C PURPOSE TO MERGE ALTIMETRY DATA WITH PRECISE SATELLITE 00697**2 C EPHEMERIS (FITTED NWL PRECISE EPHEMERIS) 00698**2 C 00699**2 C I/O DEVICES IEPH - INPUT MAGNETIC TAPE (SATELLITE EPHEMERIS) 00700**2 C 00701**2 C CALLING ROUTINES READEF 00702**2 C 00703**2 C PARAMETERS LOKG - LOCKON TIME (MINUTES) 00704**2 C NEFIL - NUMBER OF EPHEMERIS FILES TO SEARCH FOR 00705**2 C EPHEMERIS MATCH TO CURRENT SATELLITE PASS 00706**2 C IMCH - EPHEMERIS MATCH CODE 00707**2 C 1 - SUCCESSFUL RETURN 00708**2 C 2 - UNSUCCESSFUL RETURN 00709**2 C 00710**2 IMPLICIT REAL*8(A-H,O-Z) 00711**2 COMMON/ EPHEM/ IORD,LOKR,IMIN,CXYZ,IWD1,IWD2,ISWT 00712**2 DIMENSION ID(7),CXYZ(10,3) 00713**2 DATA IEPH/ 4/ 00714**2 NEFL = NEFIL + 1 00715**2 N = 0 00716**2 35 N = N + 1 00717**2 IF(N.EQ.NEFL) GO TO 50 00718**2 READ(IEPH,END = 60) IYER,ISAT,IDAY,ID,IORD 00719**2 40 CALL READEF 00720**2 GO TO (30,35) , ISWT 00721**2 30 LOKEND = LOKR + (IMIN-1) 00722**2 IF(LOKR .GT. LOKG) GO TO 50 00723**2 IF(LOKG.GT.LOKR .AND. LOKG.LT.LOKEND) GO TO 55 00724**2 GO TO 40 00725**2 55 CONTINUE 00726**2 IMCH = 1 00727**2 REWIND IEPH 00728**2 RETURN 00729**2 50 IMCH = 2 00730**2 REWIND IEPH 00731**2 60 RETURN 00732**2 END 00733**2 SUBROUTINE READEF C 00735**2 C SUBROUTINE READEF 00736**2 C 00737**2 C PURPOSE READ FITTED SATELLITE PRECISE EPHEMERIS 00738**2 C 00739**2 C I/O DEVICES IEPH - INPUT MAGNETIC TAPE (SATELLITE EPHEMERIS) 00740**2 C 00741**2 C PARAMETERS IORD - ORDER OF TCHEBYCHEV POLYNOMIAL FITTING 00742**2 C OF THE STANDARD NWL PRECISE EPHEMERIS 00743**2 C (CURRENT MAXIMUM LIMIT IORD=9) 00744**2 C LOKR - REFERENCE LOCKON TIME (MINUTES) 00745**2 C IMIN - TIME INTERVAL (FROM LOKR) OF PRECISE 00746**2 C EPHEMERIS FOR CURRENT PASS 00747**2 C CHEB - TCHEBYCHEV POLYNOMIAL COEFFICIENTS FOR 00748**2 C POSITION COMPONENTS 00749**2 C DCHEB - TCHEBYCHEV POLYNOMIAL COEFFICIENTS FOR 00750**2 C VELOCITY COMPONENTS 00751**2 C ISWT - END OF FILE INDICATOR (1 - NO , 2 - YES) 00752**2 C IWD1 , IWD2 - NOT USED 00753**2 C 00754**2 IMPLICIT REAL*8(A-H,O-Z) 00755**2 COMMON/ EPHEM/ IORD,LOKR,IMIN,CXYZ,IWD1,IWD2,ISWT 00756**2 DIMENSION CXYZ(10,3) 00757**2 DATA IEPH /4/ 00758**2 ISWT = 1 00759**2 READ(IEPH,END=10) LOKR,IMIN,((CXYZ(J,I),J=1,IORD),I=1,3),IWD1,IWD200760**2 RETURN 00761**2 10 ISWT = 2 00762**2 RETURN 00763**2 END 00764**2