IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AIR(600),WATER(600),DATA(600,2),VAL(12),C(600),P(600), @ F(600),TSNAM(2),ATA(3,4),X(3),KTYP(2) EQUIVALENCE (AIR(1),DATA(1,1)),(WATER(1),DATA(1,2)) DATA PI/3.141592653589793D0/,TSNAM/'FIRST','SECOND'/,IYES/'Y'/ HPI=PI/2.D0 108 READ (5,501,END=101) IYS,MS,IYF,MF,KTYP,KRES,EH,EDT,TOLH,TOLDT 501 FORMAT (I4,1X,I2,3X,I4,1X,I2,3X,2A1,1X,A1,6X,4F10.4) NY=IYF-IYS+1 IF (NY.GT.0.AND.NY.LE.50) GO TO 102 WRITE (6,601) IYS,MS,IYF,MF 601 FORMAT (//' ERROR IN START/FINISH DATES:',4I5) STOP 0001 102 DO 401 ID=1,2 I=1 DO 402 IY=IYS,IYF READ (5,502,END=103) KODE,JY,VAL 502 FORMAT (A1,6X,I4,12F5.1) IF (KODE.EQ.KTYP(ID).AND.JY.EQ.IY) GO TO 104 WRITE (6,602) KODE,JY 602 FORMAT (//' CARD OUT OF ORDER:',A1,I6) STOP 0002 103 WRITE (6,603) 603 FORMAT (//' INSUFFICIENT DATA CARDS') STOP 0003 104 JS=1 IF (IY.EQ.IYS) JS=MS JF=12 IF (IY.EQ.IYF) JF=MF DO 403 J=JS,JF DATA(I,ID)=VAL(J) I=I+1 403 CONTINUE 402 CONTINUE 401 CONTINUE NN=I-1 IF (MOD(NN,2).EQ.0) GO TO 105 WRITE (6,604) NN 604 FORMAT (//' WARNING: INPUT DATA TRUNCATED TO AN EVEN NUMBER OF PO @INTS =',I4) 105 N=NN/2 NN=N*2 NMO=N-1 WRITE (6,605) IYS,MS,IYF,MF,NN 605 FORMAT ('1LEAST SQUARES LINEAR TRANSFORM BETWEEN TWO TIME SERIES'/ @ 1X,54('=')//5X,'START DATE =',I5,':',I2,5X,'END DATE =',I5,':', @ I2,5X,'NUMBER OF POINTS =',I4) DO 405 ID=1,2 ITS=3-ID CALL FFS (DATA(1,ITS),NN,C) CALL POLAR (C,NN,P) WRITE (6,606) TSNAM(ITS),KTYP(ITS),C(1),P(1) 606 FORMAT (////' FOURIER TRANSFORM OF ',A6,' TIME SERIES -- TYPE CODE @ = ',A1/1X,56('-')/ @ //1X,82('-')/1X,'HARMONIC PERIOD STANDARD COEFFICIENTS @ AMPLITUDE PHASE ASC. NODE'/11X,'MONTHS COSINE',10X, @ 'SINE',21X,'RADIANS MONTHS'/1X,82('-')//5X,'0',13X,D13.6, @ 17X,F10.6) DO 404 I=1,NMO PER=NN/DFLOAT(I) FREQ=PI*2.D0/PER IC=I*2 IS=IC+1 ANODE=(P(IS)-HPI)/FREQ IF (ANODE.LT.0.D0) ANODE=ANODE+PER WRITE (6,607) I,PER,C(IC),C(IS),P(IC),P(IS),ANODE 607 FORMAT (1X,I5,F11.3,2D15.6,3F12.6) 404 CONTINUE PER=2.D0 WRITE (6,608) N,PER,C(NN),P(NN) 608 FORMAT (1X,I5,F11.3,D15.6,17X,F10.6) WRITE (6,609) 609 FORMAT (/1X,80('-')) 405 CONTINUE C *** BEGIN ITERATION LOOP ON SOLUTION *** ITER=1 C *** FORM NORMAL EQUATIONS *** 107 CALL NORM (C,NN,EDT,EH,WATER,ATA,TSQ) WRITE (6,610) ITER,EH,EDT,((ATA(I,J),J=1,4),I=1,3) 610 FORMAT ('1SOLUTION -- ITERATION #',I3/1X,26('=')//5X,'ESTIMATE OF & @HEAT TRANSFER COEFF. =',F10.6//5X,'ESTIMATE OF PHASE SHIFT',10X, @ '=',F10.6,' MONTHS'////5X,'NORMAL EQUATIONS:'/5X,16('-')///12X, @ 'DELTA(H) DELTA(T)',11X,'K',12X,'CONST.'//(8X,3F14.6,F17.6)) CALL CHOLD (ATA,3,3,DET,&106) WRITE (6,611) ((ATA(I,J),J=1,I),I=1,3),DET 611 FORMAT (////5X,'WEIGHT COEFF. MATRIX:'/5X,20('-')///8X,F14.6/ @ 8X,2F14.6/8X,3F14.6,14X,'DETERMINANT =',D14.6) C *** SOLVE BY MATRIX MULTIPLICATION *** DO 407 I=1,3 X(I)=0.D0 DO 407 J=1,3 X(I)=X(I)-ATA(I,J)*ATA(J,4) 407 CONTINUE C ***DETERMINE MINIMUM AND HENCE VARIANCE FACTOR *** DO 406 I=1,3 TSQ=TSQ+ATA(I,4)*X(I) 406 CONTINUE TSQ=TSQ/(NN-3) WRITE (6,612) X,TSQ 612 FORMAT (////5X,'SOLUTION:'/5X,8('-')///8X,3F14.6,10X,'VARIANCE FAC @TOR =',F11.6) EH=EH+X(1) EDT=EDT+X(2) EK=X(3) ITER=ITER+1 IF (ITER.LE.10) GO TO 109 WRITE (6,616) 616 FORMAT (///' SOLUTION NOT BELOW ITERATION TOLERANCE AFTER TEN ITER @ATIONS -- SOLUTION ABANDONED') GO TO 108 109 IF (TOLH.EQ.0.D0) TOLH=0.1*DABS(X(1)) IF (TOLDT.EQ.0.D0) TOLDT=0.1*DABS(X(2)) WRITE (6,617) TOLH,TOLDT 617 FORMAT (////5X,'ITERATION TOLERANCE:'/5X,19('-')///8X,2F14.6) IF (DABS(X(1)).GT.TOLH.OR.DABS(X(2)).GT.TOLDT) GO TO 107 C *** PRINT SOLUTION FOR TRANSFORM PARAMETERS *** SDEDT=DSQRT(ATA(2,2)*TSQ) SDEH=DSQRT(ATA(1,1)*TSQ) SDK=DSQRT(ATA(3,3)*TSQ) WRITE (6,613) EDT,SDEDT,EH,SDEH,EK,SDK 613 FORMAT ('1TRANSFORMATION PARAMETERS'/1X,25('=')///5X,'PHASE SHIFT @=',14X,F10.6,5X,'S.D. =',F10.6//5X,'HEAT TRANSFER COEFFICIENT =', @ F10.6,5X,'S.D. =',F10.6//5X,'MEAN TEMP. DIFFERENCE = ',F10.6, @ 5X,'S.D. =',F10.6) C *** COMPUTE RESIDUALS *** CALL FINT (C,NN,EDT,F) SR=0.D0 SSR=0.D0 RMAX=0.D0 DO 408 I=1,NN IF (I.EQ.1.AND.EDT.LT.0.D0) GO TO 111 IF (I.EQ.NN.AND.EDT.GT.0.D0) GOTO 111 C(I)=F(I)*EH+EK P(I)=WATER(I)-C(I) SR=SR+P(I) SSR=SSR+P(I)**2 RMAX=DMAX1(RMAX,DABS(P(I))) GO TO 408 111 C(I)=0.D0 P(I)=0.D0 F(I)=0.D0 408 CONTINUE IF (KRES.NE.IYES) GO TO 110 DO 409 I=1,N IF (MOD(I,50).EQ.1) WRITE (6,618) (((KTYP(K),L=1,2),K=1,2),M=1,2) 618 FORMAT ('1',2(2X,53('-'),8X)/' ',2(2X,'PT.',5X,A1,5X,A1,'(T+DT) @ ',A1,'(O)',5X,A1,'(C)',7X,'RESIDUAL',8X)/' ',2(2X,53('-'),8X)//) IPA=I IPB=I+N WRITE (6,619) IPA,AIR(IPA),F(IPA),WATER(IPA),C(IPA),P(IPA), @ IPB,AIR(IPB),F(IPB),WATER(IPB),C(IPB),P(IPB) 619 FORMAT (2(3X,I3,2(F7.1,F12.6),F12.6,7X)) 409 CONTINUE 110 SR=SR/(NN-1) SSR=DSQRT(SSR /(NN-4)) WRITE (6,620) SR,SSR,RMAX 620 FORMAT ('1ANALYSIS OF RESIDUALS'/1X,21('=')///5X,'MEAN RESIDUAL = @',F10.6//5X,'STANDARD ERROR =',F10.6//5X,'MAX. RESIDUAL =',F10.6) CALL HISTO (P,NN,SR,SSR) GO TO 108 106 WRITE (6,615) ATA,DET 615 FORMAT (/////' MATRIX IS SINGULAR -- SOLUTION ABANDONED'// @ (4D20.6)) GO TO 108 101 WRITE (6,614) 614 FORMAT (/////' ALL DATA PROCESSED'/1X,18('-')) STOP END SUBROUTINE POLAR (A,NN,P) C C******************************************************************************* C C SUBROUTINE POLAR C ***************** C C DOUBLE PRECISION ROUTINE TO CONVERT FOURIER COEFFICIENTS TO POLAR FORM C I.E. AMPLITUDE (C) AND PHASE ANGLE (D). C C INPUT C ----- C C A(NN) REAL*8 ARRAY OF FOURIER COEFFICIENTS IN ORDER: C A(0),A(1),B(1),A(2),B(2), ... A(N-1),B(N-1),A(N) C WHERE N=NN/2 C C NN INT*4 SIZE OF ARRAY A C C OUTPUT C ------ C C P(NN) REAL*8 ARRAY OF COEFFICIENTS IN POLAR FORM IN ORDER: C A(0),C(1),D(1),C(2),D(2), ... C(N-1),D(N-1),A(N) C C******************************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(1),P(1) P(1)=A(1) M=NN-2 DO 401 I=2,M,2 P(I)=DSQRT(A(I)**2+A(I+1)**2) P(I+1)=DATAN2(A(I+1),A(I)) 401 CONTINUE P(NN)=A(NN) RETURN END SUBROUTINE NORM (A,NN,EDT,EH,TAU,ATA,TSQ) C C******************************************************************************* C C SUBROUTINE NORM C *************** C C DOUBLE PRECISION ROUTINE TO FORM OBSEVATION EQUATIONS AND HENCE NORMAL C EQUATIONS. C C******************************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(1),TAU(1),ATA(3,4),C(4) DATA PI /3.141592653589793D0/ N=NN/2 OM=PI/N C(3)=1.D0 DO 400 I=1,3 DO 400 J=1,4 ATA(I,J)=0.D0 400 CONTINUE TSQ=0.D0 DO 401 I=1,NN IF (I.EQ.1.AND.EDT.LT.0.D0) GO TO 401 IF (I.EQ.NN.AND.EDT.GT.0.D0) GO TO 401 X=(I-1+EDT)*OM C(1)=A(1)/2.D0 C(2)=0.D0 INIT=1 DO 402 J=1,N K=J*2 CALL MAF (X,SINJX,COSJX,INIT) C(1)=C(1)+A(K)*COSJX Y=-A(K)*SINJX IF (J.EQ.N) GO TO 501 C(1)=C(1)+A(K+1)*SINJX Y=Y+A(K+1)*COSJX 501 C(2)=C(2)+Y*J*OM 402 CONTINUE C(2)=C(2)*EH C(4)=C(1)*EH-TAU(I) CALL SQUARE (C,ATA,3,TSQ) 401 CONTINUE RETURN END SUBROUTINE SQUARE (OBSN,ANORM,NR,TSQ) C C******************************************************************************* C C C******************************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION OBSN(1),ANORM(NR,NR) NC=NR+1 DO 401 I=1,NR DO 402 J=1,NC IF (J.LT.I) GO TO 501 ANORM(I,J)=ANORM(I,J)+OBSN(I)*OBSN(J) GO TO 402 501 ANORM(I,J)=ANORM(J,I) 402 CONTINUE 401 CONTINUE TSQ=TSQ+OBSN(NC)**2 RETURN END SUBROUTINE CHOLD(A,IRDA,NA,DETA,*) 0420 C 0430 C MATRIX INVERSION USING CHOLESKI DECOMPOSITION 0440 C 0450 C INPUT ARGUMENTS 0460 C A = ARRAY CONTAINING POSITIVE DEFINITE SYMMETRIC INPUT MATRIX 0470 C IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX 0480 C NA = SIZE OF INPUT MATRIX 0490 C OUTPUT ARGUMENTS 0500 C A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED) 0510 C DETA = DETERMINANT OF INPUT MATRIX 0520 C * = ERROR RETURN (TAKEN IF NA .LT. 1 OR IF DETA .LT. SING) 0530 C 0540 DOUBLE PRECISION A,DETA,SUM,SQRT,DSQRT,ABS,DABS,SING 0550 DIMENSION A(IRDA,NA) 0560 SQRT(SUM) = DSQRT(SUM) 0570 ABS(DETA) = DABS(DETA) 0580 DATA SING/1D-10/ 0590 C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 0600 IF(NA .LT. 1) GO TO 18 0610 DETA = A(1,1) 0620 A(1,1) = SQRT(A(1,1)) 0630 IF(NA .EQ. 1) GO TO 6 0640 DO 1 I = 2,NA 0650 1 A(I,1) = A(I,1) / A(1,1) 0660 DO 5 J = 2,NA 0670 SUM = 0. 0680 J1 = J - 1 0690 DO 2 K = 1,J1 0700 2 SUM = SUM + A(J,K) ** 2 0710 DETA = DETA * (A(J,J) - SUM) 0720 A(J,J) = SQRT(A(J,J) - SUM) 0730 IF(J .EQ. NA) GO TO 5 0740 J2 = J + 1 0750 DO 4 I = J2,NA 0760 SUM = 0. 0770 DO 3 K = 1,J1 0780 3 SUM = SUM + A(I,K) * A(J,K) 0790 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 0800 5 CONTINUE 0810 6 IF(ABS(DETA) .LT. SING) GO TO 16 0820 C INVERSION OF LOWER TRIANGULAR MATRIX 0830 DO 7 I = 1,NA 0840 7 A(I,I) = 1. / A(I,I) 0850 IF(NA .EQ. 1) GO TO 10 0860 N1 = NA - 1 0870 DO 9 J = 1,N1 0880 J2 = J + 1 0890 DO 9 I = J2,NA 0900 SUM = 0. 0910 I1 = I - 1 0920 DO 8 K = J,I1 0930 8 SUM = SUM + A(I,K) * A(K,J) 0940 9 A(I,J) = - A(I,I) * SUM 0950 C CONSTRUCTION OF INVERSE OF INPUT MATRIX 0960 10 DO 15 J = 1,NA 0970 IF(J .EQ. 1) GO TO 12 0980 J1 = J - 1 0990 DO 11 I = 1,J1 1000 11 A(I,J) = A(J,I) 1010 12 DO 14 I = J,NA 1020 SUM = 0. 1030 DO 13 K = I,NA 1040 13 SUM = SUM + A(K,I) * A(K,J) 1050 14 A(I,J) = SUM 1060 15 CONTINUE 1070 RETURN 1080 16 WRITE(6,17) DETA 1090 17 FORMAT(10X, 'SINGULAR MATRIX IN CHOLD. DET =',E20.5) 1100 RETURN 1 1110 18 WRITE(6,19) 1120 19 FORMAT(10X,'MATRIX OF DIMENSION ZERO IN CHOLD') 1130 RETURN 1 1140 END 1150 SUBROUTINE HISTO (A,NN,AMEAN,SD) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(1),NBIN(40),LINE(100) DATA IBLK,MM /' ','M'/ DO 401 I=1,40 NBIN(I)=0 401 CONTINUE DO 402 I=1,NN X=(A(I)-AMEAN)/SD IBIN=X/0.2D0 IBIN=IBIN+21 IF (IBIN.LT.1.OR.IBIN.GT.40) GO TO 402 NBIN(IBIN)=NBIN(IBIN)+1 402 CONTINUE WRITE (6,601) 601 FORMAT (///20X,10('+',9('-'))) DO 403 I=1,40 XL=(I-21)*0.2D0 XU=XL+0.2D0 DO 404 J=1,100 LINE(J)=IBLK IF (NBIN(I).GT.J-1) LINE(J)=MM 404 CONTINUE WRITE (6,602) XL,XU,LINE 602 FORMAT (6X,F4.1,' TO',F5.1,' |',100A1) 403 CONTINUE RETURN END