C * * * * * ETIDAN * * * * * 00001 C 00002 C THIS PROGRAM IS DESIGNED TO ANALYSE OBSERVED TILT DATA SEEKING UP TO00003 C 6 PERIODIC CONSTITUENTS,TEMPERATURE AND PRESSURE DEPENDENCE COEFFI- 00004 C CIENTS AND , IF DESIRED, PRECIPITATION DEPENDENCE COEFFICIENT. 00005 C IT PLOTS THE OBSERVED DATA,THE RESIDUALS AND SPECIFIED SPECTRUM 00006 C 00007 IMPLICIT REAL * 8(A-H,O-Z) 00008 REAL F2,SUM,VALMIN,VALMAX,SCALE DIMENSION COVC(65,65),A(65,65),DAT(55),PER(6),COEF(65),SEA(65), 00009 &T(4000),IIPR(41),PREC(41) DIMENSION ILT(12) 00011 DIMENSION LAG(3),XBACK(25),FX(25),FPHI(25) DIMENSION IPLOT(100) DIMENSION F(4000) INTEGER RR,S 00013 COMMON E(3,6500),TM 00014*17 COMMON SING 00015 COMMON /TILT/F2(6500) COMMON /MEANS/SUM COMMON /NUMBR/N,SCALE,IRESID FLOAT(II)=DFLOAT(II) RR = 5 00016*16 IPR=6 C C READ FLAGS FOR PLOTTING OF TILT DATA,TEMPERATURE,PRESSURE, C AND RESIDUAL. (0 - DONT PLOT, 1 - PLOT) C READ(5,200)ITILT,ITEMP,IPRESS,IRESID 200 FORMAT(4I2) C C READ SCALE= # OF DAYS PER INCH IN PLOTS. C READ(5,201)SCALE 201 FORMAT(F3.1) C 00017 C READ REQUIRED PHASE LAG (IN HOURS), FOR 00018 C TEMPERATURE, PRESSURE AND PERCIPITATION 00019 C 00020 READ( 5,66)(LAG(I),I=1,3) 66 FORMAT (3I3) 00022 DO 37 I=1,25,5 KT = I + 4 READ ( 5,82) (XBACK(L),L=I,KT) 37 CONTINUE C 00023 C READ STATION CONTROL CARD ( N = NUMBER OF HALF-DAYS,IY = STARTING 00024 C DAY,INO = STATION NUMBER) 00025 C 00026 READ(RR,1) N,IY,INO 00027 C 00028 C READ CONTROL CARD SPECIFYING THE REQUIRED SPECTRUM 00029 C FOR DESCRIPTION OF PARAMETERS SEE SUBROUTINE SPECEQ 00030 C 00031 READ (RR,4) PL,PS,NW 00032 IB = 1 00033 C 00034 C READ FORCED PERIODS IN HOURS (UP TO 6) FOLLOWED BY A BLANK CARD 00035 C 00036 NPER=0 00037 DO 5 I = 1,7 00038 READ(RR,11) PPER 00039 PER(I) = PPER 00040 IF(PER(I).EQ.0.0) GO TO 36 NPER = NPER + 1 00042 5 CONTINUE 00043 C 00044 C READ MAGNIFICATION FACTORS FOR TILT, TEMPERATURE, PRESSURE AND 00045 C PRECIPITATION (IF REQUIRED) AND VALUES FOR ZERO READING OF TEMPERATURE00046 C AND PRESSURE 00047 C 00048 36 READ(RR,9) R1,R2,R3,R4,Q2,Q3 00049 C 00050 C READ HOURLY TILT VALUES.THE DATA SET MUST START WITH HOUR 0 AND 00051 C THE SEQUENCE MUST INCLUDE BLANK CARDS FOR GAPS. ANY GAP IS AUTOMA- 00052 C TICALLY TAKEN AS BEGINNING OF A NEW DATUM. MAX. NUMBER OF GAPS = 55 00053 C 00054 SUM=0. DIVIDE=1.D-35 VALMIN=1./DIVIDE VALMAX=-VALMIN NM = 0 00055 KK = 0 00056 J = 1 00057 PREV = 1.0 00058 DO 2 I = 1,N 00059 READ (RR,12) (ILT(K), K = 1,12) 00060 DO 13 K = 1,12 00061 KK = KK + 1 00062 F2(KK)=ILT(K)*R1 IF (ILT(K).NE.0) GO TO 14 00063 IF (PREV.EQ.0.0) GO TO 13 00064 PREV = 0.0 00065 J = J + 1 00066 DAT(J) = KK 00067 GO TO 13 00068 14 PREV = 1.0 00069 NM = NM + 1 00070 F(NM) = ILT(K) * R1 00071 SUM=SUM+F(NM) IF(F(NM).LT.VALMIN) VALMIN=F(NM) IF(F(NM).GT.VALMAX) VALMAX=F(NM) T(NM) = FLOAT(KK) 00072 13 CONTINUE 00073 2 CONTINUE 00074 IF(ILT(12).EQ.0)J=J-1 PRINT 101,VALMIN,VALMAX SUM=SUM/DFLOAT(NM) NN=0.5*N+0.99 IF(ITILT.EQ.1) CALL VALPLT(IY,NN,VALMIN,VALMAX,SCALE,3) NDAT = J 00075 DAT(1) = T(1) 00076 C 00077 C PRINT HEADING 00078 C 00079 NY = 0.5*N + 0.99 00080 NLY = NY + IY - 1 00081 PRINT 7,INO 00082 PRINT 6, NM, NY, IY, NLY 00083 TM = FLOAT(NLY - IY) * 24.0 00084*17 C 00085 C CLEAR TEMPERATURE AND PRESSURE ARRAYS 00086 C 00087 NN = 12*N 00088 DO 24 I=1,2 00089 IJ=12*N+LAG(I) 00090 DO 24 J=1,IJ 00091 24 E(I,J)=0.0 00092 C 00093 C READ AND PRINT HOURLY TEMPERATURES 00094 C 00095 VALMIN=1./DIVIDE VALMAX=-VALMIN KK = 0 00096 NR = 0 00097 SUM = 0.0 00098 DO 8 I = 1,N 00099 READ (RR,12) (ILT(K), K = 1,12) 00100 DO 8 K=1,12 00101 KK = KK+1 00102 IF(ILT(K).EQ.0.) GO TO 8 IF(ILT(K).NE.0) E(1,KK+LAG(1))=ILT(K)*R2-Q2 00103 IF(E(1,KK+LAG(1)).LT.VALMIN) VALMIN=E(1,KK+LAG(1)) IF(E(1,KK+LAG(1)).GT.VALMAX) VALMAX=E(1,KK+LAG(1)) 8 CONTINUE 00104 C 00105**4 C 00106*20 IF(LAG(1).EQ.0) GOTO 98 00107**9 LAST = LAG(1) 00108**4 DO 40 J = 1,NN 00109**4 IF (E(1,J).EQ.0.0) GOTO 40 00110**4 IF (E(1,J-1).NE.0.0) GOTO 40 00111**4 ISTOP = J-1 00112**4 CALL LSAPRO (LAST,ISTOP,IY,XBACK) 40 CONTINUE 00113**4 98 CONTINUE 00114**9 KK=0 00115 DO 99 I=1,NN 00116 KK=KK+1 00117 IF (E(1,KK).EQ.0.0) GOTO 99 00118 SUM = SUM + E(1,KK) 00119 NR = NR + 1 00120 99 CONTINUE 00121 PRINT 101, VALMIN,VALMAX PRINT 25, LAG(1) 00122 CALL FIPRIN (N,1,IY) 00123 SUM = SUM/DFLOAT(NR) 00124 NN=0.5*N+0.99 IF(ITEMP.EQ.1) CALL VALPLT(IY,NN,VALMIN,VALMAX,SCALE,1) DO 19 I = 1,NN 00125 19 E(1,I) = E(1,I) - SUM 00126 PRINT 22, SUM 00127 C 00128 C READ AND PRINT HOURLY PRESSURES 00129 C 00130 VALMIN=1./DIVIDE VALMAX=-VALMIN KK = 0 00131 NR = 0 00132 SUM = 0.0 00133 DO 10 I = 1,N 00134 READ (RR,12) (ILT(K), K = 1,12) 00135 DO 10 K=1,12 00136 KK = KK+1 00137 IF(ILT(K).EQ.0.) GO TO 10 IF (ILT(K).NE.0) E(2,KK) = ILT(K)*R3+Q3 00138**2 IF(E(2,KK).LT.VALMIN) VALMIN=E(2,KK) IF(E(2,KK).GT.VALMAX) VALMAX=E(2,KK) IF (E(2,KK).EQ.0.0) GOTO 10 00139**2 SUM = SUM + E(2,KK) 00140**2 NR = NR + 1 00141 10 CONTINUE 00142 PRINT 101, VALMIN,VALMAX CALL FIPRIN (N,2,IY) 00143 SUM = SUM/DFLOAT (NR) 00144 NN=0.5*N+0.99 IF(IPRESS.EQ.1) CALL VALPLT(IY,NN,VALMIN,VALMAX,SCALE,2) DO 20 I = 1,NN 00145 20 E(2,I) = E(2,I) - SUM 00146 PRINT 22, SUM 00147 C 00148 C READ PRECIPITATION FOR EVERY 6 HOURS. FIRST CARD MUST CONTAIN AT LEAST00149 C ONE NUMBER DIFFERENT FROM ZERO. PRINT HOURLY INTERPOLATED 00150 C PRECIPITATION. IF NO PRECIPITATION IS SUPPLIED INSERT A NLANK CARD 00151 C 00152 NB=1 NNN = 0.05*N + 0.99 00154 L = 1 00155 SUM = 0.0 00156 DO 300 I=1,NNN 00157 READ (RR,301)(IIPR(J), J = 1,40) 00158 IF (I.NE.1) GO TO 16 00159 ISU = 0 00160 DO 15 KK = 1,40 00161 15 ISU = ISU + IIPR(KK) 00162 IF(ISU.NE.0) GO TO 16 00163 NB = NB-1 00164 GO TO 17 00165 16 K = 1 00166 DO 300 J = 1,40 00167 PREC (K) = IIPR(J)*R4 00168 IF(I.EQ.1.AND.J.EQ.1.AND.K.EQ.1) GO TO 303 00169 IF(K.EQ.1) PRECO = PREC(40) 00170 DO 305 JJ = 1,6 00171 IF(K.EQ.1) GO TO 304 00172 E(3,L) = (PREC(K) - 0.166667 * (6-JJ)*(PREC(K) - PREC(K-1)))/6.0 00173**2 GO TO 302 00174 304 E(3,L) = (PREC(K) - 0.166667 * (6-JJ) * (PREC(K) - PRECO))/6.0 00175**2 302 SUM = SUM + E(3,L) 00176**2 305 L = L + 1 00177 303 K = K + 1 00178 300 CONTINUE 00179 CALL FIPRIN (N,3,IY) 00180 SUM = SUM/DFLOAT(NN) 00181 DO 23 I = 1,NN 00182 23 E(3,I) = E(3,I) - SUM 00183 PRINT 22, SUM 00184 17 CONTINUE 00185 C 00186 C COMPUTATION OF SPECTRUM AND BEST FITTING COEFFFICIENTS 00187 C 00188 LT = 1 00189 M = NDAT + LT + 2*NPER +NB 00190 CALL SPECEQ (F,T,NM,PL,PS,NW,IB,DAT,NDAT,LT,PER,NPER,NB,COEF,RN, 00191 &A,IY) 00192 C 00193 C COMPUTE THE STANDARD DEVIATIONS OF COEFFICIENTS 00194 C 00195 69 RN=RN/(NM-M) 00196 DO 55 I=1,M 00197 COEF(I) = COEF(I) 00198 DO 55 J=1,M 00199 55 COVC(I,J)=A(I,J)*RN 00200 DO 33 L=1,M 00201 33 SEA(L)=DSQRT(COVC(L,L)) 00202 C 00203 C DATUM BIASES 00204 C 00205 DO 21 I = 1,NDAT 00206 21 PRINT 85,COEF(I),SEA(I) 00207 C 00208 C LINEAR TREND 00209 C 00210 NL = NDAT + 1 00211 SELT=SEA(NL)*24.0 00212 COLT=COEF(NL)*24.0 00213 IF(LT.EQ.1) PRINT 86,COLT,SELT 00214 C 00215 C 00217 C PERIODIC CONSTITUENTS 00216 PRINT 113 00218 IY = IY - 1 00219 PRINT 116,IY 00220 J = NDAT + LT + 1 00221 K = J + 2*NPER - 1 00222 DO 110 I = J,K,2 00224 L = 0 00223 CALL AMPLAG (COEF(I), COEF(I+1), COVC(I,I), COVC(I+1,I+1), COVC 00225 &(I+1,I),R,P,SR,SP) 00226 L = L + 1 00227 PRINT 112, PER(L),R,SR,P,SP 00228 110 CONTINUE 00229 C 00230 C OTHER CONSTITUENTS 00231 C 00232 PRINT 80 00233 K=K+1 00234 PRINT 83, COEF(K), SEA(K) 00235 K=K+1 00236 PRINT 84, COEF(K), SEA(K) 00237 IF (NB.EQ.2) GO TO 18 00238 K = K + 1 00239 PRINT 87, COEF(K), SEA(K) 00240 18 CONTINUE 00241 C 00242 C HIGH CORRELATINS BETWEEN COEFFICIENTS 00243 C 00244 PRINT 89 00245 MM = M - 1 00246 DO 120 I = 1,MM 00247 K = I + 1 00248 DO 120 J = K,M 00249 Q = DSQRT (COVC(I,I) * COVC (J,J)) 00250 R = COVC (I,J)/Q 00251 Q = DABS (R) 00252 IF (Q.LT.0.5)GO TO 120 00253 PRINT 88, I,J,R 00254 120 CONTINUE 00255 STOP 00256 C 00257 C FORMATS 00258 C 00259 1 FORMAT (3I4) 00260 4 FORMAT(F6.0,F4.0,I3) 00261*16 6 FORMAT (///2X,I4,'HOURLY MEANS FOR', I4,'DAYS FROM',I4, 'TO', 00262 &I4) 00263 7 FORMAT ('1', 1X, 'STATION', 3X, I4) 00264 9 FORMAT(6F8.5) 11 FORMAT(F10.6) 12 FORMAT (8X,12I6) 00267 22 FORMAT (//8X, 'MEAN =', F8.3 ///) 00268 25 FORMAT (///,10X,'INTRODUCED PHASE LAG (TEMPERATURE) = ',I2,'HRS') 00269 26 FORMAT (///,10X,'INTRODUCED PHASE LAG (PRESSURE) = ',I2,'HRS') 00270 75 FORMAT (///,10X,'INTRODUCED PHASE LAG (PERCIPATATION) = ',I2,'RS')00271 80 FORMAT (///, 2X,'OTHER CONSTITUENTS',//) 00272 82 FORMAT (5F15.11) 83 FORMAT (2X,'TEMPERATURE (IN MSEC PER DEGREE CELSIUS):' ,3X,F8.2,' 00273 &+/-',F8.2/) 00274 84 FORMAT (2X,'PRESSURE (IN MSEC PER MILLIBAR):',3X,F8.2,'+/-', 00275 &F8.2/) 00276 85 FORMAT(//2X,'DATUM BIAS = ',F14.3,' +/-',F14.3,' MSEC') 00277 86 FORMAT(//2X,'LINEAR TREND = ',F14.3,' +/-',F14.3,' MSEC./DAY.') 00278 87 FORMAT (2X, 'PRECIPITATION (IN MSEC PER MILLIMETRE):', 3X, F7.3, 00279 &'+/-', F7.3/) 00280 88 FORMAT (2X, 'CORRELATION COEFFICIENT BETWEEN',2X,I2,2X,'AND',2X, 00281 &I2,2X,'EQUALS TO ', F6.4/) 00282 89 FORMAT (///,2X,'HIGH CORRELATION BETWEEN CONSTITUENTS'//) 00283 112 FORMAT (2X, F14.6, 17X, F8.2,' +/-', F8.2, 5X, F8.3,'+/-', F8.3, 00284 &6X/) 00285 113 FORMAT (///,2X,'PERIODIC CONSTITUENTS',10X,'AMPLITUDE IN MSEC', 00286 &9X, 'PHASE IN DEGREES' //) 00287 116 FORMAT (2X,'PHASE TAKEN WITH RESPECT TO 23-RD HOUR OF DAY',I4//) 00288 301 FORMAT (40I2) 00289 101 FORMAT(' ','VALMIN= ',F8.2,'VALMAX= ',F8.2) END 00290 SUBROUTINE LSAPRO (LAST,ISTOP,IY,XBACK) C 00341 C PROGRAM 'LSAPROX' LEAST SQUARES APPROXIMATION 00342 C 00343 C C PROGRAM LSAPRO IS AN APPROXIMATION PROGRAM IN WHICH THE C MAIN PORTION OF THE ALOGRITHEM DETERMINING THE COEFF. HAS C PREVIOUSLY BEEN CALCULATED AND THE ANSWERS INPUTTED AS C 'XBACK', I.E. BACKWARD PREDICTION C FOR FURTHER REFERENCE AS HOW TO CALCULATE THIS .....SEE C H. KARABELA 'LEAST SQ. APPROXIMATION', DEPT OF SURVEYING C ENGINEERING, U.N.B. (1976) C C IMPLICIT REAL*8(A-H, O-Z) 00313 DIMENSION FPHI(25),FX(25),XBACK(25) COMMON E(3,6500),TM C C KOUNT = -1 C IBEGIN = IY + (ISTOP - LAST) / 24.0D0 IBGHR = ((IY + (ISTOP - LAST) / 24.0D0)- IBEGIN) * 24.0D0 C PRINT 55,LAST,IBEGIN,IBGHR C DO 11 I=1,25 11 FX(I) = E(1,ISTOP + I) C 99 CONTINUE C PREDIC = 0.0D0 C C DO 22 I=1,25 22 PREDIC = XBACK(I) * FX(I) + PREDIC C DO 33 I=1,24 33 FPHI ( I + 1) = FX (I) C DO 44 I=1,24 44 FX(I + 1) = FPHI(I + 1) FX(1) = PREDIC C C E (1,ISTOP + KOUNT + 1) = PREDIC C KOUNT = KOUNT - 1 C IF (KOUNT.LT.-LAST) RETURN GO TO 99 55 FORMAT (/,5X,'PREDICTED TEMP.FOR (',I3,') HRS',/,5X, @'STARTING AT DAY = ',I4,4X,'HOUR = ',I3) END SUBROUTINE AMPLAG (A,B,SA,SB,SAB,R,P,SR,SP) 00291 C 00292 C THIS SUBROUTINE COMPUTES THE AMPLITUDES AND PHASE LAGS (AND THEIR 00293 C STANDARD DEVIATIONS) OF PERIODIC CONSTITUENTS. 00294 C 00295 DOUBLE PRECISION A,B,SA,SB,SAB,R,P,SR,SP 00296 R = DSQRT (A*A +B*B) 00297 P = DATAN2 (B,A) 00298 Q = 1.0/(R*R) 00299 SR = (A*(A*SA + 2.0*B*SAB) + B*B*SB) * Q 00300 SP = (A*(A*SB - 2.0*B*SAB) + B*B*SA) * Q * Q 00301 SR = DSQRT (SR) 00302 SP = DSQRT (SP) 00303 P = P*57.29578 00304 SP = SP*57.29578 00305 RETURN 00306 END 00307 SUBROUTINE FIPRIN (N,K,IY) 00308 C 00309 C THIS SUBROUTINE PRINTS OUT OBSERVED TEMPERATURE, PRESSURE AND 00310 C PRECIPITATION (IF REQUIRED) BY DAYS 00311 C 00312 IMPLICIT REAL * 8(A-H,O-Z) DIMENSION FI(24),IH(24) 00314 COMMON E(3,6500),TM 00315*17 IF(K.EQ. 1) PRINT 5 00316 IF(K.EQ. 2) PRINT 6 00317 IF(K.EQ. 3) PRINT 7 00318 DO 9 I = 1,24 00319 9 IH(I) = I - 1 00320 PRINT 8,(IH(I),I = 1,24) 00321 NN = 0.5*N + 0.99 00322 DO 1 I = 1,NN 00323 II = IY + I - 1 00324 DO 2 J = 1,24 00325 L = I*24 - 24 + J 00326 2 FI(J) = E(K,L) 00327 IF(K.EQ.1) PRINT 3, II,(FI(J), J =1,24) 00328 IF(K.EQ.2) PRINT 4, II,(FI(J), J =1,24) 00329 IF(K.EQ.3) PRINT 10,II,(FI(J), J=1,24) 00330 1 CONTINUE 00331 RETURN 00332 3 FORMAT (3X,I3,2X,24F5.1) 00333 4 FORMAT (3X,I3,2X,24F5.0) 00334 5 FORMAT (/// 10X, 'OBSERVED TEMPERATURE IN DEGREES CELSIUS' //) 00335 6 FORMAT (/// 10X, 'OBSERVED PRESSURE IN MILLIBARS' //) 00336 7 FORMAT (/// 10X, 'INTERPOLATED PRECIPITATION IN MILLIMETRES' //) 00337 8 FORMAT (3X,'DAY',24I5//) 00338 10 FORMAT (3X,I3,2X,24F5.3) 00339 END 00340 SUBROUTINE SPECEQ(F,T,NF,PL,PS,NW,IB,DAT,NDAT,LT,PER, 00344 & NPER,NBASE,C,FNORM,A,IY) 00345 C 00346 C LEAST SQUARES SPECTRAL ANALYSIS OF PIECEWISE EQUALLY 00347 C SPACED TIME SERIES 00348 C 00349 C INPUT ARGUMENTS 00350 C F(NF) = VECTOR OF TIME SERIES VALUES 00351 C T(NF) = VECTOR OF TIME SERIES TIMES 00352 C PL = LONGEST PERIOD IN SPECTRAL BAND TO BE COMPUTED 00353 C PS = SHORTEST PERIOD IN SPECTRAL BAND TO BE COMPUTED 00354 C NW = NUMBER OF SPECTRAL VALUES IN BAND TO BE COMPUTED 00355 C IB = BAND TO BE COMPUTED (COMPUTE RESIDUAL TIME SERIES 00356 C ONLY WHEN IB = 1) 00357 C DAT(NDAT) = VECTOR OF TIMES BEGINNING NEW DATUMS 00358 C LT = LINEAR TREND SWITCH(1=FORCE,0=OMIT) 00359 C PER(NPER) = VECTOR OF PERIODS TO BE FORCED 00360 C NBASE = NUMBER OF USER-SPECIFIED CONSTITUENTS 00361 C 00362 IMPLICIT REAL * 8(A-H,O-Z) 00363 REAL F2(6500),SUM,VALMIN,VALMAX,SCALE DIMENSION MINDAT(25) DIMENSION F(NF),T(NF) 00364 DIMENSION IPLOT(100) 00365 DIMENSION XSCLE (1600) C DIMENSION FOLLOWING .GE. NDAT1 = NDAT + 1 00366 DIMENSION DAT(55) 00367 C DIMENSION FOLLOWING .GE. NPER = NUMBER OF FORCED PERIODS 00368 DIMENSION PER(6),FF(6),SM(6),CM(6),SP(6),CP(6), 00369 * CMM(6),SMM(6),SQM(6),SQP(6) 00370 C DIMENSION FOLLOWING .GE. NINT=NUMBER OF SUBINTERVALS IN F 00371 DIMENSION XN(55),XL(55),INT(55) 00372 C DIMENSION FOLLOWING .GE. (NINT,NPER) 00373 DIMENSION CLM(55,6),SLM(55,6),CNM(55,6),SNM(55,6) 00374 DIMENSION A(65,65),B(65),C(65),U(65),V(65) 00376 C DIMENSION FOLLOWING .GE. M=NDAT+LT+2*NPER+NBASE 00375 COMMON E(3,6500),TM 00377*17 COMMON / NUMBR/ N,SCALE,IRESID DATA IPR/6/, IBLANK/1H /,ISTAR/1H*/ 00378 C 00379 SQRT(AA) = DSQRT(AA) 00380 SIN(AA) = DSIN(AA) 00381 COS(AA) = DCOS(AA) 00382 ABS(AA) = DABS(AA) 00383 SIGN(A,B) = DSIGN(A,B) 00384 C MACHINE DEPENDENT CONSTANTS 00385 C DIVIDE=MINIMUM DIVISOR. SING=MINIMUM MATRIX DETERMINANT 00386 PI = 3.141592653589793D0 00387 DIVIDE = 1.D-35 00388 SING = 1.0D-55 00389*17 C 00390 C SKIP PRELIMINARY COMPUTATIONS IF IB .NE. 1 00391 IF(IB .NE. 1) GO TO 200 00392 C 00393 C INITIALIZATION 00394 C 00395 C COMPUTE TIME SERIES LENGTH 'CT' AND SPACING 'STEP' 00396 CT = T(NF) - T(1) 00397 XNF = NF 00398 STEP = T(2) - T(1) 00399 C CHECK DAT(1) .EQ. T(1) AND SET DAT(NDAT+1) .GT. T(NF) 00400 IF(NDAT .GE. 1 .AND. DAT(1) .NE. T(1)) GO TO 310 00401 NDAT1 = NDAT + 1 00402 DAT(NDAT1) = T(NF) + 1. 00403 C INITIALIZE PLOT ARRAY 00404 DO 1 I = 1,100 00405 1 IPLOT(I) = IBLANK 00406 WRITE (IPR,414) 00407 WRITE (IPR,401) 00408 C FIND MAXIMUM AND MINIMUM VALUES IN VECTOR 'F' 00409 NMIN=0 FMIN = 1. / DIVIDE 00410 FMAX = - FMIN 00411 TMIN=1./DIVIDE TMAX=-TMIN OLDMAX=0. OLDMIN=0. DO 2 I=1,NF IF(I.EQ.1) GO TO 550 GAP=(T(I)-T(I-1))/STEP-1. IF(GAP.LT.0.5) GO TO 550 IF(FMAX-FMIN.LT.OLDMAX-OLDMIN) GO TO 551 OLDMAX=FMAX OLDMIN=FMIN 551 NMIN=NMIN+1 MINDAT(NMIN)=FMIN FMIN=1./DIVIDE FMAX=-FMIN 550 IF(F(I).GT.FMAX) FMAX=F(I) IF(F(I).LT.FMIN) FMIN=F(I) IF(F(I).GT.TMAX) TMAX=F(I) IF(F(I).LT.TMIN) TMIN=F(I) 2 CONTINUE NMIN=NMIN+1 MINDAT(NMIN)=FMIN IF(FMAX-FMIN.GT.OLDMAX-OLDMIN) GO TO 552 FMAX=OLDMAX FMIN=OLDMIN 552 CONTINUE FMID=(TMAX+TMIN)/2. NMIN=0 C SCAN THROUGH TIME SERIES 00417 X = 99.0/ (FMAX-FMIN) 00418 IDAT = 1 00420 NINT = 0 00421 TA = T(1) 00422 DO 6 I = 1,NF 00423 NEWSUB = 0 00424 IF(I .EQ. NF) NEWSUB = 1 00425 C CHECK T INCREASING MONOTONICALLY 00426 IF(I .EQ. 1) GO TO 4 00427 IF(T(I) .LE. T(I-1)) GO TO 320 00428 C FIND GAPS IN TIME SERIES 00429 GAP = (T(I) - T(I-1)) / STEP - 1. 00430 IF(GAP .LT. 0.5) GO TO 4 00431 WRITE (IPR,402)GAP 00432 NEWSUB = 1 00433 C CHECK FOR NEW DATUM 00434 4 IF(DAT(IDAT) .GT. T(I)) GO TO 5 00435 NMIN=NMIN+1 Y=1.0-MINDAT(NMIN)*X WRITE (IPR,403) T(I), IDAT 00436 IDAT = IDAT + 1 00437 NEWSUB = 2 00438 C PLOT TIME SERIES 00439 5 KF = F(I)*X+Y 00440 IF(KF.LT. 1) KF = 1 00441 IF(KF.GT. 100) KF = 100 00442 IPLOT(KF) = ISTAR 00443 LL = (T(I) - 1.0)/24.0 + 0.001 00444 ID = IY + LL 00445 LL = T(I) - 24.0 * DFLOAT(LL) - 0.999 00446 IF(LL.EQ.0)WRITE(IPR,404)ID,LL,F(I),IPLOT 00447 IF(LL.NE.0)WRITE(IPR,410)LL,F(I),IPLOT 00448 IPLOT(KF) = IBLANK 00449 C NEW SUBINTERVAL IF GAP OR NEW DATUM 00450 C COMPUTE NUMBER OF SUBINTERVALS = NINT 00451 C COMPUTE VECTORS XN(NINT),XL(NINT),INT(NINT) 00452 IF(I .EQ. 1 .OR. NEWSUB .EQ. 0) GO TO 6 00453 TB = T(I-1) 00454 IF(I .EQ. NF) TB = T(NF) 00455 NINT = NINT + 1 00456 XN(NINT) = (TB - TA + STEP) / STEP 00457 XL(NINT) = (TB + TA) / STEP 00458 INT(NINT) = IDAT - NEWSUB 00459 TA = T(I) 00460 6 CONTINUE 00461 C COMPUTE QUADRATIC NORM OF ORIGINAL TIME SERIES 00462 FNORM = 0.0 00463 FBAR = 0.0 00464 DO 500 I = 1,NF 00465 FBAR = FBAR + F(I) 00466 500 FNORM = FNORM + F(I) * F(I) 00467 FNORM = FNORM - (FBAR * FBAR) / NF 00468 WRITE (IPR,501) FNORM 00469 SPP = SQRT(FNORM/NF) 00470 WRITE (IPR,600) SPP 00471 C COMPUTE FORCED FREQUENCIES 'FF' 00472 C SHORTEST PERMITTED FORCED PERIOD IS 2.*STEP 00473 IF(NPER .EQ. 0) GO TO 8 00474 DO 7 K = 1,NPER 00475 IF(PER(K) .LT. 2.*STEP) GO TO 350 00476 FF(K) = 2. * PI / PER(K) 00477 C COMPUTE TRIG(PK), TRIG(N*PK), TRIG(L*PK) 00478 PK = PI * STEP / PER(K) 00479 SM(K) = SIN(PK) 00480 CM(K) = COS(PK) 00481 DO 7 J = 1,NINT 00482 PKN = XN(J) * PK 00483 PKL = XL(J) * PK 00484 SNM(J,K) = SIN(PKN) 00485 CNM(J,K) = COS(PKN) 00486 SLM(J,K) = SIN(PKL) 00487 7 CLM(J,K) = COS(PKL) 00488 C NUMBER OF KNOWN CONSTITUENTS 00489 8 M = NDAT + LT + 2 * NPER + NBASE 00490 WRITE(IPR,405) NDAT,LT,NPER,NPER,NBASE,M 00491 IF(M .EQ. 0) GO TO 126 00492 C 00493 C SOLVE NORMAL EQUATIONS FOR KNOWN CONSTITUENTS 00494 C 00495 C CLEAR NORMAL EQUATION ARRAYS 00496 DO 100 I = 1,M 00497 B(I) = 0. 00498 DO 100 J = 1,M 00499 100 A(I,J) = 0. 00500 C CONSTRUCT MATRIX 'A' AND VECTOR 'B' 00501 DO 101 I = 1,NF 00502 DO 101 J = 1,M 00503 AA= BASE(J,T(I),DAT,NDAT,LT,FF,NPER) 00504 B(J) = B(J) + AA* F(I) 00505 DO 101 K = J,M 00506 101 A(K,J) = A(K,J) + AA * 00507 * BASE(K,T(I),DAT,NDAT,LT,FF,NPER) 00508 CALL MOUTD(A,65,M,M) C 00510*17 C C C 00542*17 C 00544 C CHOLESKI DECOMPOSITION OF MATRIX 'A' TO TRIANGULAR MATRIX 00545 DET=1.0D0 A(1,1) = SQRT(A(1,1)) 00547 IF(M .EQ. 1) GO TO 109 00548 DO 104 I = 2,M 00549 104 A(I,1)=A(I,1)/A(1,1) 00550 DO 108 J = 2,M 00551 SUM=0. 00552 J1=J-1 00553 DO 105 K = 1,J1 00554 105 SUM=SUM+A(J,K)*A(J,K) 00555 A(J,J)= SQRT(A(J,J)-SUM) 00557 IF(J.EQ.M) GO TO 108 00558 J2=J+1 00559 DO 107 I = J2,M 00560 SUM=0. 00561 DO 106 K = 1,J1 00562 106 SUM=SUM+A(I,K)*A(J,K) 00563 107 A(I,J)=(A(I,J)-SUM)/A(J,J) 00564 108 CONTINUE 00565 NDIM=65 CALL DETM(A,NDIM,M,RPART,IPART) WRITE(IPR,709)RPART,IPART 109 IF(ABS(DET) .LT. SING) GO TO 330 00567 C INVERT LOWER TRIANGULAR MATRIX 00568 DO 110 I = 1,M 00569 110 A(I,I) = 1. / A(I,I) 00570 IF(M .EQ. 1) GO TO 113 00571 N1 = M - 1 00572 DO 112 J = 1,N1 00573 J2=J+1 00574 DO 112 I = J2,M 00575 SUM=0. 00576 I1=I-1 00577 DO 111 K = J,I1 00578 111 SUM=SUM+A(I,K)*A(K,J) 00579 112 A(I,J)=-A(I,I)*SUM 00580 C CONSTRUCT INVERSE OF MATRIX 'A' 00581 113 DO 118 J = 1,M 00582 IF(J.EQ.1) GO TO 115 00583 J1=J-1 00584 DO 114 I = 1,J1 00585 114 A(I,J)=A(J,I) 00586 115 DO 117 I = J,M 00587 SUM=0. 00588 DO 116 K = I,M 00589 116 SUM=SUM+A(K,I)*A(K,J) 00590 117 A(I,J)=SUM 00591 118 CONTINUE 00592 CALL MOUTD(A,65,M,M) C 00594*17 C C C C 00601 C COMPUTE COEFFICIENTS 'C' 00602 DO 119 I = 1,M 00603 C(I) = 0. 00604 DO 119 J = 1,M 00605 119 C(I) = C(I) + A(I,J) * B(J) 00606 C WRITE PRELIMINARY AMPLITUDES OF KNOWN CONSTITUENTS 00607 IF(NDAT .GE. 1) WRITE(IPR,406) (I,C(I),I=1,NDAT) 00608 IF(LT .EQ. 1) WRITE(IPR,407) NDAT1,C(NDAT1) 00609 IF(NPER .EQ. 0) GO TO 121 00610 DO 120 I = 1,NPER 00611 K = NDAT + LT + 2 * I - 1 00612 K1 = K + 1 00613 120 WRITE(IPR,408) K,K1,PER(I),C(K),C(K1) 00614 121 IF(NBASE .EQ. 0) GO TO 123 00615 DO 122 I = 1,NBASE 00616 K = NDAT + LT + 2 * NPER + I 00617 122 WRITE(IPR,409) K,C(K) 00618 123 CONTINUE 00619 C COMPUTE RESIDUAL TIME SERIES 'F' 00620 DO 125 I = 1,NF 00621 PM = 0. 00622 DO 124 J = 1,M 00623 124 PM = PM + C(J) * BASE(J,T(I),DAT,NDAT,LT,FF,NPER) 00624 125 F(I) = F(I) - PM 00625 C PLOT RESIDUAL TIME SERIES 00626 C FIND MAX AND MIN OF RESIDUAL SERIES 00627 FMIN = 1.0/DIVIDE 00628 FMAX = - FMIN 00629 SUM=0. DO 650 I = 1, NF 00630 IF (F(I).GT.FMAX)FMAX = F(I) 00631 IF (F(I).LT. FMIN) FMIN = F(I) 00632 F2(I)=F(I) SUM=SUM+F(I) 650 CONTINUE 00633 SUM=SUM/DFLOAT(NF) VALMIN=FMIN VALMAX=FMAX NN=0.5*N+0.99 IF(IRESID.EQ.1) CALL VALPLT(IY,NN,VALMIN,VALMAX,SCALE,4) X = 99.0/(FMAX-FMIN) 00634 Y = 1.0 - FMIN*X 00635 WRITE (IPR, 415) 00636 IDAT = 1 00637 WRITE (IPR,401) 00638 DO 651 I = 1,NF 00639 IF(I.EQ.1) GO TO 652 00640 GAP = (T(I)-T(I-1))/STEP-1.0 00641 IF(GAP.LT.0.5)GO TO 652 00642 WRITE (IPR,402) GAP 00643 652 IF(DAT(IDAT).GT.T(I))GO TO 653 00644 WRITE(IPR,403)T(I),IDAT 00645 IDAT = IDAT + 1 00646 653 KF = F(I)*X+Y 00647 IF(KF.LT.1) KF = 1 00648 IF(KF. GT. 100) KF = 100 00649 IPLOT (KF) = ISTAR 00650 LL = (T(I) - 1.0)/24.0 + 0.001 00651 ID = IY + LL 00652 LL = T(I) - 24.0 * DFLOAT(LL) - 0.999 00653 IF(LL.EQ.0.)WRITE(IPR,404)ID,LL,F(I),IPLOT 00654 IF(LL.NE.0)WRITE(IPR,410)LL,F(I),IPLOT 00655 IPLOT(KF) = IBLANK 00656 651 CONTINUE 00657 C COMPUTE QUADRATIC NORM OF 'F' 00658 126 FNORM = 0. 00659 DO 127 I = 1,NF 00660 127 FNORM = FNORM + F(I) * F(I) 00661 WRITE(IPR,411) FNORM 00662 SPP = SQRT(FNORM/NF) 00663 WRITE (IPR,600) SPP 00664 C ASSUME RESIDUAL F CONSISTS OF ROUNDOFF NOISE WHEN 00665 C RMS OF RESIDUAL F .LT. 0.1% OF MID-RANGE OF ORIGINAL F 00666 IF(FNORM .LT. XNF * (0.001 * FMID)**2) GO TO 340 00667 C 00668 C ENTRY FOR SUBSEQUENT SPECTRAL BANDS TO BE COMPUTED 00669 C 00670 200 CONTINUE 00671 C DEFAULT VALUE OF PL IS CT 00672 IF(PL .EQ. 0.) PL = CT 00673 C SHORTEST PERMITTED PERIOD IS 2 STEPS 00674 IF(PS .LT. 2.*STEP) PS = 2. * STEP 00675 OM1 = 2. * PI / PL 00676 C DEFAULT VALUE OF NW IF INPUTTED AS ZERO 00677 IF(NW .EQ. 0) NW = (2. * PI / PS - OM1) * CT / 2. 00678 DOM = (2. * PI / PS - OM1) / (NW - 1.) 00679 WRITE(IPR,412) IB,NW,PL,PS 00680 DO 700 I = 1,100 00681 700 IPLOT(I) = IBLANK 00682 C 00683 C LOOP OVER SPECTRAL VALUES TO BE COMPUTED 00684 C 00685 DO 217 I = 1,NW 00686 OMEGA = OM1 + DOM * (I - 1.) 00687 FCOS = 0. 00688 FSIN = 0. 00689 S1 = 0. 00690 S2 = 0. 00691 CU = 0. 00692 CV = 0. 00693 CW = 0. 00694 C CLEAR VECTORS 'U','V' 00695 IF(M .EQ. 0) GO TO 202 00696 DO 201 J = 1,M 00697 U(J) = 0. 00698 201 V(J) = 0. 00699 C COMPUTE SCALAR PRODUCTS 'FCOS' AND 'FSIN' 00700 202 DO 204 J = 1,NF 00701 WT = OMEGA * T(J) 00702 COSWT = COS(WT) 00703 SINWT = SIN(WT) 00704 FCOS = FCOS + F(J) * COSWT 00705 FSIN = FSIN + F(J) * SINWT 00706 C COMPUTE USER-DEFINED SCALAR PRODUCTS IN 'U','V' 00707 IF(NBASE .EQ. 0) GO TO 204 00708 DO 203 L = 1,NBASE 00709 K = NDAT + LT + 2 * NPER + L 00710 AA = BASE(K,T(J),DAT,NDAT,LT,FF,NPER) 00711 U(K) = U(K) + AA * COSWT 00712 203 V(K) = V(K) + AA * SINWT 00713 204 CONTINUE 00714 C COMPUTE TRIG(Q) 00715 Q = OMEGA * STEP / 2. 00716 SQ = SIN(Q) 00717 CQ = COS(Q) 00718 SQIN = 1. / SQ 00719 CQIN = 1. / CQ 00720 COTQ = CQ * SQIN 00721 C INITIALIZE FORCED FREQUENCY PARAMETERS 00722 IF(NPER .EQ. 0) GO TO 206 00723 DO 205 J = 1,NPER 00724 X = SM(J) * CQ 00725 Y = CM(J) * SQ 00726 SQM(J) = X - Y 00727 C ENSURE SQM .GE. DIVIDE 00728 IF(ABS(SQM(J)) .LT. DIVIDE) SQM(J) = 00729 * SIGN(DIVIDE,SQM(J)) 00730 SQP(J) = X + Y 00731 CP(J) = 0. 00732 SP(J) = 0. 00733 CMM(J) = 0. 00734 205 SMM(J) = 0. 00735 C COMPUTE TRIG(N*Q), TRIG(L*Q) AND THEIR PRODUCTS 00736 206 DO 210 J = 1,NINT 00737 QN = XN(J) * Q 00738 QL = XL(J) * Q 00739 SNQ = SIN(QN) 00740 CNQ = COS(QN) 00741 SLQ = SIN(QL) 00742 CLQ = COS(QL) 00743 X1 = SNQ * SLQ 00744 X2 = CNQ * CLQ 00745 X3 = SNQ * CLQ 00746 X4 = CNQ * SLQ 00747 C SUMMATIONS FOR SCALAR PRODUCTS 'S1', 'S2' , 'S3' 00748 S2 = S2 + X3 * X4 00750 S1 = S1 + X2 * X3 - X1 * X4 00749 C COMPUTE DATUM BIAS SCALAR PRODUCTS IN 'U','V' 00751 IF(M .EQ. 0) GO TO 210 00752 IF(NDAT .EQ. 0) GO TO 207 00753 K = INT(J) 00754 U(K) = U(K) + SQIN * X3 00755 V(K) = V(K) + SQIN * X1 00756 C SUMMATIONS FOR LINEAR TREND SCALAR PRODUCTS IN 'U','V' 00757 207 IF(LT .EQ. 0) GO TO 208 00758 U(NDAT1) = U(NDAT1) - COTQ*X1+XN(J)*X4+XL(J)*X3 00759 V(NDAT1) = V(NDAT1) + COTQ*X3-XN(J)*X2+XL(J)*X1 00760 C COMPUTE FORCED FREQUENCY PARAMETERS 00761 208 IF(NPER .EQ. 0) GO TO 210 00762 DO 209 K = 1,NPER 00763 X = SNM(J,K) * CNQ 00764 Y = CNM(J,K) * SNQ 00765 SNQM = X - Y 00766 SNQP = X + Y 00767 X = SLM(J,K) * CLQ 00768 Y = CLM(J,K) * SLQ 00769 SLQM = X - Y 00770 SLQP = X + Y 00771 X = CLM(J,K) * CLQ 00772 Y = SLM(J,K) * SLQ 00773 CLQM = X + Y 00774 CLQP = X - Y 00775 CMM(K) = CMM(K) + SNQM * CLQM 00776 SMM(K) = SMM(K) + SNQM * SLQM 00777 CP(K) = CP(K) + SNQP * CLQP 00778 209 SP(K) = SP(K) + SNQP * SLQP 00779 210 CONTINUE 00780 C COMPUTE SCALAR PRODUCTS 'S1','S2', AND 'S3' 00781 AA = S1 * SQIN * CQIN 00782 S1 = (XNF + AA) / 2. 00783 S3 = (XNF - AA) / 2. 00784 S2 = S2 * SQIN * CQIN 00785 IF(M .EQ. 0) GO TO 215 00786 C COMPUTE LINEAR TREND SCALAR PRODUCT IN 'U','V' 00787 IF(LT .EQ. 0) GO TO 211 00788 AA = SQIN * STEP / 2. 00789 U(NDAT1) = U(NDAT1) * AA 00790 V(NDAT1) = V(NDAT1) * AA 00791 C COMPUTE FORCED FREQUENCY SCALAR PRODUCTS IN 'U','V' 00792 211 IF(NPER .EQ. 0) GO TO 213 00793 DO 212 J = 1,NPER 00794 SQMIN = 0.5 / SQM(J) 00795 SQPIN = 0.5 / SQP(J) 00796 K = NDAT + LT + 2 * J - 1 00797 U(K ) = SQMIN * CMM(J) + SQPIN * CP(J) 00798 V(K ) = -SQMIN * SMM(J) + SQPIN * SP(J) 00799 U(K+1) = SQMIN * SMM(J) + SQPIN * SP(J) 00800 212 V(K+1) = SQMIN * CMM(J) - SQPIN * CP(J) 00801 C COMPUTE BILINEAR FORMS 'CU','CV','CW' 00802 213 DO 214 J = 1,M 00803 DO 214 K = 1,M 00804 CU = CU + A(J,K) * U(J) * U(K) 00805 CV = CV + A(J,K) * V(J) * V(K) 00806 214 CW = CW + A(J,K) * U(J) * V(K) 00807 C COMPUTE PERIOD AND SPECTRAL VALUE (PERCENTAGE VARIANCE) 00808 215 PERIOD = 2. * PI / OMEGA 00809 SIG = 0. 00810 D = (S1-CU) * (S3-CV) - (S2-CW) * (S2-CW) 00811 IF(ABS(D) .LT. DIVIDE) GO TO 216 00812 SIG = ((S3-CV)*FCOS*FCOS -2.*(S2-CW)*FCOS*FSIN + 00813 * (S1-CU)*FSIN*FSIN) / (D*FNORM) 00814 SIG = SIG * 100. 00815 C WRITE PERIOD AND SPECTRAL VALUE 00816 216 KSIG = SIG 00817 IF(KSIG .LT. 1) KSIG = 1 00818 IF(KSIG .GT. 100) KSIG = 100 00819 IPLOT(KSIG) = ISTAR 00820 WRITE(IPR,413) PERIOD,SIG,IPLOT 00821 IPLOT(KSIG) = IBLANK 00822 217 CONTINUE 00823 Q = 1.0/(NF - M) 00824 WRITE (IPR,418)Q 00825 418 FORMAT (///,2X,'LEVEL OF RANDOM NOISE SPECTRUM =',F6.4) 00826 RETURN 00827 C 00828 C ERROR STOPS 00829 310 WRITE(IPR,311) 00830 311 FORMAT(9X,26HINPUT ERROR-DAT(1).NE.T(1)) 00831 STOP 00832 320 WRITE(IPR,321) 00833 321 FORMAT(9X,34HINPUT ERROR-STEPNUMBER SEQUENCE AT,2F9.0) 00834 STOP 00835 330 WRITE(IPR,331) 00836 331 FORMAT(9X,40HERROR-NORMAL EQUATION MATRIX IS SINGULAR) 00837 STOP 00838 340 WRITE(IPR,341) 00839 341 FORMAT(9X,40HRESIDUAL TIME SERIES=ROUNDOFF NOISE.STOP) 00840 STOP 00841 350 WRITE(IPR,351) PER(I) 00842 351 FORMAT(9X,35HINPUT ERROR-FORCED PERIOD TOO SHORT,F9.5) 00843 STOP 00844 C FORMAT STATEMENTS 00845 401 FORMAT(2X, 'DAY HOUR TILT(MSEC)' //) 00846 402 FORMAT (/65X,6HGAP OF,F6.0,6H STEPS) 00847 403 FORMAT(2X,E11.4,2X,33(1H-) , 18HBEGINNING OF DATUM, 00848 &I5,39(1H-)/) 00849 404 FORMAT (2X,I4,I2,F10.2,100A1) 00850 405 FORMAT(33H1 SOLUTION FOR KNOWN CONSTITUENTS/19X, 00851 * 5HDATUM,8X,6HLINEAR,8X,6HFORCED,8X,6HCOSINE,10X, 00852 * 4HSINE,10X,4HUSER,5X,5HTOTAL/20X,4HBIAS,9X,5HTREND, 00853 * 8X,6HPERIOD,35X,7HDEFINED//2X,7HNUMBER ,2I14,14X, 00854 * 3I14,I9//2X,9HAMPLITUDE//) 00855 406 FORMAT(7X,I3,E14.6) 00856 407 FORMAT(7X,I3,14X,E14.6) 00857 408 FORMAT(3X,I3,1H-,I3,28X,F14.6,2E14.6) 00858 409 FORMAT(7X,I3,70X,E14.6) 00859 410 FORMAT(6X,I2,F10.2,100A1) 00860 411 FORMAT (///9X,35HRESIDUAL TIME SERIES QUADRATIC NORM,E9.4,/) 00861 412 FORMAT(1H1,10X,13HSPECTRAL BAND,I5/10X,I5, 00862 * 24H SPECTRAL VALUES BETWEEN,2F14.6//9X,6HPERIOD,6X, 00863 * 19HPERCENTAGE VARIANCE) 00864 413 FORMAT(5X,F10.5,F12.3,100A1) 00865 414 FORMAT(///,2X,'PLOT OF THE ANALYSED TIME SERIES',////) 00866 415 FORMAT(///,2X, 'PLOT OF RESIDUAL TIME SERIES',////) 00867 501 FORMAT (///9X,26HTIME SERIES QUADRATIC NORM, E9.4/) 00868 600 FORMAT(9X, 'STANDARD DEVIATION = ', E9.4,//) 00869 709 FORMAT(//5X,'DETERMINANT= ',F7.5,'D',I4) C 00871 END 00872 FUNCTION BASE(I,T,DAT,NDAT,LT,FF,NPER) 00873 C 00874 C COMPUTE BASE FUNCTIONS FOR NORMAL EQUATIONS 00875 C 00876 IMPLICIT REAL * 8(A-H,O-Z) 00877 COMMON E(3,6500),TM 00878*17 DIMENSION DAT(55),FF(6) 00879 SIN(X) = DSIN(X) 00880 COS (X) = DCOS(X) 00881 IF(I - NDAT) 1,1,2 00882 C BASE FUNCTION = DATUM BIAS 00883 1 BASE = 0. 00884 IF(T.GE.DAT(I).AND.T.LT.DAT(I + 1)) BASE = 1.0 00885*17 RETURN 00886 2 IF(I - NDAT - LT) 3,3,4 00887 C BASE FUNCTION = LINEAR TREND 00888 3 BASE = (T - TM) 00889*17 RETURN 00890 4 IF(I - NDAT - LT - 2*NPER) 5,5,8 00891 C BASE FUNCTION = TRIGONOMETRIC FUNCTION OF FORCED FREQUENCY 00892 5 L = I - NDAT - LT + 1 00893 IF(MOD(L,2)) 6,6,7 00894 6 L = L / 2 00895 BASE = COS(FF(L) * T) 00896 RETURN 00897 7 L = L / 2 00898 BASE = SIN(FF(L) * T) 00899 RETURN 00900 C BASE FUNCTIONS = TEMPERATURE, PRESSURE, RIVER DISCHARGES 00901 8 IT = T + 0.1 00902 J = I - NDAT - LT - 2*NPER 00903 BASE = E(J,IT) 00904 RETURN 00905 END 00906 SUBROUTINE DETM(A,NDIM,M,RPART,IPART) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(NDIM,NDIM),VD(40) DO 1 I=1,M 1 VD(I)=A(I,I)**2 DET=0.D0 DO 2 I=1,M DET=DET+DLOG10(VD(I)) 2 CONTINUE IPART=DET Y=DET-DFLOAT(IPART) RPART=10**Y RETURN END SUBROUTINE VALPLT(IY,NN,VALMIN,VALMAX,SCALE,NFLAG) IMPLICIT REAL (A-H,O-Z) REAL*8 E(3,6500),TM DIMENSION FI(24),SMALL(4),BIG(4) COMMON E,TM COMMON /MEANS/SUM COMMON /TILT/F2(6500) DATA SMALL/-40.,950.,0.,0./,BIG/40.,1050.,6000.,200./ ISCALE=SCALE START=IY ENDAY=START+NN WIDTH=NN/SCALE+1. HEIGHT=16. CALL AREA(WIDTH,HEIGHT) CALL SETPLT(START-2.,SMALL(NFLAG),ENDAY,BIG(NFLAG)) CALL SHIFT(START+2.,SMALL(NFLAG)+5.) CALL AXS2(START,SUM,ENDAY,SUM,'PEND NO 2 AT FREDERICTON DAYS 198 &/74 - 85/75\',START,SCALE,-NN/ISCALE) IF(NFLAG.EQ.2) GO TO 2 IF(NFLAG.EQ.3) GO TO 5 IF(NFLAG.EQ.4) GO TO 6 CALL AXS2(START,SMALL(NFLAG),START,BIG(NFLAG),'TEMP IN DEGREES CEL &SIUS\',-40.,5.,16) GO TO 1 2 CALL AXS2(START,SMALL(NFLAG),START,BIG(NFLAG),'PRESSURE IN MILLIBA &RS\',950.,5.,20) GO TO 1 5 CALL AXS2(START,SMALL(NFLAG),START,BIG(NFLAG),'TILT IN MSEC*10\', & 0.,25.,24) GO TO 1 6 CALL AXS2(START,SMALL(NFLAG),START,BIG(NFLAG),'TILT IN MSEC*10\', & 0.,10.,20) 1 DAY =IY NM=0 IPREV=0 DO 3 I=1,NN DO 4 J=1,24 L=I*24-24+J IF(NFLAG.EQ.3)FI(J)=F2(L) 4 IF(NFLAG.NE.3) FI(J)=E(NFLAG,L) DO 100 JJ=1,24 NM=NM+1 IF(NM.NE.1) DAY=DAY+1./24. IFLAG=1 IF(FI(JJ).EQ.0.) GO TO 101 IF(IPREV.EQ.0) IFLAG=0 CALL NOWPLT(IFLAG,DAY,FI(JJ)) IPREV=1 GO TO 100 101 IPREV=0 100 CONTINUE 3 CONTINUE CALL ENDPLT RETURN END