DOUBLE PRECISION A,COVC,SEA,DSQRT,COLT,SELT,COEF,RN,R,P,SR,SP,Q DIMENSION COVC(17,17),A(17,17),DAT(12),PER(5),COEF(17),SEA(17), &F(700),T(700) DIMENSION NAME(5) INTEGER RR,S COMMON E(5,700) RR = 5 C READ THE CONTROL CARD C FOR DESCRIPTION OF PARAMETERS SEE SUBROUTINE SPECEQ READ (RR,4)PL,PS,NW,IB,LT,(PER(I),I=1,5),NDAT,(DAT(I),I=1,5),S 4 FORMAT(2F5.0,2I5,I1,5F6.2,I1,5F3.0,I2) NPER=0 DO 5 I=1,5 IF(PER(I). EQ. 0.0) GO TO 5 NPER = NPER + 1 5 CONTINUE C READ STATION INFORMATION FROM DISC C INO IS STATION NUMBER C NAME IS STATION NAME C IY IS STARTING YEAR C IM IS STARTING MONTH C N IS NUMBER OF MONTHLY VALUES INCLUDING THE MISSING ONES REWIND S READ(S,101) INO 101 FORMAT(I6) DO 102 I=1,5 102 READ(S,103) NAME(I) 103 FORMAT(A4) READ(S,104) IY 104 FORMAT(I4) READ(S,105) IM 105 FORMAT(I2) READ(S,106) N 106 FORMAT(I3) DAT(1)=IM C READ MONTHLY SEA LEVEL VALUES NM=0 DO 2 I = 1,N READ(S,3) Q IF (Q.EQ. 0.0) GO TO 2 NM=NM + 1 F(NM) = Q * 30.48 T(NM)=FLOAT(I) 2 CONTINUE 3 FORMAT (F6.2) NY =N/12 + 0.99 NLY = NY + IY PRINT 7, INO,(NAME(I),I=1,5) PRINT 6, NM, NY, IY, NLY 6 FORMAT(///2X,I6,' MONTHLY MEANS FOR',I4,' YEARS FROM',I5,' TO',I5) 7 FORMAT ('1', 1X, 'STATION', 3X, I6, 2X, 5A4) NL = NDAT + 1 NB = 0 J = 1 C READ MONTHLY TEMPERATURES DO 8 I = 1,N READ (S,3) E(J,I) 8 CONTINUE DO 9 I=1,12 IF(E(J,I))40,9,40 9 CONTINUE GO TO 41 40 NB=NB + 1 J = J + 1 41 CONTINUE C READ MONTHLY PRESSURES DO 10 I = 1,N READ (S,3) E(J,I) 10 CONTINUE DO 11 I = 1,12 IF(E(J,I))42,11,42 11 CONTINUE GO TO 43 42 NB = NB + 1 J = J + 1 43 CONTINUE C READ MONTHLY DISCHARGES C THIS IS THE LAST READ STATEMENT DO 16 K = 1,3 DO 13 I = 1,N READ(S,15) IDISCH E(J,I)=FLOAT(IDISCH) 13 CONTINUE DO 12 I = 1,12 IF(E(J,I))44,12,44 12 CONTINUE GO TO 16 44 NB = NB + 1 J = J + 1 16 CONTINUE IF (J.EQ.2) PRINT 17 17 FORMAT (///5X,'NO RIVER DISCHARGE CONSIDERED'///) 15 FORMAT (I7) M = NDAT + LT + 2*NPER +NB CALL SPECEQ (F,T,NM,PL,PS,NW,IB,DAT,NDAT,LT,PER,NPER,NB,COEF,RN,A) C COMPUTE THE STANDARD DEVIATIONS OF COEFFICIENTS 69 RN=RN/(NM-M) PRINT 19 DO 55 I=1,M DO 55 J=1,M 55 COVC(I,J)=A(I,J)*RN DO 33 L=1,M 33 SEA(L)=DSQRT(COVC(L,L)) SELT=SEA(NL)*1200. COLT=COEF(NL)*1200. DO 21 I = 1,NDAT 21 PRINT 85,COEF(I),SEA(I) IF(LT.EQ.1) PRINT 86,COLT,SELT PRINT 113 113 FORMAT (///,2X,'PERIODIC CONSTITUENTS',15X,'AMPLITUDE IN CM', 9X, &'PHASE IN DEGREES', 8X, 'ASCENDING NODE DATE',//) IY = IY - 1 PRINT 116,IY 116 FORMAT (2X,'PHASE TAKEN WITH RESPECT TO MIDDLE OF DECEMBER',I5//) J = 2 + LT K = J + 2*NPER - 1 DO 110 I = J,K,2 CALL AMPLAG (COEF(I), COEF(I+1), COVC(I,I), COVC(I+1,I+1), COVC &(I+1,I),R,P,SR,SP) L = (I - LT)/2 PY = P*2.3148 E-4*PER(L) DAS = FLOAT (IY) + 0.958333 + PY - 0.0208333 * PER (L) SPY = SP * PER (L) * 2.3148 E-4 PRINT 112, PER(L),R,SR,P,SP,DAS,SPY 110 CONTINUE 112 FORMAT (2X, F6.2, 25X, F8.5,' +/-', F8.5, 5X, F8.3,' +/-', F8.3, &6X, F9.4,' +/- ', F7.4,/) PRINT 80 80 FORMAT (///, 2X,'OTHER CONSTITUENTS',//) K=K+1 COEF(K) = COEF(K) * 1.8 SEA (K) = SEA (K) * 1.8 PRINT 83, COEF(K), SEA(K) 83 FORMAT (2X,'TEMPERATURE (IN CM PER DEGREE CELSIUS):',15X,F8.5, ' &+/-',F8.5/) K=K+1 COEF (K) = COEF (K)/33.863 SEA (K) = SEA (K)/33.863 PRINT 84, COEF(K), SEA(K) 84 FORMAT (2X,'PRESSURE (IN CM PER MILLIBAR):',24X,F8.5,' +/-',F8.5/) K = K + 1 IF(K.GT.M) GO TO 114 DO 115 I = K,M COEF(I) = COEF(I)*35.315 SEA(I) = SEA(I)*35.315 PRINT 87, COEF (I), SEA (I) 115 CONTINUE 114 CONTINUE 87 FORMAT (2X,'RIVER DISCHARGE (IN CM PER CUBIC METRE PER SECOND):', &3X,F8.5,' +/-',F8.5/) PRINT 89 89 FORMAT (///,2X,'HIGH CORRELATION BETWEEN CONSTITUENTS'//) MM = M - 1 DO 120 I = 1,MM K = I + 1 DO 120 J = K,M Q = DSQRT (COVC(I,I) * COVC (J,J)) R = COVC (I,J)/Q Q = DABS (R) IF (Q.LT.0.5)GO TO 120 PRINT 88, I,J,R 120 CONTINUE 88 FORMAT (2X, 'CORRELATION COEFFICIENT BETWEEN',2X,I2,2X,'AND',2X, &I2,2X,'EQUALS TO ', F6.4/) 98 PRINT 19 19 FORMAT('1') 82 FORMAT(2(10D10.2/)) 85 FORMAT(//2X,'DATUM BIAS = ',F8.3,' +/-',F8.3,' CM') 86 FORMAT(//2X,'LINEAR TREND = ',F8.3,' +/-',F8.3,' CM/CENTURY') STOP END SUBROUTINE AMPLAG (A,B,SA,SB,SAB,R,P,SR,SP) DOUBLE PRECISION A,B,SA,SB,SAB,R,P,SR,SP R = DSQRT (A*A +B*B) P = DATAN2 (B,A) Q = 1.0/(R*R) SR = (A*(A*SA + 2.0*B*SAB) + B*B*SB) * Q SP = (A*(A*SB - 2.0*B*SAB) + B*B*SA) * Q * Q SR = DSQRT (SR) SP = DSQRT (SP) P = P*57.29578 SP = SP*57.29578 RETURN END SUBROUTINE SPECEQ(F,T,NF,PL,PS,NW,IB,DAT,NDAT,LT,PER, 0010 * NPER,NBASE,C,FNORM,A) C 0030 C LEAST SQUARES SPECTRAL ANALYSIS OF PIECEWISE EQUALLY 0040 C SPACED TIME SERIES 0050 C 0060 C INPUT ARGUMENTS 0070 C F(NF) = VECTOR OF TIME SERIES VALUES 0080 C T(NF) = VECTOR OF TIME SERIES TIMES 0090 C PL = LONGEST PERIOD IN SPECTRAL BAND TO BE COMPUTED 0100 C PS = SHORTEST PERIOD IN SPECTRAL BAND TO BE COMPUTED 0110 C NW = NUMBER OF SPECTRAL VALUES IN BAND TO BE COMPUTED 0120 C IB = BAND TO BE COMPUTED (COMPUTE RESIDUAL TIME SERIES 0130 C ONLY WHEN IB = 1) 0140 C DAT(NDAT) = VECTOR OF TIMES BEGINNING NEW DATUMS 0150 C LT = LINEAR TREND SWITCH(1=FORCE,0=OMIT) 0160 C PER(NPER) = VECTOR OF PERIODS TO BE FORCED 0170 C NBASE = NUMBER OF USER-SPECIFIED CONSTITUENTS 0180 C 0190 DIMENSION F(NF),T(NF) DIMENSION FS(700) DIMENSION IPLOT(100) 0210 C DIMENSION FOLLOWING .GE. NDAT1 = NDAT + 1 0220 DIMENSION DAT(12) C DIMENSION FOLLOWING .GE. NPER = NUMBER OF FORCED PERIODS 0240 DIMENSION PER(5),FF(5),SM(5),CM(5),SP(5),CP(5), 0250 * CMM(5),SMM(5),SQM(5),SQP(5) 0260 C DIMENSION FOLLOWING .GE. NINT=NUMBER OF SUBINTERVALS IN F 0270 DIMENSION XN(12),XL(12),INT(12) C DIMENSION FOLLOWING .GE. (NINT,NPER) 0290 DIMENSION CLM(12,5),SLM(12,5),CNM(12,5),SNM(12,5) C DIMENSION FOLLOWING .GE. M=NDAT+LT+2*NPER+NBASE 0310 DIMENSION A(17,17),B(17),C(17),U(17),V(17) DATA IPR/6/, IBLANK/1H /,ISTAR/1H*/ C 0340 DOUBLE PRECISION A,AA,B,C,FNORM,PM,SUM,SQRT,DSQRT, 0350 * FCOS,FSIN,S1,S2,S3,U,V,CU,CV,CW,COSWT,SINWT 0360 SQRT(AA) = DSQRT(AA) 0370 C MACHINE DEPENDENT CONSTANTS 0380 C DIVIDE=MINIMUM DIVISOR. SING=MINIMUM MATRIX DETERMINANT 0390 PI = 3.141592653589793D0 0400 DIVIDE = 1.D-35 0410 C 0430 SING = 1.D-10 0420 C SKIP PRELIMINARY COMPUTATIONS IF IB .NE. 1 0440 IF(IB .NE. 1) GO TO 200 0450 C 0460 C INITIALIZATION 0470 C 0480 C COMPUTE TIME SERIES LENGTH 'CT' AND SPACING 'STEP' 0490 CT = T(NF) - T(1) 0500 XNF = NF 0510 STEP = T(2) - T(1) 0520 C CHECK DAT(1) .EQ. T(1) AND SET DAT(NDAT+1) .GT. T(NF) 0530 IF(NDAT .GE. 1 .AND. DAT(1) .NE. T(1)) GO TO 310 0540 NDAT1 = NDAT + 1 0550 DAT(NDAT1) = T(NF) + 1. 0560 WRITE(IPR,401) STEP 0600 C FIND MAXIMUM AND MINIMUM VALUES IN VECTOR 'F' 0610 FMIN = 1. / DIVIDE 0620 FMAX = - FMIN 0630 DO 2 I = 1,NF 0640 IF(F(I) .GT. FMAX) FMAX = F(I) 0650 IF(F(I) .LT. FMIN) FMIN = F(I) 0660 2 CONTINUE 0670 FMID = (FMAX + FMIN) / 2. 0680 C SCAN THROUGH TIME SERIES 0690 X = 29.0/ (FMAX-FMIN) SCALE = 20.0 IDAT = 1 0720 NINT = 0 0730 TA = T(1) 0740 DO 6 I = 1,NF 0750 NEWSUB = 0 0760 IF(I .EQ. NF) NEWSUB = 1 0770 C CHECK T INCREASING MONOTONICALLY 0780 IF(I .EQ. 1) GO TO 4 0790 IF(T(I) .LE. T(I-1)) GO TO 320 0800 C FIND GAPS IN TIME SERIES 0810 GAP = (T(I) - T(I-1)) / STEP - 1. 0820 IF(GAP .LT. 0.5) GO TO 4 0830 NEWSUB = 1 0850 C CHECK FOR NEW DATUM 0860 4 IF(DAT(IDAT) .GT. T(I)) GO TO 5 0870 WRITE (IPR,403) I,T(I),F(I),IDAT IDAT = IDAT + 1 0890 NEWSUB = 2 0900 5 CONTINUE C COMPUTE NUMBER OF SUBINTERVALS = NINT 0990 C COMPUTE VECTORS XN(NINT),XL(NINT),INT(NINT) 1000 IF(I .EQ. 1 .OR. NEWSUB .EQ. 0) GO TO 6 1010 TB = T(I-1) 1020 IF(I .EQ. NF) TB = T(NF) 1030 NINT = NINT + 1 1040 XN(NINT) = (TB - TA + STEP) / STEP 1050 XL(NINT) = (TB + TA) / STEP 1060 INT(NINT) = IDAT - NEWSUB 1070 TA = T(I) 1080 6 CONTINUE 1090 C PLOT TIME SERIES WRITE (IPR,414) 414 FORMAT(///,2X,'PLOT OF THE ANALYSED TIME SERIES',////) CALL TSPLOT (F,T,NF,X,SCALE) C COMPUTE QUADRATIC NORM OF ORIGINAL TIME SERIES FNORM = 0.0 FBAR = 0.0 DO 500 I = 1,NF FBAR = FBAR + F(I) 500 FNORM = FNORM + F(I) * F(I) FNORM = FNORM - (FBAR * FBAR) / NF WRITE (IPR,501) FNORM 501 FORMAT (///9X,29HTIME SERIES QUADRATIC NORM = ,E9.4/) SPP = SQRT(FNORM/NF) WRITE (IPR,600) SPP 600 FORMAT(9X, 'STANDARD DEVIATION = ', E9.4,//) C COMPUTE FORCED FREQUENCIES 'FF' 1100 C SHORTEST PERMITTED FORCED PERIOD IS 2.*STEP 1110 IF(NPER .EQ. 0) GO TO 8 1120 DO 7 K = 1,NPER 1130 IF(PER(K) .LT. 2.*STEP) GO TO 350 1140 FF(K) = 2. * PI / PER(K) 1150 C COMPUTE TRIG(PK), TRIG(N*PK), TRIG(L*PK) 1160 PK = PI * STEP / PER(K) 1170 SM(K) = SIN(PK) 1180 CM(K) = COS(PK) 1190 DO 7 J = 1,NINT 1200 PKN = XN(J) * PK 1210 PKL = XL(J) * PK 1220 SNM(J,K) = SIN(PKN) 1230 CNM(J,K) = COS(PKN) 1240 SLM(J,K) = SIN(PKL) 1250 7 CLM(J,K) = COS(PKL) 1260 C NUMBER OF KNOWN CONSTITUENTS 1270 8 M = NDAT + LT + 2 * NPER + NBASE 1280 WRITE(IPR,405) NDAT,LT,NPER,NPER,NBASE,M 1290 IF(M .EQ. 0) GO TO 126 1300 C 1310 C SOLVE NORMAL EQUATIONS FOR KNOWN CONSTITUENTS 1320 C 1330 C CLEAR NORMAL EQUATION ARRAYS 1340 DO 100 I = 1,M 1350 B(I) = 0. 1360 DO 100 J = 1,M 1370 100 A(I,J) = 0. 1380 C CONSTRUCT MATRIX 'A' AND VECTOR 'B' 1390 DO 101 I = 1,NF 1400 DO 101 J = 1,M 1410 AA= BASE(J,T(I),DAT,NDAT,LT,FF,NPER) 1420 B(J) = B(J) + AA* F(I) 1430 DO 101 K = J,M 1440 101 A(K,J) = A(K,J) + AA * 1450 * BASE(K,T(I),DAT,NDAT,LT,FF,NPER) 1460 C NK=M C AUTOMATIC SCALING FOR SUBROUTINE SPECEQ C C DETERMINE THE MEAN OF THE DIAGONAL ELEMENTS OF THE 'A' MATRIX C AMN = 0.0D0 XSCLE = 1.0D0 ITRMS = 0 DO 4001 I = 1,NK AMN = AMN + A(I,I) ITRMS = ITRMS + 1 4001 CONTINUE AMN = AMN / ITRMS GAMA = 1.0D0 DO 4005 I = 1,NK GAMA = GAMA * (A(I,I) / AMN) 4005 CONTINUE C C DETERMINE SCALE FACTOR TO DIVIDE ALL THE ELEMENTS C OF THE 'A' MATRIX BY C C SCALE FACTOR = (GAMA ** (1/NK)) TIMES AMN C XSCLE = (GAMA ** (1/NK)) * AMN C C SCALE THE 'A' MATRIX C DO 4010 I = 1,NK DO 4010 J = 1,NK A(I,J) = A(I,J) / XSCLE 4010 CONTINUE C C CHOLESKI DECOMPOSITION OF MATRIX 'A' TO TRIANGULAR MATRIX 1470 DET = A(1,1) 1480 A(1,1) = SQRT(A(1,1)) 1490 IF(M .EQ. 1) GO TO 109 1500 DO 104 I = 2,M 1510 104 A(I,1)=A(I,1)/A(1,1) 1520 DO 108 J = 2,M 1530 SUM=0. 1540 J1=J-1 1550 DO 105 K = 1,J1 1560 105 SUM=SUM+A(J,K)*A(J,K) 1570 DET =DET *(A(J,J)-SUM) 1580 A(J,J)= SQRT(A(J,J)-SUM) 1590 IF(J.EQ.M) GO TO 108 1600 J2=J+1 1610 DO 107 I = J2,M 1620 SUM=0. 1630 DO 106 K = 1,J1 1640 106 SUM=SUM+A(I,K)*A(J,K) 1650 107 A(I,J)=(A(I,J)-SUM)/A(J,J) 1660 108 CONTINUE 1670 109 IF(ABS(DET) .LT. SING) GO TO 330 1680 C INVERT LOWER TRIANGULAR MATRIX 1690 DO 110 I = 1,M 1700 110 A(I,I) = 1. / A(I,I) 1710 IF(M .EQ. 1) GO TO 113 1720 N1 = M - 1 1730 DO 112 J = 1,N1 1740 J2=J+1 1750 DO 112 I = J2,M 1760 SUM=0. 1770 I1=I-1 1780 DO 111 K = J,I1 1790 111 SUM=SUM+A(I,K)*A(K,J) 1800 112 A(I,J)=-A(I,I)*SUM 1810 C CONSTRUCT INVERSE OF MATRIX 'A' 1820 113 DO 118 J = 1,M 1830 IF(J.EQ.1) GO TO 115 1840 J1=J-1 1850 DO 114 I = 1,J1 1860 114 A(I,J)=A(J,I) 1870 115 DO 117 I = J,M 1880 SUM=0. 1890 DO 116 K = I,M 1900 116 SUM=SUM+A(K,I)*A(K,J) 1910 117 A(I,J)=SUM 1920 118 CONTINUE 1930 C C DELETE THE SCALE FACTOR FROM THE INVERTED MATRIX 'A' C DO 4110 I = 1,NK DO 4110 J = 1,NK A(I,J) = A(I,J) / XSCLE 4110 CONTINUE C C COMPUTE COEFFICIENTS 'C' 1940 DO 119 I = 1,M 1950 C(I) = 0. 1960 DO 119 J = 1,M 1970 119 C(I) = C(I) + A(I,J) * B(J) 1980 C WRITE PRELIMINARY AMPLITUDES OF KNOWN CONSTITUENTS 1990 IF(NDAT .GE. 1) WRITE(IPR,406) (I,C(I),I=1,NDAT) 2000 IF(LT .EQ. 1) WRITE(IPR,407) NDAT1,C(NDAT1) 2010 IF(NPER .EQ. 0) GO TO 121 2020 DO 120 I = 1,NPER 2030 K = NDAT + LT + 2 * I - 1 2040 K1 = K + 1 2050 120 WRITE(IPR,408) K,K1,PER(I),C(K),C(K1) 2060 121 IF(NBASE .EQ. 0) GO TO 123 2070 DO 122 I = 1,NBASE 2080 K = NDAT + LT + 2 * NPER + I 2090 122 WRITE(IPR,409) K,C(K) 2100 123 CONTINUE 2110 C COMPUTE RESIDUAL TIME SERIES 'F' 2120 DO 125 I = 1,NF 2130 PM = 0. 2140 DO 124 J = 1,M 2150 124 PM = PM + C(J) * BASE(J,T(I),DAT,NDAT,LT,FF,NPER) 2160 125 F(I) = F(I) - PM 2170 C PLOT RESIDUAL TIME SERIES WRITE (IPR, 415) 415 FORMAT(///,2X, 'PLOT OF RESIDUAL TIME SERIES',////) CALL TSPLOT (F,T,NF,X,SCALE) C COMPUTE QUADRATIC NORM OF 'F' 2190 126 FNORM = 0. 2200 DO 127 I = 1,NF 2210 127 FNORM = FNORM + F(I) * F(I) 2220 WRITE(IPR,411) FNORM 2230 SPP = SQRT(FNORM/NF) WRITE (IPR,600) SPP C PLOT SMOOTHED RESIDUAL TIME SERIES NFS=NF-2 FS(1)=0.0 FS(2)=0.0 DO 511 I=3,NFS FS(I) =0.1*(F(I-2)+F(I+2)+2.0*(F(I-1)+F(I+1))+4.0*F(I)) 511 CONTINUE WRITE (IPR, 416) 416 FORMAT (///,2X,'PLOT OF SMOOTHED RESIDUAL TIME SERIES',////) CALL TSPLOT (FS,T,NFS,X,SCALE) C COMPUTE QUADRATIC NORM OF SMOOTHED RESIDUAL TIME SERIES FNOR = 0.0 DO 502 I = 3,NFS 502 FNOR = FNOR + FS(I)*FS(I) WRITE (IPR, 417 ) FNOR SPP = SQRT(FNOR/(NF-4*NINT)) WRITE (IPR,600) SPP 417 FORMAT(///,9X,'SMOOTHED RESIDUAL SERIES QUADRATIC NORM = ',E9.4,/) C ASSUME RESIDUAL F CONSISTS OF ROUNDOFF NOISE WHEN 2240 C RMS OF RESIDUAL F .LT. 0.1% OF MID-RANGE OF ORIGINAL F 2250 IF(FNORM .LT. XNF * (0.001 * FMID)**2) GO TO 340 2260 C 2270 C ENTRY FOR SUBSEQUENT SPECTRAL BANDS TO BE COMPUTED 2280 C 2290 200 CONTINUE 2300 C DEFAULT VALUE OF PL IS CT 2310 IF(PL .EQ. 0.) PL = CT 2320 C SHORTEST PERMITTED PERIOD IS 2 STEPS 2330 IF(PS .LT. 2.*STEP) PS = 2. * STEP 2340 OM1 = 2. * PI / PL 2350 C DEFAULT VALUE OF NW IF INPUTTED AS ZERO 2360 IF(NW .EQ. 0) NW = (2. * PI / PS - OM1) * CT / 2. 2370 DOM = (2. * PI / PS - OM1) / (NW - 1.) 2380 WRITE(IPR,412) IB,NW,PL,PS 2390 DO 700 I = 1,100 700 IPLOT(I) = IBLANK C 2400 C LOOP OVER SPECTRAL VALUES TO BE COMPUTED 2410 C 2420 DO 217 I = 1,NW 2430 OMEGA = OM1 + DOM * (I - 1.) 2440 FCOS = 0. 2450 FSIN = 0. 2460 S1 = 0. 2470 S2 = 0. 2480 CU = 0. 2490 CV = 0. 2500 CW = 0. 2510 C CLEAR VECTORS 'U','V' 2520 IF(M .EQ. 0) GO TO 202 2530 DO 201 J = 1,M 2540 U(J) = 0. 2550 201 V(J) = 0. 2560 C COMPUTE SCALAR PRODUCTS 'FCOS' AND 'FSIN' 2570 202 DO 204 J = 1,NF 2580 WT = OMEGA * T(J) 2590 COSWT = COS(WT) 2600 SINWT = SIN(WT) 2610 FCOS = FCOS + F(J) * COSWT 2620 FSIN = FSIN + F(J) * SINWT 2630 C COMPUTE USER-DEFINED SCALAR PRODUCTS IN 'U','V' 2640 IF(NBASE .EQ. 0) GO TO 204 2650 DO 203 L = 1,NBASE 2660 K = NDAT + LT + 2 * NPER + L 2670 AA = BASE(K,T(J),DAT,NDAT,LT,FF,NPER) 2680 U(K) = U(K) + AA * COSWT 2690 203 V(K) = V(K) + AA * SINWT 2700 204 CONTINUE 2710 C COMPUTE TRIG(Q) 2720 Q = OMEGA * STEP / 2. 2730 SQ = SIN(Q) 2740 CQ = COS(Q) 2750 SQIN = 1. / SQ 2760 CQIN = 1. / CQ 2770 COTQ = CQ * SQIN 2780 C INITIALIZE FORCED FREQUENCY PARAMETERS 2790 IF(NPER .EQ. 0) GO TO 206 2800 DO 205 J = 1,NPER 2810 X = SM(J) * CQ 2820 Y = CM(J) * SQ 2830 SQM(J) = X - Y 2840 C ENSURE SQM .GE. DIVIDE 2850 IF(ABS(SQM(J)) .LT. DIVIDE) SQM(J) = 2860 * SIGN(DIVIDE,SQM(J)) 2870 SQP(J) = X + Y 2880 CP(J) = 0. 2890 SP(J) = 0. 2900 CMM(J) = 0. 2910 205 SMM(J) = 0. 2920 C COMPUTE TRIG(N*Q), TRIG(L*Q) AND THEIR PRODUCTS 2930 206 DO 210 J = 1,NINT 2940 QN = XN(J) * Q 2950 QL = XL(J) * Q 2960 SNQ = SIN(QN) 2970 CNQ = COS(QN) 2980 SLQ = SIN(QL) 2990 CLQ = COS(QL) 3000 X1 = SNQ * SLQ 3010 X2 = CNQ * CLQ 3020 X3 = SNQ * CLQ 3030 X4 = CNQ * SLQ 3040 C SUMMATIONS FOR SCALAR PRODUCTS 'S1', 'S2' , 'S3' 3050 S1 = S1 + X2 * X3 - X1 * X4 3060 S2 = S2 + X3 * X4 3070 C COMPUTE DATUM BIAS SCALAR PRODUCTS IN 'U','V' 3080 IF(M .EQ. 0) GO TO 210 3090 IF(NDAT .EQ. 0) GO TO 207 3100 K = INT(J) 3110 U(K) = U(K) + SQIN * X3 3120 V(K) = V(K) + SQIN * X1 3130 C SUMMATIONS FOR LINEAR TREND SCALAR PRODUCTS IN 'U','V' 3140 207 IF(LT .EQ. 0) GO TO 208 3150 U(NDAT1) = U(NDAT1) - COTQ*X1+XN(J)*X4+XL(J)*X3 3160 V(NDAT1) = V(NDAT1) + COTQ*X3-XN(J)*X2+XL(J)*X1 3170 C COMPUTE FORCED FREQUENCY PARAMETERS 3180 208 IF(NPER .EQ. 0) GO TO 210 3190 DO 209 K = 1,NPER 3200 X = SNM(J,K) * CNQ 3210 Y = CNM(J,K) * SNQ 3220 SNQM = X - Y 3230 SNQP = X + Y 3240 X = SLM(J,K) * CLQ 3250 Y = CLM(J,K) * SLQ 3260 SLQM = X - Y 3270 SLQP = X + Y 3280 X = CLM(J,K) * CLQ 3290 Y = SLM(J,K) * SLQ 3300 CLQM = X + Y 3310 CLQP = X - Y 3320 CMM(K) = CMM(K) + SNQM * CLQM 3330 SMM(K) = SMM(K) + SNQM * SLQM 3340 CP(K) = CP(K) + SNQP * CLQP 3350 209 SP(K) = SP(K) + SNQP * SLQP 3360 210 CONTINUE 3370 C COMPUTE SCALAR PRODUCTS 'S1','S2', AND 'S3' 3380 AA = S1 * SQIN * CQIN 3390 S1 = (XNF + AA) / 2. 3400 S3 = (XNF - AA) / 2. 3410 S2 = S2 * SQIN * CQIN 3420 IF(M .EQ. 0) GO TO 215 3430 C COMPUTE LINEAR TREND SCALAR PRODUCT IN 'U','V' 3440 IF(LT .EQ. 0) GO TO 211 3450 AA = SQIN * STEP / 2. 3460 U(NDAT1) = U(NDAT1) * AA 3470 V(NDAT1) = V(NDAT1) * AA 3480 C COMPUTE FORCED FREQUENCY SCALAR PRODUCTS IN 'U','V' 3490 211 IF(NPER .EQ. 0) GO TO 213 3500 DO 212 J = 1,NPER 3510 SQMIN = 0.5 / SQM(J) 3520 SQPIN = 0.5 / SQP(J) 3530 K = NDAT + LT + 2 * J - 1 3540 U(K ) = SQMIN * CMM(J) + SQPIN * CP(J) 3550 V(K ) = -SQMIN * SMM(J) + SQPIN * SP(J) 3560 U(K+1) = SQMIN * SMM(J) + SQPIN * SP(J) 3570 212 V(K+1) = SQMIN * CMM(J) - SQPIN * CP(J) 3580 C COMPUTE BILINEAR FORMS 'CU','CV','CW' 3590 213 DO 214 J = 1,M 3600 DO 214 K = 1,M 3610 CU = CU + A(J,K) * U(J) * U(K) 3620 CV = CV + A(J,K) * V(J) * V(K) 3630 214 CW = CW + A(J,K) * U(J) * V(K) 3640 C COMPUTE PERIOD AND SPECTRAL VALUE (PERCENTAGE VARIANCE) 3650 215 PERIOD = 2. * PI / OMEGA 3660 SIG = 0. 3670 D = (S1-CU) * (S3-CV) - (S2-CW) * (S2-CW) 3680 IF(ABS(D) .LT. DIVIDE) GO TO 216 3690 SIG = ((S3-CV)*FCOS*FCOS -2.*(S2-CW)*FCOS*FSIN + 3700 * (S1-CU)*FSIN*FSIN) / (D*FNORM) 3710 SIG = SIG * 100. 3720 C WRITE PERIOD AND SPECTRAL VALUE 3730 216 KSIG = SIG 3740 IF(KSIG .LT. 1) KSIG = 1 3750 IF(KSIG .GT. 100) KSIG = 100 3760 IPLOT(KSIG) = ISTAR 3770 WRITE(IPR,413) PERIOD,SIG,IPLOT 3780 IPLOT(KSIG) = IBLANK 3790 217 CONTINUE 3800 Q = 1.0/(NF - M) WRITE (IPR,418)Q 418 FORMAT (///,2X,'LEVEL OF RANDOM NOISE SPECTRUM =',F6.4) RETURN 3810 C 3820 C ERROR STOPS 3830 310 WRITE(IPR,311) 3840 311 FORMAT(9X,26HINPUT ERROR-DAT(1).NE.T(1)) 3850 STOP 3860 320 WRITE(IPR,321) 3870 321 FORMAT(9X,34HINPUT ERROR-STEPNUMBER SEQUENCE AT,2F9.0) 3880 STOP 3890 330 WRITE(IPR,331) 3900 331 FORMAT(9X,40HERROR-NORMAL EQUATION MATRIX IS SINGULAR) 3910 STOP 3920 340 WRITE(IPR,341) 3930 341 FORMAT(9X,40HRESIDUAL TIME SERIES=ROUNDOFF NOISE.STOP) 3940 STOP 3950 350 WRITE(IPR,351) PER(I) 3960 351 FORMAT(9X,35HINPUT ERROR-FORCED PERIOD TOO SHORT,F9.5) 3970 STOP 3980 C FORMAT STATEMENTS 3990 401 FORMAT(1H1,9X,11HTIME SERIES//5X,1HI,8X,4HT(I),8X, 4000 * 4HF(I),10X,24HCONSTANT TIME INCREMENT=,E12.4) 4010 402 FORMAT(65X,6HGAP OF,E14.6,6H STEPS) 4020 403 FORMAT(2X,I4,2X,E11.4,2X,E11.4,33(1H-) , 18HBEGINNING OF DATUM, &I5,39(1H-)) 404 FORMAT(2X,I4,2E12.4,100A1) 4040 405 FORMAT(33H1 SOLUTION FOR KNOWN CONSTITUENTS/19X, 4050 * 5HDATUM,8X,6HLINEAR,8X,6HFORCED,8X,6HCOSINE,10X, 4060 * 4HSINE,10X,4HUSER,5X,5HTOTAL/20X,4HBIAS,9X,5HTREND, 4070 * 8X,6HPERIOD,35X,7HDEFINED//2X,7HNUMBER ,2I14,14X, 4080 * 3I14,I9//2X,9HAMPLITUDE//) 4090 406 FORMAT(7X,I3,E14.6) 4100 407 FORMAT(7X,I3,14X,E14.6) 4110 408 FORMAT(3X,I3,1H-,I3,28X,F14.6,2E14.6) 4120 409 FORMAT(7X,I3,70X,E14.6) 4130 410 FORMAT (///9X,20HRESIDUAL TIME SERIES/99(10E10.2/)) 411 FORMAT (///9X,38HRESIDUAL TIME SERIES QUADRATIC NORM = ,E9.4,/) 412 FORMAT(1H1,10X,13HSPECTRAL BAND,I5/10X,I5, 4160 * 24H SPECTRAL VALUES BETWEEN,2F14.6//9X,6HPERIOD,6X, 4170 * 19HPERCENTAGE VARIANCE) 4180 413 FORMAT(5X,F10.5,F12.3,100A1) 4190 C 4200 END 4210 FUNCTION BASE(I,T,DAT,NDAT,LT,FF,NPER) 4220 C 4230 C COMPUTE BASE FUNCTIONS FOR NORMAL EQUATIONS 4240 C 4250 COMMON E(5,700) DIMENSION DAT(12),FF(5) IF(I - NDAT) 1,1,2 4270 C BASE FUNCTION = DATUM BIAS 4280 1 BASE = 0. 4290 IF(T .GE. DAT(I) .AND. T .LT. DAT(I+1)) BASE = 1. 4300 RETURN 4310 2 IF(I - NDAT - LT) 3,3,4 4320 C BASE FUNCTION = LINEAR TREND 4330 3 BASE = T 4340 RETURN 4350 4 IF(I - NDAT - LT - 2*NPER) 5,5,8 4360 C BASE FUNCTION = TRIGONOMETRIC FUNCTION OF FORCED FREQUENCY 4370 5 L = I - NDAT - LT + 1 4380 IF(MOD(L,2)) 6,6,7 4390 6 L = L / 2 4400 BASE = COS(FF(L) * T) 4410 RETURN 4420 7 L = L / 2 4430 BASE = SIN(FF(L) * T) 4440 RETURN 4450 C BASE FUNCTIONS = TEMPERATURE, PRESSURE, RIVER DISCHARGES 8 IT = T + 0.1 J = I - NDAT - LT - 2*NPER BASE = E(J,IT) RETURN 4490 END 4500 C C******************************************************************************* C C SUBROUTINE 'T S P L O T': C FOR COMPUTING AND PRINTING THE MAXIMUM, C MINIMUM AND MEAN VALUES OF A GIVEN TIME SERIES. ALSO IT PLOTS THE C CORRESPONDING TIME SERIES ALONG WITH THE FUNCTIONAL VALUES AND C THE SCALE OF THE PLOT. C C******************************************************************************* C C MODIFIED VERSION FOR THE M.S.L. PLOTS. C SEPT. 25, 1975. C C******************************************************************************* SUBROUTINE TSPLOT(F,T,N,X,SCALE) DIMENSION F(N),T(N) DIMENSION IPLOT(30),IFRAME(30) DATA IBLANK/' '/,ISTAR/'*'/,IPLUS/'+'/,IDASH/'|'/ INTEGER K C C SPECIFY THE DESIRED LATERAL SPACE FOR THE PLOT ,IN TERMS OF THE C EQUIVALENT NUMBER OF CHARACTERS ON THE LINE PRINTER 'NCOL'. C NCOL=30 C C SET-UP VALUES OF '+' FOR THE ELEMENTS OF THE ARRAY 'IFRAME', C TO BE USED IN FRAMING THE PLOT. C DO 9 I=1,NCOL IFRAME(I)=IPLUS 9 CONTINUE C C COMPUTE THE MEAN, MAXIMUM AND MINIMUM FUNCTIONAL VALUES FOR THE C GIVEN TIME-SERIES (REQUIRED TO BE PLOTTED). C C NOTE:- THE ZERO VALUES WHICH REPRESENT MISSING DATA ARE EXCLUDED C FROM THE COMPUTATIONS. C NZEROF=0 SUMF=0.0 DO 1 I=1,N IF(F(I).EQ.0.0) NZEROF=NZEROF+1 SUMF=SUMF+F(I) 1 CONTINUE NSUMF=N-NZEROF FBAR=SUMF/FLOAT(NSUMF) C FMAX=0.0 FMIN=1.0E 10 DO 2 I=1,N FMAX=AMAX1(FMAX,F(I)) IF(F(I).EQ.0.0) GO TO 2 FMIN=AMIN1(FMIN,F(I)) 2 CONTINUE C C INITIALIZE THE PLOTTING ARRAY. C DO 10 I=1,NCOL IPLOT(I)=IBLANK 10 CONTINUE IPLOT(1)=IPLUS IPLOT(NCOL)=IPLUS C Y=FMIN*X-1.0 KBAR=FBAR*X-Y IPLOT(KBAR)=IDASH C C PRINT-OUT THE MEAN, MAX. & MIN. VALUES OF THE TIME SERIES. C WRITE(6,60) FMAX,FMIN,FBAR 60 FORMAT(5X,'MAX. FUNCTIONAL VALUE IN THE SERIES = ',E11.4,/,5X,'MIN 1. FUNCTIONAL VALUE IN THE SERIES = ',E11.4,/,5X,'MEAN FUNCTIONAL V 2ALUE IN THE SERIES = ',E11.4/) C C PRINT-OUT THE SCALE OF THE PLOT, BY TAKING INTO ACCOUNT THE LINE C PRINTER # 1403 (AT U.N.B.), WHICH PRINTS 10 CHARACTERS PER INCH. C WRITE(6,70) SCALE 70 FORMAT(3X,'NO.',7X,'T',10X,'VALUE',6X,'SCALE OF THE PLOT IS: 1 INC 1H REPRESENTS',E10.3,' UNITS') C C PLOT THE GIVEN TIME SERIES, WITH ITS FUNCTIONAL VALUES. C WRITE(6,40) (IFRAME(J),J=1,NCOL) K=IFIX(T(1)-0.6) DO 30 I=1,N IF (I.EQ.1) GO TO 4 IF ((T(I) - T(I-1)).GT.1.1) GO TO 3 4 K = K + 1 IF (K. GT. 12) K = 1 KF=F(I)*X-Y IF(KF.LT.1) KF=1 IF(KF.GT.NCOL) KF=NCOL IPLOT(KF)=ISTAR WRITE(6,20) K,T(I),F(I),(IPLOT(J),J = 1,NCOL) C C RE-INITIALIZE THE PLOTTING ARRAY. C IPLOT(KF)=IBLANK IPLOT(1)=IPLUS IPLOT(NCOL)=IPLUS IPLOT(KBAR)=IDASH 30 CONTINUE WRITE(6,40) (IFRAME(J),J=1,NCOL) 6 FORMAT ('1') 20 FORMAT(2X,I4,2X,E11.4,2X,E11.4,3X,30A1) 40 FORMAT(35X,30A1) RETURN 3 NN = T(I) - T(I-1) - 0.9 TT = T(I-1) DO 5 JJ = 1,NN K = K + 1 IF (K.GT.12)K = 1 TT = TT + 1 FF = 0.0 KF = 1 IPLOT(KF) = ISTAR WRITE (6,20) K,TT,FF,(IPLOT(J),J = 1,NCOL) IPLOT (KF) = IBLANK IPLOT (1) = IPLUS IPLOT (NCOL) = IPLUS IPLOT(KBAR) = IDASH 5 CONTINUE GO TO 4 END