C * * * * * S E A H T * * * * * C C C 00004**2 C TYPE PROGRAM 00005**2 C 00006**2 C PURPOSE LEAST SQUARES ADJUSTMENT OF INTERSECTING ARCS FOR 00007**2 C SEA SURFACE REFINEMENT 00008**2 C 00009**2 C CALLING ROUTINES TIDE , ELAPST , CARTE , PREDCT , YRMDAY , 00010**2 C HHMMDD , ZERO , CLEAR , DELETE , SELECT , 00011**2 C MATMPY , SOLVEQ , SPIN , TAURE , SORT , 00012**2 C LPRNT , PBLANK , PLOT 00013**2 C 00014**2 C I/O DEVICES 00015**2 C 00016**2 C LOGICAL UNIT USE 00017**2 C 00018**2 C INSSH - 1 INPUT SSH DATA 00019**2 C INTIDE - 2 INPUT DATA FOR TIDAL MODEL 00020**2 C ICDTA - 3 OUTPUT DISK FILE OF CORRECTED DATA 00021**2 C ICR - 5 CARD INPUT 00022**2 C IPRT - 6 PRINT OUTPUT 00023**2 C IPCH - 7 PUNCH OUTPUT 00024**2 C IOUT - 8 PRINT OUTPUT 00025**2 C ITAPE - 9 OUTPUT CLEAN TAPE 00026**3 C 00027**2 C ITF1 - ITF4 TEMPORARY DISK STORAGE 00028**2 C ISF1 - ISF3 TEMPORARY DISK STORAGE 00029**2 C 00030**2 DEFINE FILE 97(400,55,L,LOCREC) , 98(25000,90,L,LPRSUM) 00031 C 00032 INTEGER SN , DELETE , SELECT , HVEC(20) , IAREA(4) 00033 C 00034 COMMON NEDR , IREC , IFORM , LSMALT , LSIG , LAT , LONG , 00035 # LTROP ,LGEOID , LTIDE , LSEAH , IRES , LH 00036 C 00037 COMMON /TIDEM/ C , COVAR , CA , COA , CB , COB , 00038 # YLAT , YLON , ALON , XLAT , XLON , NCONS , 00039 # L , N , T , HTIDE 00040 COMMON /KPAR/ K1 , K2 , K3 00041 COMMON /REFELL/ AE , BE 00042**2 COMMON /AREA/ XLAT1,XLAT2,XLON1,XLON2 /CUTOFF/ SSHMN , SSHMX 00043**3 C 00044 COMMON /SOLVE/ ASN , ANS , F , TMSN , TMNS , SN , NS , NERP , 00045 # B1 , B2 , ISUM , XPX , KK , KL , NSNP , NNSP 00046 # ,IWRITP 00047**2 COMMON /BIASES/ SIGDH , SIGDV , SIGDA , SIGALT 00048**2 EQUIVALENCE (LDATE1,LDATE(1)) , (LDATE2,LDATE(2)) , 00049 # (LTIME1,LTIME(1)) , (LTIME2,LTIME(2)) , 00050 # (FRCT1,FRCT(1)) , (FRCT2,FRCT(2)) 00051 C 00052 EQUIVALENCE (K1,KC(1)) , (IAREA(1),IHEAD(1)) , (IRCF(1),HVEC(1)) ,00053**2 # (NEDR,IRCD(1)) 00054**2 C 00055 DIMENSION C(4,16,2) , COVAR(4,16,2) , CA(4,16,2) , COA(4,16,2) , 00056 # CB(4,16,2) , COB(4,16,2) , YLAT(4,2) , YLON(4,2) , 00057 # ALON(4) 00058 DIMENSION IHEAD(10) , IRCD(13) , LREV(2) , LFILE(2) , LDATE(2) , 00059 # LTIME(2) , FRCT(2) , XSSH(2) , IRCF(20) 00060**2 C 00061**2 DIMENSION BLAT(55) , BLON(55) , B1(3) , KC(3) 00062**4 C 00063**2 C THE DIMENSIONS OF THE FOLLOWING ARRAYS AND VECTORS LIMIT THE 00064**2 C MAXIMUM SIZE OF THE ADJUSTED NETWORK OF SATELLITE ARCS 00065**2 C 00066**2 C WORKSPACE IS LIMITED BY : 00067**2 C MXOBS - THE MAXIMUM NUMBER OF TOTAL OBSERVATIONS 00068**2 C MXOPP - THE MAXIMUM NUMBER OF OBERVATIONS PER PASS 00069**2 C MNUP - THE MAXIMUM TOTAL NUMBER OF ORBITAL ERROR 00070**2 C PARAMETERS CORRESPONDING TO THE N/S PASSES 00071**2 C 00072**2 C ALL THESE VALUES MUST BE SPECIFIED BEFORE CALLING SOLVEQ 00073**2 C CURRENT VALUES ARE : MXOBS = 700 , MXOPP = 50 , MNUP = 150 00074**2 C 00075**2 C TO INCREASE THE REQUIRED WORKSPACE CHANGE THE FOLLOWING 00076**2 C DIMENSIONS BELOW AND IN SUBROUTINE SOLVEQ 00077**2 C 00078**2 C CHANGE MXOBS BELOW 00079**2 DIMENSION F(700) , ASN(700) , ANS(700) 00080**2 C 00081**2 C CHANGE MXOPP BELOW 00082**2 DIMENSION SN(50,4) , NS(50,4) , TMSN(50) , TMNS(50) 00083**2 C 00084**2 C CHANGE MNUP BELOW 00085**2 DIMENSION B2(150) 00086**2 C 00087**2 C 00088 DATA IM2/4HM2 /, IS2/4HS2 /, IN2/4HN2 /, IK1/4HK1 / 00089 DATA MNERP , NXL , IRTN , NP , NXOBS / 3 , 4*0 / 00090**2 DATA INSSH, ICDTA, ICR, IPRT, IPCH, IOUT, ITF1, ITF2, ITF3 00091**2 # / 1 , 3 , 5 , 6 , 7 , 8 , 15 , 16 , 17 / 00092**2 C 00093 NSNP = 0 00094 NNSP = 0 00095 C 00096**2 C DEFINE WORKSPACE ALLOCATION PARAMETERS 00097**3 C 00098**3 MXOBS = 700 00099 MNUP = 150 00100 MXOPP = 50 00101 C 00102 C 00103 READ(ICR,1003) IAREA , XLAT1 , XLAT2 , XLON1 , XLON2,IBDRY 00104**4 READ(ICR,1004) NFILE, IDEL, ISEL, IPRINT, IWRITP, IORBC, 00105 # N, L, NCONS, NDAYP 00106 C 00107 IDEL = IDEL + 1 00108 ISEL = ISEL + 1 00109 C 00110 READ(ICR,1005) NERP, MINX, SSHMN, SSHMX, AE, BE, 00111 # SIGDH, SIGDV, SIGDA , SIGALT 00112**2 READ(ICR,1006) ALPHA, NUM, RMEAN, RMS 00113 GO TO (5,6) , IDEL 00114 6 N1 = 1 00115 N2 = 10 00116 READ(ICR,1007)(IRCF(I) , I=N1,N2) 00117 IF(IRCF(N2).NE.0) GO TO 7 00118 8 N2 = N2 - 1 00119 IF(IRCF(N2).EQ.0) GO TO 8 00120 7 NDEL = N2 00121 GO TO 9 00122 5 NDEL = 0 00123 9 CONTINUE 00124 C 00125 GO TO (18,19) , ISEL 00126 19 N1 = 1 00127 N2 = 10 00128 1 READ(ICR,1008) (ANS(I) , I=N1,N2) 00129 IF(ANS(N2) .EQ. 0) GO TO 2 00130 N1 = N2 + 1 00131 N2 = N2 + 10 00132 GO TO 1 00133 2 N2 = N2 - 1 00134 IF(ANS(N2).EQ.0) GO TO 2 00135 NSEL = N2 00136 GO TO 29 00137 18 NSEL = 0 00138 29 CONTINUE 00139 C 00140**4 IF(IBDRY.EQ.0) GO TO 43 00141**4 C 00142**4 C 00143**4 C THE FOLLOWING INFORMATION PERTAINS TO A POLYGON OF POINTS 00144**4 C SPECIFYING THE AREA FROM WHICH DATA IS DESIRED 00145**4 C (E.G. COASTLINE BOUNDARIES CAN BE DEFINED IN THIS MANNER) 00146**4 C 00147**4 C NPOINT - THE NUMBER OF POINTS IN THE POLYGON. THIS MUST BE A 00148**4 C CLOSED POLYGON , THEREFORE THE FIRST AND THE LAST SETS 00149**4 C OF COORDINATES MUST BE THE SAME. 00150**4 C BLAT - LATITUDE OF POINTS IN DECIMAL DEGREES 00151**4 C BLON - LONGIDUTE OF POINTS IN DECIMAL DEGREES 00152**4 C 00153**4 C (IF IBDRY.EQ.0 SELECTION IS MADE WITHIN THE BOUNDS 00154**4 C XLAT1 , XLAT2 AND XLON1 , XLON2) 00155**4 C 00156**4 I = 0 00157**4 41 READ(ICR,1009,END=42) GLAT1,GLAT2,GLON1,GLON2 00158**4 I = I + 1 00159**4 1009 FORMAT(8X,2F2.0,5X,F3.0,F2.0) 00160**4 BLAT(I) = GLAT1 + GLAT2/60. 00161**4 BLON(I) = GLON1 + GLON2/60. 00162**4 GO TO 41 00163**4 42 BLAT(I) = BLAT(1) 00164**4 BLON(I) = BLON(1) 00165**4 43 CONTINUE 00166**4 NPOINT = I 00167**4 C 00168**4 WRITE(IOUT,2000) IAREA ,XLAT1,XLAT2,XLON1,XLON2 00169**2 WRITE(IOUT,2001) IM2,IS2,IN2,IK1 00170**2 WRITE(IOUT,2002) SIGDH 00171**2 IF(NERP.GE.2) WRITE(IOUT,2003) SIGDV 00172**2 IF(NERP.EQ.3) WRITE(IOUT,2004) SIGDA 00173**2 WRITE(IOUT,2006) SIGALT 00174**2 WRITE(IOUT,2005) AE,BE 00175**2 C 00176**2 C READ COEFFICIENTS OF TIDAL MODEL 00177**2 C 00178**2 CALL TIDE 00179**2 C 00180**2 SIGDH = 1./SIGDH**2 00181**2 SIGDV = 1./SIGDV**2 00182**2 SIGDA = 1./SIGDA**2 00183**2 C 00184**2 W = 1. 00185**2 IF(SIGALT.NE.0.) W = 1./SQRT(SIGALT**2+SIGALT**2) 00186**2 C 00187**2 C 00188 DO 10 I = 1,NFILE 00189 LOCREC = I 00190 READ(97'LOCREC,1000) IFILE,LSREC,LEREC,IREV,IDATE,NX,LEQXT,XLNODE 00191 IF(NXL.EQ.0) N1 = NDEL 00192 IF (NX.GT.MINX) GO TO 11 00193 IF(NX) 10,10,14 00194 14 NXL = NXL + 1 00195 N1 = NDEL + NXL 00196 IRCF(N1) = IFILE 00197 11 CONTINUE 00198 NUMB = 1 00199 IF(N1.NE.0) NUMB = DELETE(IFILE,IRCF,N1) 00200 IF(NUMB.EQ.2) GO TO 10 00201 I1 = LSREC 00202 I2 = LEREC 00203 READ(98'I1,1001) (IRCD(J),J=1,3) , FODAY , (IRCD(J),J=4,13) 00204 T1 = FODAY 00205 SLAT = FLOAT(LAT) / 1.E2 00206 READ(98'I2,1001) (IRCD(J),J=1,3) , FODAY , (IRCD(J),J=4,13) 00207 T2 = FODAY 00208 ELAT = FLOAT(LAT) / 1.E2 00209 TM = T1 + (T2-T1) / 2.E0 00210 C 00211**3 C CHECK WHETHER PASS IS ASCENDING OR DESCENDING 00212**3 C 00213**3 IF(SLAT .LT. ELAT) GO TO 12 00214 NNSP = NNSP + 1 00215 NS(NNSP,1) = IFILE 00216 NS(NNSP,3) = IREV 00217 NS(NNSP,4) = 0 00218 TMNS(NNSP) = TM 00219 GO TO 13 00220 12 NSNP = NSNP + 1 00221 SN(NSNP,1) = IFILE 00222 SN(NSNP,3) = IREV 00223 SN(NSNP,4) = 0 00224 TMSN(NSNP) = TM 00225 13 CONTINUE 00226 10 CONTINUE 00227 ALF = (1.-ALPHA)*100. 00228 NDF = NUM - 1 00229 A1 = ALPHA/2. 00230 CALL TAURE(NUM,NDF,A1,C1) 00231**4 RMS = RMS/W 00232**4 FRC = C1 * RMS 00233**4 C 00234 KK = NERP * NNSP 00235 KL = NERP * NSNP 00236 20 IC = 0 00237 C 00238 21 READ(INSSH,1002,END=800) 00239**2 # (LREV(I),LFILE(I),I=1,2),CRLAT,CRLON, 00240**2 # (LDATE(I),LTIME(I),FRCT(I),I=1,2) , 00241 # (XSSH(I),I=1,2) 00242 IF(LFILE(1).GT.NFILE) GO TO 800 00243 IF(LFILE(2).GT.NFILE) GO TO 21 00244 IF(IC) 23,22,23 00245 22 IC = 1 00246 NOBSP = 0 00247 GO TO 31 00248 23 IF (LREV(1)-IRV) 30,31,30 00249 31 CONTINUE 00250 XLAT = CRLAT 00251 XLON = CRLON 00252 NP = NP + 1 00253 C 00254 IF(IPRINT.EQ.1) WRITE(IPRT,2010) NP,XLAT,XLON 00255**2 C 00256**4 C CHECK WHETHER CROSSOVER POINT IS WITHIN SPECIFIED AREA 00257**4 C 00258**4 IF(IBDRY) 36,36,37 00259**4 36 IF(XLAT.LT.XLAT1 .OR. XLAT.GT.XLAT2) GO TO 226 00260**4 IF(XLON.LT.XLON1 .OR. XLON.GT.XLON2) GO TO 226 00261**4 GO TO 44 00262**4 37 CALL INPOLY(BLAT,BLON,NPOINT,XLAT,XLON,IPOINT) 00263**4 IF(IPOINT.EQ.0) GO TO 226 00264**4 44 CONTINUE 00265**4 C 00266**4 DO 24 I = 1,2 00267 C 00268 C DELETE ASSIGNED REVOLUTIONS FROM THIS RUN 00269 C 00270 IF(N1.EQ.0) GO TO 32 00271 NUMB1 = DELETE(LFILE(1),IRCF,N1) 00272 NUMB2 = DELETE(LFILE(2),IRCF,N1) 00273 NUMB = NUMB1 * NUMB2 00274 IF(NUMB.NE.1) GO TO 25 00275 32 CONTINUE 00276 C 00277 C 00278 C PASS CROSSOVER SSH THROUGH THE INNER AND OUTER CUTOFF WINDOWS 00279 C 00280 IF(XSSH(1).LT.SSHMN .OR. XSSH(1).GT.SSHMX .OR. 00281 # XSSH(2).LT.SSHMN .OR. XSSH(2).GT.SSHMX) GO TO 225 00282**4 C 00283 C 00284 NSELCT = 1 00285 IF(I.EQ.1 .AND. NSEL.NE.0) 00286 # NSELCT = SELECT(LREV(1),LREV(2),ANS,NSEL) 00287 IF(NSELCT.EQ.2) GO TO 25 00288 C 00289 C 00290 C TEST OBSERVATIONS FOR WITHIN CONTEXT OUTLIERS 00291**4 C 00292 DSSH = ABS(XSSH(1)-XSSH(2)-RMEAN) 00293 IF(DSSH - FRC) 26,26,27 00294 27 IF(IPRINT.EQ.1) WRITE(IPRT,2014) DSSH,ALF 00295**2 GO TO 25 00296 C 00297 26 CONTINUE 00298 GO TO (50,51) , I 00299 50 NXOBS = NXOBS + 1 00300 NOBSP = NOBSP + 1 00301 C 00302 C 00303 C SEPARATE S/N AND N/S PASSES 00304 C 00305 DO 52 J = 1,NSNP 00306 IF(LFILE(I).NE.SN(J,1)) GO TO 52 00307 KPASS = J 00308 TM = TMSN(J) 00309 IREV1 = SN(KPASS,3) 00310 IFILE1 = SN(KPASS,1) 00311 SN(KPASS,4) = SN(KPASS,4) + 1 00312 GO TO 53 00313 52 CONTINUE 00314 51 CONTINUE 00315 C 00316 DO 54 J = 1,NNSP 00317 IF(LFILE(I).NE.NS(J,1)) GO TO 54 00318 KPASS = J 00319 TM = TMNS(J) 00320 IREV2 = NS(KPASS,3) 00321 IFILE2 = NS(KPASS,1) 00322 NS(KPASS,4) = NS(KPASS,4) + 1 00323 GO TO 53 00324 54 CONTINUE 00325 53 CONTINUE 00326 C 00327 C CORRECT SSH FOR TIDAL EFFECT 00328 C 00329 C 00330 CALL YRMDAY (LDATE(I),LYR,LMON,LDAY) 00331 CALL HHMMSS (TM,LHR,LMIN,SECL) 00332 CALL ELAPST (NDAYP,LYR,LMON,LDAY,LHR,LMIN,SECL,T) 00333**2 CALL TIDEIN 00334 C 00335 SSHC = XSSH(I) - HTIDE/1.E2 00336 C 00337 GO TO (15,16) , I 00338 15 SSHC1 = SSHC 00339 HTIDE1 = HTIDE / 1.E2 00340 IF(IPRINT.EQ.1) WRITE(IPRT,2013) 00341**2 IF(IPRINT.EQ.1) WRITE(IPRT,2011) LREV(I),LDATE(I),LTIME(I), 00342**2 # FRCT(I),SSHC1,HTIDE1 00343**2 GO TO 17 00344 C 00345 16 SSHC2 = SSHC 00346 HTIDE2 = HTIDE / 1.E2 00347 IF(IPRINT.EQ.1) WRITE(IPRT,2012) LREV(I),LDATE(I),LTIME(I), 00348**2 # FRCT(I),SSHC2,HTIDE2,NXOBS 00349**2 17 CONTINUE 00350 C 00351 NERP1 = NERP - 1 00352 KC(1) = NERP * (KPASS-1) + 1 00353 IF(NERP1) 60,60,61 00354 61 CONTINUE 00355 DO 28 M = 1,NERP1 00356 MP1 = M + 1 00357 28 KC(MP1) = KC(1) + M 00358 60 CONTINUE 00359 GO TO (58,59) , I 00360 58 CALL YRMDAY (LTIME(I),JHR,JMIN,JSEC) 00361 TIM = FLOAT(JHR)*60. + FLOAT(JMIN) + (FLOAT(JSEC) + FRCT(I))/60. 00362 DT = TIM - TM/6.E7 00363 C 00364 C 00365 C FORM DESIGN MATRIX B PARTISITIONING FOR S/N AND N/S PASSES 00366 C OUTPUT PARTISIONED MATRICES IN TEMPORARY FILES ITF1 , ITF2 00367 C (ROW BY ROW STORAGE MODE) 00368 CALL CLEAR (B1,KL) 00369 GO TO (35,34,33) , NERP 00370 33 B1(3) = DT ** 2 00371 34 B1(2) = DT 00372 35 B1(1) = 1.E0 00373 WRITE(ITF1) (B1(IB) , IB=1,NERP) 00374 SN(KPASS,2) = NOBSP 00375 GO TO 24 00376 59 CALL YRMDAY (LTIME(I),JHR,JMIN,JSEC) 00377 TIM = FLOAT(JHR)*60. + FLOAT(JMIN) + (FLOAT(JSEC) + FRCT(I))/60. 00378 DT = TIM - TM/6.E7 00379 CALL CLEAR (B2,KK) 00380 GO TO (40,39,38) , NERP 00381 38 B2(K3) = - DT ** 2 00382 39 B2(K2) = - DT 00383 40 B2(K1) = - 1.E0 00384 WRITE(ITF2) (B2(IB0) , IB0 = 1,KK) 00385 IREV = IREV1 * 10000 + IREV2 00386 IFILE = IFILE1 * 10000 + IFILE2 00387 C 00388 WRITE(ITF3) IREV,IFILE,XLAT,XLON,LDATE1,LTIME1,FRCT1, 00389 # LDATE2,LTIME2,FRCT2,SSHC1,SSHC2 00390 C 00391 C 00392 C FORM MISCLOSURE VECTOR F 00393 C 00394 F(NXOBS) = SSHC1 - SSHC2 00395 24 CONTINUE 00396 GO TO 25 00397**4 225 IF(IPRINT.EQ.1) 00398**4 # WRITE(IPRT,2021) (XSSH(I),LDATE(I),LTIME(I),FRCT(I),I=1,2) 00399**4 GO TO 25 00400**4 226 IF(IPRINT.EQ.1) WRITE(IPRT,2022) 00401**4 25 CONTINUE 00402 C 00403 IC = IC + 1 00404 IRV = LREV(1) 00405 GO TO 21 00406 30 IC = 0 00407 IRV = LREV(1) 00408 GO TO 22 00409 800 CONTINUE 00410 VSUM = 0.0 00411 SUM = 0.0 00412 REWIND ITF1 00413 REWIND ITF2 00414 REWIND ITF3 00415 C 00416 C 00417 C CHECK WHETHER THERE ARE ENOUGH CROSSOVER POINTS FOR EACH PASS 00418 C (THERE SHOULD BE AT LEAST 'NEPR' CROSSOVER POINTS) 00419 00420 DO 3 I = 1,NSNP 00421 IF(SN(I,4).GE.NERP) GO TO 3 00422 WRITE(IPRT,2009) NERP,SN(I,3),SN(I,4) 00423**2 IRTN = 1 00424 3 CONTINUE 00425 C 00426 DO 4 I = 1,NNSP 00427 IF(NS(I,4).GE.NERP) GO TO 4 00428 WRITE(IPRT,2009) NERP,NS(I,3),NS(I,4) 00429**2 IRTN = 1 00430 4 CONTINUE 00431 IF(IRTN .EQ. 1) STOP 00432 C 00433 C 00434 C BEGIN ADJUSTMENT OF INTERSECTING ARCS 00435 CALL SOLVEQ(MXOBS,MXOPP,MNERP,MNUP) 00436 C 00437 C 00438 C APPLY CORRECTIONS TO THE DATA USING THE RECOVERED ORBITAL 00439 C ERROR BIASES. COMPUTE ADJUSTED SSH AT THE CROSSOVER POINTS 00440 C 00441 DO 55 I = 1,ISUM 00442 READ(ITF3) IREV,IFILE,XLAT,XLON,LDATE1,LTIME1,FRCT1, 00443 # LDATE2,LTIME2,FRCT2,SSHC1,SSHC2 00444 VP1 = SSHC1 - ASN(I) 00445 VP2 = SSHC2 + ANS(I) 00446 ASN(I) = ASN(I) + ANS(I) - F(I) 00447 ANS(I) = IREV1*10000. + IREV2 00448 FRCT1 = FRCT1 + FLOAT(LTIME1) 00449 FRCT2 = FRCT2 + FLOAT(LTIME2) 00450 VSUM = VSUM + ASN(I) 00451 WRITE(ITF2) IREV ,IFILE,XLAT,XLON,LDATE1,FRCT1,LDATE2,FRCT2, 00452 # SSHC1,SSHC2,VP1,VP2 00453 IREV1 = IREV/10000 00454 IREV2 = IREV - IREV1*10000 00455 55 CONTINUE 00456 C 00457 C COMPUTE MEAN AND RMS DIFFERENCE OF THE ADJUSTED 00458 C SEA SURFACE HEIGHTS AT THE CROSSOVER POINTS 00459 C 00460 VTPV = 0. 00461**2 VM = VSUM / ISUM 00462 DO 56 I = 1,ISUM 00463 VTPV = VTPV + (ASN(I)*W)**2 00464**3 56 SUM = SUM + (ASN(I))**2 00465**3 RMS = SQRT(SUM/ISUM) 00466**3 SIG = (VTPV+XPX)/(2.*ISUM) 00467**3 SVI = SQRT(SIG)/W 00468**3 SIGMSS = SQRT(SIG/2.) 00469**4 C 00470 C 00471 C CHECK FOR OUTLIERS WITHIN THE CONTEXT OF THE ADJUSTED 00472 C OBSERVATION SERIES (ADJUSTED DSSH) 00473 C 00474 A1 = A1 / ISUM 00475**3 CALL TAURE (1,ISUM,A1,C1) 00476**3 TESTF = C1*SVI 00477**3 NSEL = 0 00478 AVM = VM 00479 NN = ISUM 00480 WRITE(IOUT,2016) 00481**2 REWIND ITF2 00482 WRITE(IPRT,2020) 00483**2 DO 100 I = 1,ISUM 00484 F(I) = 0. 00485 IF( ABS(ASN(I) ) - TESTF) 109,109,101 00486**3 101 NSEL = NSEL + 1 00487 ANS(NSEL) = ANS(I) 00488 WRITE(IPRT,2015) I,ALF,ASN(I) 00489**2 F(I) = FLOAT(I) 00490 AAZ = NN*AVM - ASN(I) 00491 NN = NN - 1 00492 AVM = AAZ / NN 00493 109 READ(ITF2) IREV ,IFILE,XLAT,XLON,LDATE1,FRCT1,LDATE2,FRCT2, 00494 # SSHC1,SSHC2,VP1,VP2 00495 AVG = (VP1 + VP2) / 2. 00496 RES1 = VP1 - AVG 00497 RES2 = VP2 - AVG 00498 IREV1 = IREV/10000 00499 IREV2 = IREV - IREV1*10000 00500 WRITE(IOUT,2017) IREV1,IREV2,XLAT,XLON,SSHC1,SSHC2,VP1,VP2, 00501**2 # AVG,RES1,RES2 00502 IFILE1 = IFILE/10000 00503 IFILE2 = IFILE - IFILE1*10000 00504 IF(IWRITP.EQ.1) 00505**2 # WRITE(ICDTA) IREV1,IFILE1,IREV2,IFILE2,XLAT,XLON,LDATE1, 00506**2 # FRCT1,LDATE2,FRCT2,VP1,VP2,AVG 00507 100 CONTINUE 00508 C 00509 END FILE ICDTA 00510**2 C 00511 WRITE(IPRT,2018) VTPV,XPX,ISUM,SIG,SIGMSS 00512**4 IF(NSEL.EQ.0) GO TO 901 00513 N1 = 1 00514 N2 = 10 00515 IF(NSEL.LT.N2) N2 = NSEL 00516 104 CONTINUE 00517 DO 102 I = N1,N2 00518 IHEAD(I) = ANS(I) 00519 102 CONTINUE 00520 IF(N2.GE.NSEL) GO TO 103 00521 WRITE(IPCH,1004) (IHEAD(I),I=N1,N2) 00522**2 N1 = N2 + 1 00523 N2 = N2 + 10 00524 IF(NSEL.LT.N2) N2 = NSEL 00525 GO TO 104 00526 103 WRITE(IPCH,1004) (IHEAD(I),I=N1,N2) 00527**2 VV = 0. 00528 DO 105 I = 1,ISUM 00529 ICHECK = INT(F(I)) 00530 IF(I.NE.ICHECK) GO TO 106 00531 GO TO 105 00532 106 VN = AVM - ASN(I) 00533 V2 = VN*VN 00534 VV = VV + V2 00535 105 CONTINUE 00536 STD = SQRT( VV/FLOAT(NN) ) 00537 901 CALL SORT (ASN,MXOBS,ISUM) 00538**4 NCLASS = 20 00539 DO 120 J = 1,NCLASS 00540 120 HVEC(J) = 0 00541 XX = -10 00542**4 J = 1 00543 DO 121 K = 1,NCLASS 00544 124 IF(ASN(J).GT.XX .AND. K.LT.NCLASS) GO TO 122 00545 HVEC(K) = HVEC(K) + 1 00546 J = J + 1 00547 IF(J.GT.ISUM) GO TO 123 00548 GO TO 124 00549 122 XX = XX + 1. 00550 121 CONTINUE 00551 123 CONTINUE 00552 CALL PLOT (NCLASS,HVEC) 00553 WRITE(IPRT,2019) VM,AVM,RMS,STD 00554**2 C 00555 1000 FORMAT(7I6,F13.6) 00556**2 1001 FORMAT(2I4,I3,F12.0,I9,I6,2I8,6I6) 00557**2 1002 FORMAT(4I5,2F10.6,2(I8,I6,F7.6),2F11.6) 00558**2 1003 FORMAT(4A4,4F4.0,I4) 00559**4 1004 FORMAT(10I8) 00560 1005 FORMAT(2I4,2F6.0,2F12.2,4F8.4) 00561**2 1006 FORMAT(F6.3,I6,2F6.2) 00562 1007 FORMAT(10I4) 00563 1008 FORMAT(10F8.0) 00564 2000 FORMAT('1',46X,'MEAN SEA SURFACE COMPUTATION',//,47X,'USING GEOS-300565 # ALTIMETER DATA',///,36X,'LEAST SQUARES ADJUSTMENT OF SEA SURFACE'00566 #,' HEIGHTS',//,36X,'COMPUTED AT INTERSECTING POINTS OF CROSSING ',00567 #'ARCS ',////,53X,4A4,///,41X,'LATITUDE',4X,F5.2,2X, 00568 #'TO',3X,F5.2,2X,'DEGREES',//,41X,'LONGITUDE',2X,F6.2,2X,'TO',2X, 00569 #F6.2,2X,'DEGREES',/////,10X,'OPTIONS USED FOR THIS RUN',//) 00570 2001 FORMAT(20X,'ANALYTICAL TIDAL MODELLING FOR ',4A4 ,////,20X, 00571 # 'ORBITAL ERROR MODELLING FOR ',30X,'ORBITAL CONSTRAINTS',//) 00572 2002 FORMAT(30X,'ABSOLUTE BIAS ',6X,'DH (METERS)',12X,F10.2, 00573 # 4X,'(METERS)' ) 00574 2003 FORMAT(/,30X,'TILT',16X,'DV (METERS/MIN)' ,8X,F10.3,3X, 00575 # ' (METERS/MIN)',/) 00576 2004 FORMAT(30X,'SECOND ORDER EFFECT DA (METERS/MIN**2)', 5X,F10.3,3X, 00577 # ' (METERS/MIN**2)',//) 00578 2005 FORMAT( ///,20X,'REFERENCE ELLIPSOID',//,30X,'A = ',F12.2,5X, 00579 # 'B = ',F12.2,//) 00580 2006 FORMAT(///,20X,'ALTIMETRY MEASUREMENT WEIGHTING (SIGALT = 0 - IDE00581**2 #NTITY WEIGHTING)',//,30X,'SIGALT = ',F5.2,' (METERS)') 00582**2 2009 FORMAT(10X,'FEWER THAN ',I2,' ARCS INTERSECT REV #',I5, 00583 # ' NX = ',I2,/) 00584 2010 FORMAT(//,100X,'CROSSOVER POINT # ',I4,/,5X,'LAT : ',F10.6,10X, 00585 # 'LONG : ',F10.6,/) 00586 2011 FORMAT(5X,'S/N REV : ',I4,10X,I6,16X,I6,F5.4,16X,F10.4,8X, 00587 # F10.4,/) 00588 2012 FORMAT(5X,'N/S REV : ',I4,10X,I6,16X,I6,F5.4,16X,F10.4,8X,F10.4, 00589 # //,5X,45('-'),' ACCEPTED CROSSOVER POINT # ',I4,1X,45('-') ) 00590 2013 FORMAT(30X,'DATE',20X,'TIME',20X,'A-PRIORI SSH' 00591 # ,10X,'TIDE',/,28X,'(YYMMDD)',14X,'(HHMMSS.SSSS)',17X,'(METERS)'00592 # ,10X,'(METERS)',/) 00593 2014 FORMAT(/,5X,'DSSH (METERS) : ',F10.4,' TAU-TEST FAILED AT ',F3.0, 00594 # ' % - OBSERVATION REJECTED' , /) 00595 2015 FORMAT(10X,'*** DSSH AT CROS. POINT # ',I4,' FAILED MAX-TAU TEST 00596 #AT ',F3.0,' % - DSSH (METERS) : ',F10.4,' ***') 00597 2016 FORMAT('1',35X,'A D J U S T E D S E A S U R F A C E H E I G H T00598 # S',///,30X,'CROSSING POINT',8X,'A-PRIORI SSH',8X, 00599 # 'ADJUSTED SSH',8X,'AVERAGE',5X,2('RESIDUAL',3X), 00600 #/,5X,'S/N REV. N/S REV.',3X,'LATITUDE',3X,'LONGITUDE',4X,2(2X,'S/00601 #N',6X,'N/S',6X),'ADJUSTED SSH',4X,'S/N',8X,'N/S',/,26X, 00602 # 2('(DEGREES)',2X),1X,4('(METERS)',2X),3X,'(METERS)',4X, 00603 # 2('(METERS)',3X),//) 00604 2017 FORMAT(7X,I4,6X,I4,3X,2(2X,F9.4),1X,4(2X,F8.4),4X,F8.4,4X, 00605**3 # 2(F8.4,3X)) 00606 2018 FORMAT(///,10X,'VTPV = ',F12.4,//,10X,'XTPX = ',F12.4,//,10X, 00607 # 'DEGREES OF FREEDOM = ',I4,//,10X,'ESTIMATED VARIANCE FACTOR ='00608 # ,F12.4,//,10X,'(ALL OBSERVATIONS INCLUDED)',///, 00609**4 # 10X,'ESTIMATED ACCURACY OF MEAN SEA SURFACE (METERS) = ',F8.4) 00610**4 2019 FORMAT(///,30X,'HISTOGRAM OF ADJUSTED CROSSOVER SEA SURFACE HEIGHT00611 # DIFFERENCES',///,10X,'MEAN DIFF. = ',F12.4,30X,'MEAN DIFF. = ' 00612 # ,F12.4,//,10X,'RMS DIFF. = ',F12.4,30X,'RMS DIFF. = ',F12.4, 00613 # //,10X,'(ALL OBSERVATIONS INCLUDED)',28X, 00614 # '(REJECTED OBSERVATIONS EXCLUDED)',//) 00615 2020 FORMAT('1',35X,'ADJUSTMENT STATISTICS SUMMARY',///) 00616 2021 FORMAT(/,5X,'SSH OUTSIDE CUTOFF LIMITS (SSH , DATE , TIME) : ', 00617**4 # 2(T60,F10.4,5X,I8,2X,I8,F5.4,/)) 00618**4 2022 FORMAT(/,5X,'CROSSOVER POINT OUTSIDE SPECIFIED COASTLINE BOUNDARIE00619**4 #S',/) 00620**4 C 00621 900 STOP 00622 END 00623 SUBROUTINE TIDE C 00625**2 C NAME TIDE 00626**2 C 00627**2 C PURPOSE ANALYTICAL EVALUATION OF THE TIDAL EFFECT OF 00628**2 C THE HUDSON AND JAMES BAY 00629**2 C 00630**2 C I/O DEVICES INTIDE 00631**3 C 00632**3 C CALLING ROUTINES CARTE , PREDCT 00633**3 C 00634**3 C 00635**2 C SUBROUTINE 'TIDE' MUST BE CALLED ONCE TO ESTABLISH THE COEFFICIENTS 00636**2 C FOR THE POLYNOMIAL APPROXIMATION OF AMPLITUDES AND PHASES OF THE 00637**2 C M2 , S2 , N2 , K1 TIDAL CONSTITUENTS. THEREAFTER 'TIDEIN' 00638**2 C MAY BE USED. 00639**2 C 00640**2 COMMON /TIDEM/ C , COVAR , CA , COA , CB , COB , 00641 # YLAT , YLON , ALON , XLAT , XLON , NCONS , 00642 # L , N , T , HTIDE 00643 DIMENSION C(4,16,2) , COVAR(4,16,2) , CA(4,16,2) , COA(4,16,2) , 00645 # CB(4,16,2) , COB(4,16,2) , YLAT(4,2) , YLON(4,2) , 00646 # ALON(4) , LINE (4) , WFREQ(4) , FK(4) 00647 DIMENSION PM(16) , H(16) , GA(16) , GB(16) 00648 COMMON /REFELL/ RA , RB 00649**2 COMMON /BLOCK/ WFREQ , FK 00650**3 DATA IN2 /4HN2 / , IS2 /4HS2 / , IK1 /4HK1 / ,IW /4HWEST/ , 00651**2 # IBLK / 4H / , CK / 0.872664626E-2 / , INTIDE / 2 / 00652**2 C 00644 C 00653**2 PI = ARCOS(-1.E0) 00654 C 00655 200 READ(INTIDE,1000,END=210) KTIDE , LINE , XLAT1 , XLAT2 , XLON1 , 00656**2 # XLON2 , ORLAT , ORLON , OMEGA , FNODE 00657 C 00658 M = 1 00659 IF (KTIDE .EQ. IN2) M = 2 00660 IF (KTIDE .EQ. IS2) M = 3 00661 IF (KTIDE .EQ. IK1) M = 4 00662 K = 2 00663 IF (LINE(4) .EQ. IW .OR. LINE(4) .EQ. IBLK) K = 1 00664 READ(INTIDE,1001,END=210) (C(M,J,K),J=1,L) , (COVAR(M,J,K),J=1,L) 00665**2 # , (CA(M,J,K),J=1,L) , (COA(M,J,K),J=1,L) 00666 # , (CB(M,J,K),J=1,L) , (COB(M,J,K),J=1,L) 00667 GO TO (220,230) , K 00668 220 ALON(M) = XLON2 00669 YLAT(M,1) = ORLAT 00670 YLON(M,1) = ORLON 00671 WFREQ(M) = OMEGA 00672 FK(M) = FNODE 00673 GO TO 240 00674 230 YLAT(M,2) = ORLAT 00675 YLON(M,2) = ORLON 00676 240 CONTINUE 00677 GO TO 200 00678 210 CONTINUE 00679 1000 FORMAT(5A4,4F6.0,2F10.3,F10.6,F6.3) 00680 1001 FORMAT(4E20.10) 00681 RETURN 00682 C 00683 C 00684 ENTRY TIDEIN 00685 C 00686 HTIDE = 0.E0 00687 IDP = N + 1 00688 DO 10 I = 1,NCONS 00689 K = 2 00690 IF(XLON .LE. ALON(I)) K = 1 00691 CALL CARTE (XLAT,XLON,YLAT(I,K),YLON(I,K),XP,YP) 00692 C 00693 DO 20 J = 1,L 00694 20 PM(J) = 0.E0 00695 ICDP = 0 00696 DO 21 JK = 1,IDP 00697 KA = JK-1 00698 DO 22 J = 1,IDP 00699 JA = J-1 00700 ICDP = ICDP + 1 00701 PM(ICDP) = XP**KA*YP**JA 00702 22 CONTINUE 00703 21 CONTINUE 00704 C 00705 DO 30 JJ = 1,L 00706 H(JJ) = C(I,JJ,K) 00707 GA(JJ) = CA(I,JJ,K) 00708 GB(JJ) = CB(I,JJ,K) 00709 30 CONTINUE 00710 C 00711 CALL PREDCT (L,PM,H,PH) 00712 CALL PREDCT (L,PM,GA,PGA) 00713 CALL PREDCT (L,PM,GB,PGB) 00714 PG = ATAN (PGB/PGA)/CK 00715 IF (PG. LT. 0.E0) PG = 360.E0 + PG 00716 ARG = WFREQ(I)*T + PG 00717 ARGR = (AMOD(ARG,360.E0))*PI/180.E0 00718 HTIDE = HTIDE + PH*FK(I)*COS(ARGR) 00719 10 CONTINUE 00720 RETURN 00721 END 00722 SUBROUTINE ELAPST (DAYP,IYR,MON,IDAY,IHR,MIN,SEC,T) C 00724**3 C NAME ELAPST 00725**3 C 00726**3 C PURPOSE GIVEN THE DATE IN YEAR, MONTH AND DAY AND THE 00727**3 C TIME IN HOURS, MINUTES AND SECONDS COMPUTE THE 00728**3 C TIME ELAPSED SINCE 0 HOURS OF 'DAYP' (FIRST DAY 00729**3 C OF PREDICTIONS) 00730**3 C 00731**3 INTEGER DAYR(12) , DAYP 00732 C 00733**3 C THE DATA STATEMENT BELOW ACCOUNTS FOR ONE LEAP YEAR - 1976 - ONLY 00734**3 C VECTOR DAYR CONTAINS FOR EACH MONTH THE NUMBER OF 00735**3 C DAYS ELAPSED SINCE THE BEGINNING OF THE YEAR 00736**3 C 00737**3 DATA DAYR /31,60,91,121,152,182,213,244,274,305,335,366 / 00738**2 JDAY = DAYR(MON-1) + IDAY - DAYP 00739 HRS = FLOAT(JDAY)*24.E0 + FLOAT(IHR) + (FLOAT(MIN)+SEC/60.E0)/60. 00740 T = HRS 00741 RETURN 00742 END 00743 SUBROUTINE CARTE (ALAT,ALON,PHIO,ALONO,X,Y) C 00745**2 C NAME CARTE 00746**2 C 00747**2 C PURPOSE CONVERT LATITUDE AND LONGITUDE INTO X AND Y 00748**2 C LOCAL CARTECIAN COORDINATES REFERRED TO A LOCAL 00749**2 C SYSTEM WITH ORIGIN (PHIO,ALONO) TO BE SPECIFIED 00750**2 C 00751**2 COMMON /REFELL/ RA , RB 00752**2 PI = ARCOS(-1.E0) 00753 EC = (RA**2-RB**2)/RA**2 00754 PHIR = PHIO*PI/180.E0 00755 ALONR = ALONO*PI/180.E0 00756 XM = ((1.E0-EC)*RA)/( SQRT(1.E0-EC*(SIN(PHIR)**2))**3) 00757 XN = RA/SQRT(1.E0-EC*(SIN(PHIR)**2)) 00758 R = SQRT(XM*XN) 00759 X = R*(ALAT-PHIO)*PI/180.E0 00760 Y = R*COS(PHIR)*(ALON-ALONO)*PI/180.E0 00761 RETURN 00762 END 00763 SUBROUTINE PREDCT (L,PM,COEF,FUNC) C 00765**2 C NAME PREDCT 00766**2 C 00767**2 C PURPOSE PREDICT THE APPROXIMANT 'FUNC' OF A GIVEN FUNCTION 00768**2 C FROM A GIVEN SET OF POLYNOMIAL COEFFICIENTS 'COEF' 00769**2 C AND THE DEGREE OF THE APPROXIMATION 'L' 00770**2 C 00771**2 DIMENSION PM(L) , COEF(L) 00772 SUM = 0.E0 00773 DO 10 I = 1,L 00774 FUNC = PM(I) * COEF(I) 00775 SUM = SUM + FUNC 00776 10 CONTINUE 00777 FUNC = SUM 00778 RETURN 00779 END 00780 SUBROUTINE YRMDAY (IDATE,IYR,MON,IDAY) C 00782 C 00783**2 C NAME YRMDAY 00784**2 C 00785**2 C PURPOSE CONVERT A DATE GIVEN IN THE FORM 'YYMMDD' INTO 00786**2 C YEAR (IYR) , MONTH (MON) AND DAY (IDAY) 00787**2 C 00788**2 IYR = IDATE / 10000 00789 IX = MOD (IDATE,10000) 00790 MON = IX / 100 00791 IDAY = MOD (IX,100) 00792 C 00793 RETURN 00794 END 00795 SUBROUTINE HHMMSS (TMSEC,IH,IM,SEC) C 00797**3 C NAME HHMMSS 00798**3 C 00799**3 C PURPOSE CONVERT TIME 'TMSEC' IN MICROSECONDS INTO 00800**3 C HOURS , MINUTES AND SECONDS 00801**3 C 00802**3 S = TMSEC / 1.E6 00803 IH = S / 3600.E0 00804 IM = AMOD (S,3600.) / 60. 00805 SEC = AMOD (S,60.) 00806 RETURN 00807 END 00808 SUBROUTINE ZERO(A,N,M) C 00810**2 C NAME ZERO 00811**2 C 00812**2 C PUPROSE CLEAR AN ARRAY AA FOR N TIMES M VALUES 00813**2 C STARTING AT LOCATION A(1) 00814**2 C 00815**2 DIMENSION A(1) 00816 I = N * M 00817 DO 1 J = 1,I 00818 1 A(J) = 0.0 00819 RETURN 00820 END 00821 SUBROUTINE CLEAR (AA,N) C 00823**2 C NAME CLEAR 00824**2 C 00825**2 C PURPOSE CLEAR A VECTOR AA OF SIZE N 00826**2 C 00827**2 DIMENSION AA(1) 00828 DO 1 M = 1 , N 00829 1 AA(M) = 0.0 00830 RETURN 00831 END 00832 INTEGER FUNCTION DELETE (IFILE,IRCF,N) C 00834**3 C NAME DELETE 00835**3 C 00836**3 C TYPE INTEGER FUNCTION 00837**3 C 00838**3 C PURPOSE DELETE AN ASSIGNED REVOLUTION IN FILE 'IFILE' 00839**3 C FROM LEAST SQUARES ADJUSTMENT 00840**3 C 00841**3 DIMENSION IRCF(N) 00842 DELETE = 1 00843 DO 1 I = 1,N 00844 IF (IFILE.EQ.IRCF(I)) GO TO 2 00845 1 CONTINUE 00846 RETURN 00847 2 DELETE = 2 00848 RETURN 00849 END 00850 SUBROUTINE MATMPY(M1,M2,M3,L,M,N,JL,JM,JN,ICODE) 00851 C 00852**2 C NAME MATMPY 00853**2 C 00854**2 C PURPOSE TO COMPUTE THE PRODUCT OF TWO MATRICES IN ANY 00855**2 C ALLOWABLE TRANSPOSE COMBINATION AS FOLLOWS: 00856**2 C 00857**2 C OPTION ICODE PRODUCT M3 00858**2 C 00859**2 C 1 M1*M2 00860**2 C 2 (M1)T*M2 00861**2 C 3 M1*(M2)T 00862**2 C 4 (M1)T*(M2)T 00863**2 C 00864**2 C 00865**2 C (L,M),(M,N),(L,N) ARE THE DIMENSIONS OF THE 00866**2 C PRE- , POST- AND PRODUCT-MATRICES 00867**2 C RESPECTIVELY 00868**2 C JL,JM,JN ARE CORRESPONDING DECLARED ROW 00869**2 C DIMENSIONS AT THE CALLINE PROGRAM 00870**2 C 00871**2 C 00872 REAL M1,M2,M3 00873**2 DIMENSION M1(JL,1) , M2(JM,1) , M3(JN,1) 00874 C 00875 DO 11 I = 1,L 00876 DO 11 J = 1,N 00877 M3(I,J) = 0. 00878 DO 11 K = 1,M 00879 GO TO(1,2,3,4),ICODE 00880 1 CONTINUE 00881 C M3 = M1 * M2 00882 M3(I,J) = M3(I,J) + M1(I,K)*M2(K,J) 00883 GO TO 11 00884 2 CONTINUE 00885 C M3 = M1 TRANSPOSE * M2 00886 M3(I,J) = M3(I,J) + M1(K,I)*M2(K,J) 00887 GO TO 11 00888 3 CONTINUE 00889 C M3 = M1 * M2 TRANSPOSE 00890 M3(I,J) = M3(I,J) + M1(I,K)*M2(J,K) 00891 GO TO 11 00892 4 CONTINUE 00893 C M3 = M1 TRANSPOSE * M2 TRANSPOSE 00894 M3(I,J) = M3(I,J) + M1(K,I)*M2(J,K) 00895 11 CONTINUE 00896 RETURN 00897 END 00898 SUBROUTINE SOLVEQ (MXOBS,MXOPP,MNERP,MNUP) 00899 C 00900**3 C NAME SOLVEQ 00901**3 C 00902**3 C PURPOSE TO CONTROL THE SET-UP OF THE NORMAL EQUATIONS , 00903**3 C THE APPLICATION OF THE ORBITAL CONSTRAINTS AND 00904**3 C THE LEAST SQUARES ESTIMATION PROCESS 00905**3 C 00906**3 C 00907 INTEGER SN 00908 COMMON /SOLVE/ ASN , ANS , F , TMSN , TMNS , SN , NS , NERP , 00909 # B1 , B2 , ISUM , XPX , KK , KL , NSNP , NNSP 00910 # ,IWRITP 00911**2 COMMON /BIASES/ SIGDH , SIGDV , SIGDA , SIGALT 00912**2 C 00913**3 DIMENSION YSN(3) , BTF(3) , B1(3) , DF(3) , ORBIAS(3) , 00914**3 # DG(3) , BN(3,3) , BCOV(3,3) 00915**3 C 00916**3 C IF MAXIMUM DIMENSIONS SPECIFIED BELOW ARE FOUND INSUFFICIENT , 00917**3 C CHANGE THE DIMENSIONS BELOW AND IN THE MAIN PROGRAM 00918**3 C 00919**3 C 00920 C CHANGE MXOBS BELOW 00921 DIMENSION ASN(700) , ANS(700) , F(700) 00922 C 00923 C CHANGE MXOPP BELOW 00924 DIMENSION TMSN(50) , TMNS(50) , SN(50,4) , NS(50,4) , FY(50) 00925 C 00926 C CHANGE MNUP BELOW 00927 DIMENSION BCD(150,150) , DN(150,150) 00928 DIMENSION XSN(150) , XNS(150) , DFN(150) , BTY(150) , 00929 # UU(150) , B2(150) , DN1(150) 00930 C CHANGE MXOPP , MNERP BELOW 00931 DIMENSION BA(50,3) 00932 C 00933 C CHANGE MNERP , MNUP BELOW 00934 DIMENSION BAT(3,150) , BBT(150,3) 00935 C 00936 C CHANGE MXOPP , MNUP BELOW 00937 DIMENSION BB(50,150) 00938 C 00939 C 00940 EQUIVALENCE (BCD(1,1),BN(1,1)) , (DN(1,1),BAT(1,1)) 00941 EQUIVALENCE (BBT(1,1),XSN(1)), (BBT(1,2),XNS(1)), (BBT(1,3),FY(1))00942 EQUIVALENCE (DFN(1),ASN(1)) , (DN1(1),ANS(1),YSN(1)) , 00943 # (ORBIAS(1),SIGDH) 00944 C 00945 C THE EQUIVALENCE STATEMENT BELOW SHOULD BE (UU(1),ANS(MNUP+1)) 00946 EQUIVALENCE (UU(1),ANS(151)) 00947 C 00948 DATA ITF1 , ITF2 , ITF3 , ITF4 / 15 , 16 , 17 , 18 / 00949 DATA ISF1 , ISF2 , ISF3 / 10 , 11 , 12 / 00950**2 DATA IOUT / 8 / 00951**3 C 00952**2 W = 1. 00953**2 IF(SIGALT.NE.0.) W = 1./SQRT(SIGALT**2+SIGALT**2) 00954**2 NERP1 = NERP - 1 00955 MI = 0 00956 MJ = 0 00957 DO 32 IK = 1,NSNP 00958 L = SN(IK,2) 00959 MI = MI + 1 00960 MJ = MJ + L 00961 I = 0 00962 DO 34 IL = MI,MJ 00963 I = I + 1 00964 FY(I) = F(IL)*W 00965**2 J = 0 00966 READ(ITF1) (B1(IB) , IB=1,NERP) 00967 DO 34 JL = 1,KL 00968 J = J + 1 00969 34 BA(I,J) = B1(JL)*W 00970**2 C 00971 CALL MATMPY (BA,BA,BN,NERP,L,NERP,MXOPP,MXOPP,MNERP,2) 00972 CALL MATMPY (BA,FY,BTF,NERP,L,1,MXOPP,MXOPP,MNERP,2) 00973 C 00974 GO TO (23,22,21) , NERP 00975 21 BN(3,3) = BN(3,3) + ORBIAS(3) 00976 22 BN(2,2) = BN(2,2) + ORBIAS(2) 00977 23 BN(1,1) = BN(1,1) + ORBIAS(1) 00978 C 00979 CALL SPIN (BN,NERP,MNERP,DET,IDEXP) 00980 I = 0 00981 DO 35 IL = MI,MJ 00982 I = I + 1 00983 J = 0 00984 READ(ITF2) (B2(IB0) , IB0=1,KK) 00985 DO 35 JL = 1,KK 00986 J = J + 1 00987 35 BB(I,J) = B2(JL)*W 00988**2 MI = MJ 00989 CALL MATMPY (BB,FY,BTY,KK,L,1,MXOPP,MXOPP,MNUP,2) 00990 CALL MATMPY (BA,BB,BAT,NERP,L,KK,MXOPP,MXOPP,MNERP,2) 00991 CALL MATMPY (BAT,BN,BBT,KK,NERP,NERP,MNERP,MNERP,MNUP,2) 00992 DO 41 II = 1,NERP 00993 WRITE(ISF1) (BN(II,JJ) , JJ=1,NERP) 00994 41 WRITE(ISF2) (BAT(II,JK),JK=1,KK) 00995 CALL MATMPY (BBT,BAT,BCD,KK,NERP,KK,MNUP,MNERP,MNUP,1) 00996 CALL MATMPY (BB,BB,DN,KK,L,KK,MXOPP,MXOPP,MNUP,2) 00997 C 00998 DO 47 I1 = 1,KK,NERP 00999 J1 = I1 01000 J2 = J1 + 1 01001 J3 = J1 + 2 01002 GO TO (48,49,50) , NERP 01003 50 IF(DN(J3,J3).NE.0.0) DN(J3,J3) = DN(J3,J3) + ORBIAS(3) 01004 49 IF(DN(J2,J2).NE.0.0) DN(J2,J2) = DN(J2,J2) + ORBIAS(2) 01005 48 IF(DN(J1,J1).NE.0.0) DN(J1,J1) = DN(J1,J1) + ORBIAS(1) 01006 47 CONTINUE 01007 C 01008 CALL MATMPY (BBT,BTF,DFN,KK,NERP,1,MNUP,MNERP,MNUP,1) 01009 C 01010 DO 39 I = 1,KK 01011 BTY(I) = BTY(I) - DFN(I) 01012 IF(IK.EQ.1) UU(I) = BTY(I) 01013 DO 31 J = 1,KK 01014 31 DN(I,J) = DN(I,J) - BCD(I,J) 01015 IF(IK.EQ.1) WRITE(ITF4) (DN(I,J),J=1,KK) 01016 39 CONTINUE 01017 C 01018 IF(IK.EQ.1) GO TO 33 01019 END FILE ITF4 01020 REWIND ITF4 01021 DO 40 I = 1,KK 01022 UU(I) = UU(I) + BTY(I) 01023 READ(ITF4) (DN1(J),J=1,KK) 01024 DO 40 J = 1,KK 01025 40 DN(I,J) = DN(I,J) + DN1(J) 01026 REWIND ITF4 01027 DO 46 I = 1,KK 01028 46 WRITE(ITF4) (DN(I,J),J=1,KK) 01029 C 01030 33 CONTINUE 01031 WRITE(ISF3) (BTF(II) , II=1,NERP) 01032 32 CONTINUE 01033 C 01034 CALL SPIN (DN,KK,MNUP,DET,IDEXP) 01035 CALL MATMPY (DN,UU,XNS,KK,KK,1,MNUP,MNUP,MNUP,1) 01036 WRITE(IOUT,2020) 01037**2 2020 FORMAT('1',40X,'O R B I T A L B I A S E S S O L U T I O N ',///,01038 # 20X,'N/S REV.',5X,'ABS. BIAS',15X,'TILT',15X, 01039 #'SECOND ORDER EFFECT',/,36X,'(M)',16X,'(M/MIN)',18X,'(M/MIN**2)' 01040 # ,//) 01041 IORB = 0 01042 XPX = 0. 01043 DO 67 IPX = 1,KK,NERP 01044 IORB = IORB + 1 01045 IX1 = IPX 01046 IX2 = IPX + NERP - 1 01047 K1 = 0 01048**3 DO 69 KI = IX1,IX2 01049**3 K1 = K1 + 1 01050**3 YSN(K1) = XNS(KI) 01051**3 K2 = 0 01052**3 DO 69 KJ = IX1,KI 01053**3 K2 = K2 + 1 01054**3 BCOV(K1,K2) = DN(KI,KJ) 01055**3 IF(K2.NE.K1) BCOV(K2,K1) = BCOV(K1,K2) 01056**3 69 CONTINUE 01057**3 C 01058**3 IF(IWRITP.GT.5) 01059**3 # CALL WRITP (NS(IORB,3),NS(IORB,1),TMNS(IORB),NERP,YSN,BCOV) 01060**3 C 01061**3 WRITE(IOUT,2021) NS(IORB,3) , (XNS(I),I=IX1,IX2) 01062**2 2021 FORMAT(22X,I4,6X,E11.4,10X,E11.4,16X,E11.4,/) 01063 XPX = XPX + XNS(IPX)**2 * ORBIAS(1) 01064 IF(NERP.LT.2) GO TO 67 01065 XPX = XPX + XNS(IPX+1)**2 * ORBIAS(2) 01066 IF(NERP.LT.3) GO TO 67 01067 XPX = XPX + XNS(IPX+2)**2 * ORBIAS(3) 01068 67 CONTINUE 01069 C 01070 END FILE ISF1 01071 END FILE ISF2 01072 END FILE ISF3 01073 REWIND ISF1 01074 REWIND ISF2 01075 REWIND ISF3 01076 C 01077 CALL CLEAR (XSN,MNUP) 01078 ISUM = 0 01079 WRITE(IOUT,2022) 01080**2 2022 FORMAT(//,20X,'S/N REV.',5X,'ABS. BIAS',15X,'TILT',15X, 01081 #'SECOND ORDER EFFECT',/,36X,'(M)',16X,'(M/MIN)',18X,'(M/MIN**2)' 01082 # ,//) 01083 DO 42 I = 1,NSNP 01084 I1 = (I-1)*NERP+1 01085 I2 = I*NERP 01086 DO 43 II = 1,NERP 01087 READ(ISF1) (BN(II,JJ) , JJ=1,NERP) 01088 READ(ISF2) (BAT(II,JK),JK=1,KK) 01089 43 CONTINUE 01090 READ(ISF3) (BTF(II) , II=1,NERP) 01091 CALL MATMPY (BAT,XNS,DF,NERP,KK,1,MNERP,MNUP,MNERP,1) 01092 DO 44 M=1,NERP 01093 44 DG(M) = BTF(M) - DF(M) 01094 CALL MATMPY (BN,DG,YSN,NERP,NERP,1,MNERP,MNERP,MNERP,1) 01095 WRITE(IOUT,2021) SN(I,3) , (YSN(IORB),IORB=1,NERP) 01096**2 IF(IWRITP.GT.5) CALL WRITP (SN(I,3),SN(I,1),TMSN(I),NERP,YSN,BN) 01097**3 MI = 0 01098 DO 45 M = I1,I2 01099 MI = MI + 1 01100 45 XSN(M) = XSN(M) + YSN(MI) 01101 C 01102 ISUM = ISUM + SN(I,2) 01103 42 CONTINUE 01104 DO 68 IPX = 1,KL,NERP 01105 XPX = XPX + XSN(IPX)**2 * ORBIAS(1) 01106 IF(NERP.LT.2) GO TO 68 01107 XPX = XPX + XSN(IPX+1)**2 * ORBIAS(2) 01108 IF(NERP.LT.3) GO TO 68 01109 XPX = XPX + XSN(IPX+2)**2 * ORBIAS(3) 01110 68 CONTINUE 01111 REWIND ITF1 01112 REWIND ITF2 01113 C 01114 CALL CLEAR (ASN,MXOBS) 01115 MI = 0 01116 MJ = 0 01117 DO 51 IK = 1,NSNP 01118 I1 = (IK-1)*NERP+1 01119 I2 = IK*NERP 01120 L = SN(IK,2) 01121 MI = MI + 1 01122 MJ = MJ + L 01123 I = 0 01124 C 01125 DO 52 IL = MI,MJ 01126 I = I + 1 01127 J = 0 01128 READ(ITF1) (B1(IB) , IB=1,NERP) 01129 C 01130 DO 53 JL = 1,KL 01131 J = J + 1 01132 53 BA(I,J) = B1(JL) 01133 C 01134 M = 0 01135 READ(ITF2) (B2(IB0),IB0=1,KK) 01136 DO 54 ML = 1,KK 01137 M = M + 1 01138 54 BB(I,M) = B2(ML) 01139 C 01140 52 CONTINUE 01141 C 01142 MK = 0 01143 DO 56 M = I1,I2 01144 MK = MK + 1 01145 56 B1(MK) = XSN(M) 01146 CALL MATMPY (BA,B1,FY,L,NERP,1,MXOPP,MNERP,MXOPP,1) 01147 CALL MATMPY (BB,XNS,DG,L,KK,1,MXOPP,MNUP,MNUP,1) 01148 C 01149 I = 0 01150 DO 55 IL = MI,MJ 01151 I = I + 1 01152 ASN(IL) = FY(I) 01153 55 ANS(IL) = DG(I) 01154 MI = MJ 01155 51 CONTINUE 01156 REWIND ITF2 01157 RETURN 01158 END 01159 SUBROUTINE SPIN(Q,N,MM,DET,IDEXP) 01160 C 01161 C SUBROUTINE SPIN IS A MATRIX INVERSION ROUTINE FOR SYMMETRIC 01162 C POSITIVE-DEFINITE MATRICES. THE MATRIX INVERTED IS THE UPPER 01163 C N BY N PORTION OF THE MATRIX Q WHICH IS DIMENSIONED MM BY MM 01164 C IN THE CALLING ROUTINE. 01165 C 01166 C INPUT: 01167 C Q - THE MATRIX DIMENSIONED MM BY MM WHICH CONTAINS THE 01168 C MATRIX TO BE INVERTED. 01169 C 01170 C N - THE DIMENSION OF THE ACTUAL PART( UPPER LEFT CORNER) 01171 C OF Q WHICH IS TO BE INVERTED. ( N MAY BE EQUAL BUT MUST 01172 C NOT BE LARGER THAN MM) . 01173 C MM- DIMENSIONED SIZE OF Q IN THE CALLING ROUTINE. 01174 C 01175 C 01176 C OUTPUT: 01177 C 01178 C Q - THE UPPER LEFT N BY N PORTION CONTAINS THE INVERSE OF 01179 C THE INPUT UPPER LEFT N BY N PORTION. 01180 C 01181 C DET - THE NON-EXPONENT PORTION OF THE DETERMINANT OF THE 01182 C INPUT N BY N (UPPER LEFT PORTION OF Q) MATRIX. SEE 01183 C IDEXP BELOW. 01184 C 01185 C IDEXP - THE EXPONENT (OF 10) PART OF THE DETERMINANT DESCRIBED 01186 C UNDER DET ABOVE. THUS THE DETERMINANT IS RETURNED IN 01187 C TWO PARTS CORRESPONDING TO 01188 C DETERMINANT = DET * 10 ** IDEXP . 01189 C THIS IS DONE TO AVOID UNDER OR OVERFLOW IN THE 01190 C COMPUTATION OF THE DETERMINANT. TO PRINT THE DETERM- 01191 C INANT THE USER SHOULD PRINT BOTH NUMBERS AS FOLLOWS; 01192 C (FOR EXAMPLE) 01193 C PRINT 10,DET,IDEXP 01194 C 10 FORMAT(' ','DETERMINANT= ',F7.4,'D',I4) 01195 C 01196 C R.R. STEEVES 01197 C SEPT., 1979 01198 C 01199 C 01200 REAL*8 Q(MM,MM),DET,SUM,DSQRT,DABS,RPART,APART 01201 DET=0.D0 01202 DO 4 J=1,N 01203 DO 4 I=1,J 01204 IF(I.EQ.1) GO TO 2 01205 M=I-1 01206 SUM=0.0D0 01207 DO 1 K=1,M 01208 1 SUM=SUM+Q(K,I)*Q(K,J) 01209 Q(I,J)=Q(I,J)-SUM 01210 2 IF(I.EQ.J) GO TO 3 01211 Q(I,J)=Q(I,J)/Q(I,I) 01212 GO TO 4 01213 3 CONTINUE 01214 DET=DET+DLOG10(Q(I,I)) 01215 Q(I,I)=DSQRT(Q(I,I)) 01216 4 CONTINUE 01217 IDEXP=DET 01218 RPART=DET-IDEXP 01219 APART=DABS(RPART) 01220 IF(APART.LT.1.D-20) DET=1.D0 01221 IF(APART.LT.1.D-20) GO TO 10 01222 DET=10**RPART 01223 10 CONTINUE 01224 DO 7 J=1,N 01225 DO 7 I=1,J 01226 IF(I.LT.J) GO TO 5 01227 Q(J,J)=1.0D0/Q(J,J) 01228 GO TO 7 01229 5 SUM=0.0D0 01230 M=J-1 01231 DO 6 K=I,M 01232 6 SUM=SUM-Q(I,K)*Q(K,J) 01233 Q(I,J)=SUM/Q(J,J) 01234 7 CONTINUE 01235 DO 9 J=1,N 01236 DO 9 I=1,J 01237 SUM=0.0D0 01238 DO 8 K=J,N 01239 8 SUM=SUM+Q(I,K)*Q(J,K) 01240 Q(I,J)=SUM 01241 IF(I.EQ.J) GO TO 9 01242 Q(J,I)=SUM 01243 9 CONTINUE 01244 RETURN 01245 END 01246 INTEGER FUNCTION SELECT (IREV1,IREV2,SLCT,N) 01247 DIMENSION SLCT(N) 01248 SELECT = 1 01249 DO 10 I = 1,N 01250 ISELCT = INT(SLCT(I)) 01251 LREV1 = ISELCT/10000 01252 LREV2 = MOD(ISELCT,10000) 01253 IF(IREV1.EQ.LREV1 .AND. IREV2.EQ.LREV2) GO TO 2 01254 10 CONTINUE 01255 RETURN 01256 2 SELECT = 2 01257 RETURN 01258 END 01259 SUBROUTINE PLOT(WINT,HVEC) 01260 C***********************************************************************01261 C* 01262 C* PLOT PLOTS THE STANDARD NORMAL CURVE OVERLAYED WITH THE HISTOGRAM OF01263 C* STANDARD RESIDUALS. 01264 C* 01265 C* 01266 C* INPUT: 01267 C* WINT- NUMBER OF HISTOGRAM INTERVALS 01268 C* HVEC- VECTOR CONTAINING THE NUMBER OF RESIDUALS IN EACH HISTOG01269 C* INTERVAL 01270 C* 01271 C* 01272 C* WRITTEN BY: 01273 C* LAURIE PACH, JULY, 1978 01274 C* 01275 C***********************************************************************01276 INTEGER WINT,NVEC,RVEC,HVEC,STAR,HLINE,VLINE,PLOTV,PLOTL,PLOTH 01277 DIMENSION PLOTV(110),NVEC(53),PLOTL(110),PLOTH(110),RVEC(20), 01278 @ HVEC(20) 01279 DATA STAR,HLINE,VLINE/'.','-','|'/ 01280 DATA NVEC/22*0,4*1,2,2,3,4,4,5,6,8,9,10,12,14,16,17,19,21,23,25, 01281 @27,28,29,31,31,32,32,32,31/ 01282 MAX=50 01283 A=.4 01284 K=2 01285 I=52 01286 NUMR=0 01287 KK=1 01288 NVAL=32 01289 NFLG=0 01290 IF(WINT.EQ.20)INT=5 01291 IF(WINT.EQ.10.OR.WINT.EQ.9)INT=10 01292 IF(WINT.EQ.2)INT=50 01293 IF(WINT.EQ.3)INT=30 01294 IF(WINT.EQ.5)INT=20 01295 IF(WINT.EQ.4)INT=25 01296 WIDF = 10./FLOAT(INT) 01297**5 CALL PBLANK(PLOTV) 01298 CALL PBLANK(PLOTH) 01299 CALL PBLANK(PLOTL) 01300 DO 2 JJ=1,WINT 01301 NUMR=NUMR+HVEC(JJ) 01302 2 CONTINUE 01303 DO 33 JJ=1,WINT 01304 RVEC(JJ)=(80*HVEC(JJ)/NUMR)*WIDF+.5 01305**5 33 CONTINUE 01306 PRINT103 01307 DO 29 JJ=1,50 01308 DO 28 N=1,WINT 01309 IF(RVEC(N).GE.MAX)GOTO19 01310 28 CONTINUE 01311 PRINT104 01312 IF(N.EQ.(WINT+1))N=WINT 01313 IF(MAX.EQ.32)GOTO21 01314 MAX=MAX-1 01315 29 CONTINUE 01316 19 RVEC(N)=0 01317 MM=(INT*(N-1))+1 01318 IF(WINT.EQ.9.OR.WINT.EQ.3)MM=(INT*(N-1))+6 01319 PLOTL(MM)=VLINE 01320 III=INT-1 01321 DO 12 JJ=1,III 01322 PLOTH(JJ+MM)=HLINE 01323 12 CONTINUE 01324 PLOTL(MM+INT)=VLINE 01325 DO 32 N=1,WINT 01326 IF(RVEC(N).EQ.MAX.OR.RVEC(N).GT.50)GOTO19 01327 32 CONTINUE 01328 IF(NFLG.EQ.1)GOTO25 01329 PRINT101,(PLOTH(L),L=1,110) 01330 CALL PBLANK(PLOTH) 01331 CALL LPRNT(NVAL,MAX,PLOTL,RVEC,WINT,N,KK) 01332 IF(RVEC(N).EQ.MAX.AND.MAX.NE.32)GOTO19 01333 21 DO 31 L=2,100 01334 PLOTV(I)=STAR 01335 IF(NVEC(I).NE.NVEC(I-1))GOTO4 01336 K=K+1 01337 I=I-1 01338 31 CONTINUE 01339 4 I=I-1 01340 NFLG=1 01341 IF(RVEC(N).EQ.MAX)GOTO19 01342 25 DO 26 L=1,110 01343 IF(PLOTH(L).EQ.HLINE)GOTO36 01344 26 CONTINUE 01345 GOTO38 01346 36 DO 39 L=1,110 01347 IF(PLOTH(L).EQ.HLINE)PLOTV(L)=PLOTH(L) 01348 39 CONTINUE 01349 38 CALL PBLANK(PLOTH) 01350 11 PRINT102,(PLOTV(L),L=1,110) 01351 IF(MAX.EQ.32)PRINT113,A 01352 IF(I.EQ.1)GOTO20 01353 NVAL=NVEC(I) 01354 IF(NVAL.EQ.0)GOTO20 01355 CALL PBLANK(PLOTV) 01356 DO 37 L=1,110 01357 IF(PLOTL(L).EQ.VLINE)PLOTV(L)=PLOTL(L) 01358 37 CONTINUE 01359 CALL LPRNT(NVAL,MAX,PLOTV,RVEC,WINT,N,KK) 01360 CALL PBLANK(PLOTV) 01361 DO 3 L=2,100 01362 PLOTV(I)=STAR 01363 PLOTV(I+K)=STAR 01364 K=K+2 01365 IF(I.EQ.1)GOTO11 01366 IF(NVEC(I).NE.NVEC(I-1))GOTO4 01367 I=I-1 01368 3 CONTINUE 01369 GOTO4 01370 20 PRINT111 01371 PRINT112 01372 IF(WINT.EQ.20)PRINT107,(HVEC(L),L=1,20) 01373 IF(WINT.EQ.10)PRINT106,(HVEC(L),L=1,10) 01374 IF(WINT.EQ.4)PRINT108,(HVEC(L),L=1,4) 01375 IF(WINT.EQ.2)PRINT109,(HVEC(L),L=1,2) 01376 IF(WINT.EQ.9)PRINT116,(HVEC(L),L=1,9) 01377 IF(WINT.EQ.5)PRINT114,(HVEC(L),L=1,5) 01378 IF(WINT.EQ.3)PRINT115,(HVEC(L),L=1,3) 01379 101 FORMAT(' ',6X,110A1) 01380 102 FORMAT('+',6X,110A1) 01381 103 FORMAT('1') 01382 104 FORMAT(' ') 01383 106 FORMAT(' ',8X,10(I4,6X),/) 01384 108 FORMAT(' ',15X,I4,3(21X,I4),/) 01385 107 FORMAT(' ',7X,20(I4,1X),/) 01386 109 FORMAT(' ',34X,I4,36X,I4,/) 01387 111 FORMAT(' ',6X,20('|----'),'|') 01388 112 FORMAT(' ',5X,'-10',7X,'- 8',7X,'- 6',7X,'-4',8X,'-2',9X,'0',9X, 01389 # '2',9X,'4',8X,' 6',8X,' 8',8X,'10',/) 01390 113 FORMAT('+',3X,F3.1,'-') 01391 114 FORMAT(' ',13X,I4,4(16X,I4),/) 01392 115 FORMAT(' ',23X,3(I4,26X),/) 01393 116 FORMAT(' ',14X,9(I4,6X),/) 01394 RETURN 01395 END 01396 SUBROUTINE WRITP (IREV,IFILE,TM,NERP,COEF,COV) 01397**3 C 01398**4 C NAME WRITP 01399**4 C 01400**4 C TYPE SUBROUTINE 01401**4 C 01402**4 C PURPOSE WRITE ON TAPE A DATA SET PER SATELLITE ARC 01403**4 C CONTAINING INPUT DATA , RECOVERED ORBITAL ERRORS 01404**4 C AND ADJUSTED SEA SURFACE HEIGHTS 01405**4 C 01406**4 C I/O DEVICES ITAPE - OUTPUT MAGNETIC TAPE LOGICAL UNIT (UNIT 9) 01407**4 C 01408**4 C CALLING ROUTINES YRMDAY , ELAPST , TIDEIN , HHMMSS 01409**4 C 01410**4 C 01411**3 DEFINE FILE 97(400,55,L,LOCREC) , 98(25000,90,L,LPRSUM) 01412**3 DOUBLE PRECISION FODAY 01413**3 C 01414**3 COMMON NEDR , IREC , IFORM , LSMALT , LSIG , LAT , LONG , 01415**3 # LTROP , LGEOID , LTIDE , LSEAH , IRES , LH 01416**3 COMMON /TIDEM/ C , COVAR , CA , COA, CB , COB , 01417**3 # YLAT , YLON , ALON , XLAT , XLON , NCONS , 01418**3 # L , N , T , HTIDE 01419**3 COMMON /BLOCK/ WFREQ(4) , FK(4) 01420**3 COMMON /AREA/ XLAT1,XLAT2,XLON1,XLON2 /CUTOFF/ SSHMN , SSHMX 01421**3 EQUIVALENCE (NEDR,IRCD(1)) 01422**3 DIMENSION C(4,16,2) , COVAR(4,16,2) , CA(4,16,2) , COA(4,16,2) , 01423**3 # CB(4,16,2) , COB(4,16,2) , YLAT(4,2) , YLON(4,2) , 01424**3 # ALON(4) 01425**3 DIMENSION IRCD(13) , COEF(3) , B1(3) , BT(3) , COV(3,3) 01426**3 DATA ITAPE / 9 / 01427**3 C 01428**3 LOCREC = IFILE 01429**3 READ(97'LOCREC,1002,END=901) LFILE,LSREC,LEREC,LREV, 01430**3 # IDATE,NX,LEQXT,XLNODE 01431**3 IF(LREV.NE.IREV) GO TO 999 01432**3 C 01433**3 CALL YRMDAY (IDATE,IYR,IMON,IDAY) 01434**3 CALL HHMMSS (TM,IHR,IMIN,SEC) 01435**3 CALL ELAPST (244,IYR,IMON,IDAY,IHR,IMIN,SEC,T) 01436**3 C 01437**3 I1 = LSREC 01438**3 I2 = LEREC 01439**3 IND = 0 01440**3 C 01441**3 DO 1 I = I1,I2 01442**3 LOCR = I 01443**3 READ(98'LOCR,1003,END=901) (IRCD(J),J=1,3) , FODAY , 01444**3 # (IRCD(J),J=4,13) 01445**3 XLAT = FLOAT(LAT)/1.E4 01446**3 XLON = FLOAT(LONG)/1.E4 01447**3 C 01448**3 C CHECK WHETHER ALTIMETRY MEASUREMENT IS OVER SPECIFIED AREA 01449**3 C 01450**3 IF(XLAT.LT.XLAT1 .OR. XLAT.GT.XLAT2 .OR. 01451**3 # XLON.LT.XLON1 .OR. XLON.GT.XLON2) GO TO 1 01452**3 C 01453**3 SSH = FLOAT(LSEAH)/1.E2 01454**3 C 01455**3 C PASS SSH THROUGTH THE INNER AND OUTER CUTOFF WINDOWS 01456**3 C 01457**3 IF(SSH.LT.SSHMN .OR. SSH.GT.SSHMX) GO TO 1 01458**3 IND = IND + 1 01459**3 C 01460**3 C CORRECT MEASUREMENT FOR TIDAL EFFECT 01461**3 C 01462**3 CALL TIDEIN 01463**3 C 01464**3 SSHC = (FLOAT(LSEAH)-HTIDE)/1.E2 01465**3 C 01466**3 C CORRECT FOR ORBIT BIAS 01467**3 C 01468**3 DT = (FODAY-TM)/6.E7 01469**3 GO TO (2,3,4) , NERP 01470**3 4 B1(3) = DT**2 01471**3 3 B1(2) = DT 01472**3 2 B1(1) = 1. 01473**3 C 01474**3 CBIAS = 0.0 01475**3 C 01476**3 C 01477**3 DO 5 K = 1,NERP 01478**3 IF(IND-1) 12,13,12 01479**3 12 CBIAS = CBIAS + COEF(K)*B1(K)/1.E6 01480**3 GO TO 5 01481**3 13 CBIAS = CBIAS + COEF(K)*B1(K) 01482**3 5 CONTINUE 01483**3 C 01484**3 DO 6 JJ = 1,NERP 01485**3 SUM = 0. 01486**3 DO 7 K = 1,NERP 01487**3 IF(IND-1) 14,15,14 01488**3 14 SUM = SUM + COV(JJ,K)*B1(K)/1.E6 01489**3 GO TO 7 01490**3 15 SUM = SUM + COV(JJ,K)*B1(K) 01491**3 7 CONTINUE 01492**3 6 BT(JJ) = SUM 01493**3 SUM = 0. 01494**3 DO 8 J = 1,NERP 01495**3 8 SUM = SUM + B1(J)*BT(J) 01496**3 VAR = SUM 01497**3 C 01498**3 SSHC1 = SSHC - CBIAS 01499**3 C 01500**3 ITIDE = INT(HTIDE) 01501**3 IBIAS = INT(CBIAS*100.) 01502**3 IBIASD = INT(SQRT(VAR)*100.) 01503**3 ISEAH = INT(SSHC1*100.) 01504**3 ISIG = INT( SQRT((FLOAT(LSIG)/1.E2)**2 + VAR)*100.) 01505**3 C 01506**3 IF(IND.NE.1) GO TO 11 01507**3 DO 10 II = 1,NERP 01508**3 COEF(II) = COEF(II)*1.E6 01509**3 DO 10 JJ = 1,II 01510**3 COV(II,JJ) = COV(II,JJ)*1.E6 01511**3 10 CONTINUE 01512**3 11 CONTINUE 01513**3 C 01514**3 C 01515**3 WRITE(ITAPE,1004) IREV,IFILE,NEDR,IFORM,LSMALT,LSIG,FODAY, 01516**3 # LAT,LONG,LTROP,LGEOID,LTIDE,LSEAH,IRES,LH,NERP,TM, 01517**3 # (COEF(II),II=1,3) , ((COV(II,JJ),JJ=1,II),II=1,3) , 01518**3 # ITIDE,IBIAS,IBIASD,ISEAH,ISIG 01519**3 C 01520**3 1 CONTINUE 01521**3 GO TO 16 01522**5 999 WRITE(6,1000) LREV , IREV 01523**3 1000 FORMAT(5X,'***** WARNING : INCONSISTENT REV NUMBERS - LREV = ', 01524**3 # I4,' IREV = ',I4,' *****') 01525**3 1001 FORMAT(2I5,2F16.0,7F10.4) 01526**3 1002 FORMAT(7I6,F13.6) 01527**3 1003 FORMAT(2I4,I3,F12.0,I9,I6,2I8,6I6) 01528**3 1004 FORMAT(3I5,I2,I9,I6,F12.0,2I8,6I6,I2,F12.0,9F10.0,5I6) 01529**3 900 STOP 01530**3 901 CONTINUE 01531**3 16 CONTINUE 01532**5 DO 17 II = 1,NERP 01533**5 COEF(II) = COEF(II)/1.E6 01534**5 DO 17 JJ = 1,II 01535**5 COV(II,JJ) = COV(II,JJ)/1.E6 01536**5 17 CONTINUE 01537**5 END FILE ITAPE 01538**3 RETURN 01539**3 END 01540**3