C C THIS IS PART 'A' OF THE PROGRAM. EACH RECORD CONTAINS LATITUDE AND C LONGITUDE OF THE STATION AND THE YEARLY MEANS. ALL YEARLY MEANS ARE C KEPT IN THE MATRIX 'A' . THE LATITUDES AND LONGITUDES ARE KEPT IN THE C VECTORS---> PHID,PHIM,LAMD,LAMM. C IMPLICIT REAL (A-H,O-Z) REAL*4 A(71,3)/213*8888./,AA(71,3)/213*8888./,DEL(71)/71*0.0/, #PHID(10),PHIM(10),LAMD(10),LAMM(10),B(100) N=3 M=71 KK=0 1 KK=KK+1 READ(5,5,END=120) (B(L),L=1,16) 5 FORMAT(16F5.2) DO11 J=1,4 JJ=(J*16)+1 JK=JJ+15 READ(5,5,END=120) (B(L),L=JJ,JK) DO10 L=JJ,JK IF(B(L).EQ.999.)GOTO 14 NO=L 10 CONTINUE 11 CONTINUE 14 DO13 I=1,NO AA(I,KK)=B(I) IF(I.EQ.NO)GOTO 1 13 CONTINUE 120 CONTINUE DO15 I=1,N PHID(I)=AA(1,I) PHIM(I)=AA(2,I) LAMD(I)=AA(3,I) LAMM(I)=AA(4,I) 15 CONTINUE CALL MOUTS(AA,M,M,N) I=0 25 I=I+1 DO26 J=1,M A(J,I)=AA(J+4,I) IF(I.EQ.N.AND.J.EQ.(M-4)) GOTO27 IF(J.EQ.(M-4))GOTO 25 26 CONTINUE 27 CONTINUE DO30 I=1,N DO30 J=1,M AA(J,I)=0.0 30 CONTINUE C C THIS IS PART 'B' OF THE PROGRAM. IT IS DESIGNED SO THAT THE DATA FOR C ALL GUAGES IS IN THE CORRECT ORDER SO THAT THE 1974,1973,1972....MEANS C FOR ALL THE GUAGES ARE IN THE SAME POSITION IN EACH ROW. THIS DATA IS C IN IN THE MATRIX 'AA' . C C I=1 DO35 J=1,N KNT=0 DO36 K=1,M IF(A(K,J).GT.888.) KNT=KNT+1 IF(K.EQ.M) GOTO 40 36 CONTINUE 40 MM=M-KNT DO45 JJ=1,MM AA(JJ+KNT,I)=A(JJ,I) 45 CONTINUE I=I+1 35 CONTINUE C C C PART 'C' OF THE PROGRAM. CREATES ALL THE POSSIBLE PAIRS OF GUAGES AND C FORMS THEIR CORRESPONDING SERIES ' DELTA ' (WHICH EGUALS L2-L1),WHERE C L1,L2,L3...... ARE THE SERIES OF YEARLY MEANS FOR EACH GUAGE. C J=0 50 LL=0 J=J+1 IF(J.EQ.N)GOTO 65 55 K=1 LL=LL+1 IF((J+LL).GT.N)GOTO 50 60 DEL(K)=AA(K,J+LL)-AA(K,J) IF((AA(K,J+LL)).EQ.0.0.AND.(AA(K,J)).EQ.0.0) DEL(K)=0.0 IF((AA(K,J+LL)).EQ.0.0.AND.(AA(K,J)).NE.0.0) DEL(K)=0.0 IF((AA(K,J+LL)).NE.0.0.AND.(AA(K,J)).EQ.0.0) DEL(K)=0.0 K=K+1 IF(K.LE.M)GOTO 60 WRITE(6,67)PHID(J),PHIM(J),LAMD(J),LAMM(J) WRITE(6,66)PHID(J+LL),PHIM(J+LL),LAMD(J+LL),LAMM(J+LL) 66 FORMAT('+',T50,2(F5.0,3X,F5.2,3X),5(/)) 67 FORMAT(120('*'),10(/),1X,T22,'PHI(A)',T37,'LAMDA(A)',9X,'PHI(B)',T @67,'LAMDA(B)',//,1X,T15,2(F5.0,3X,F5.2,3X)) C C CALCULATION OF DISTANCE BETWEEN GUAGES (KM.) C ( USING MEAN RADIUS OF EARTH = 6371 KM. ) C PHIM(J)=PHIM(J)/60.0 PHIM(J+LL)=PHIM(J+LL)/60.0 LAMM(J)=LAMM(J)/60.0 LAMM(J+LL)=LAMM(J+LL)/60.0 RHID1=PHID(J)+PHIM(J) RHID2=PHID(J+LL)+PHIM(J+LL) RAMD1=LAMD(J)+LAMM(J) RAMD2=LAMD(J+LL)+LAMM(J+LL) DP=ABS(RHID1-RHID2) DL=ABS(RAMD1-RAMD2) DPL=SQRT(DP**2.0+DL**2.0) DIST=DPL*111.11 PHIM(J)=PHIM(J)*60.0 PHIM(J+LL)=PHIM(J+LL)*60.0 LAMM(J)=LAMM(J)*60.0 LAMM(J+LL)=LAMM(J+LL)*60.0 WRITE(6,122)DIST 122 FORMAT(1X,T20,'DIST=',F8.3) CALL TREND(2,1,DEL,DIST) GOTO 55 65 CONTINUE C C C END OF PART 'C' OF THE PROGRAM. C C STOP END C PROGRAM 'LSAPROX' LEAST SQUARES APPROXIMATION 3 C 4 C INPUT: CARD 1: JOB TITLE FORMAT(20A4) 5 C CARD 2: N=NUMBER OF BASE FUNCTIONS 6 C IW=WEIGHTING CODE ==> IW=0 READ WEIGHT ON CARD 3 7 C IW=1 COMPUTE WEIGHT IN FUNCTION 8 C WEIGHT 9 C FORMAT(2I4) 10 C CARD TYPE 3: X - VALUE OF BASE POINTS 11 C FUNCTIONAL VALUES OF X 12 C WEIGHT OF X (IF IW=1 LEAVE BLANK) 13 C FORMAT(3F20.0) 14 C NOTE: CARD TYPE 3 IS REPEATED ONCE FOR EACH DATA POINT 15 C CARD TYPE 4: MARKS END OF CARD TYPE 3 -- TYPE -99999.0 IN 16 C COLUMNS 1 - 8 17 C CARD TYPE 5: VALUES OF 'X' AT WHICH APPROXIMATED VALUE IS 18 C DESIRED 19 C NOTE: IF ONLY COEFFICIENTS ARE DESIRED CARD TYPE 5 IS NOT 20 C REQUIRED 21 C 22 C NOTES: MAXIMUM DEGREE OF APPROXIMATION FUNCTION IS 50 - TO INCREASE 23 C CHANGE 'NMAX' AND DIMENSION STATEMENT 24 C MAXIMUM NUMBER OF DATA POINTS IS 100 - TO INCREASE CHANGE 25 C 'MMAX' AND DIMENSION STATEMENT 26 C 27 C THE FOLLOWING CHECKS ARE PROVIDED 28 C (1) N.LT.1 29 C (2) N.GT.NMAX 30 C (3) M.GT.MMAX 31 C (4) N.GT.M 32 C (5) DETERMINANT OF GRAMM'S MATRIX .LT.SING IE. BASE FUNCTIONS ARE 33 C NOT LINEARLY INDEPENDENT 34 C 35 C THE USER MUST SUPPLY TWO DOUBLE PRECISION FUNCTIONS 36 C (1) FUNCTION WEIGHT -- DESCRIBES THE WEIGHT 37 C (2) FUNCTION PHI -- DESCRIBES THE BASE FUNCTIONS 38 C 39 C THE FOLLOWING MINIMUM REQUIREMENTS ARE SUPPLIED 40 C 41 C FUNCTION WEIGHT(X) 42 C REAL*8 WEIGHT,X 43 C WEIGHT=1.0D0 44 C RETURN 45 C END 46 C 47 C FUNCTION PHI(I,X) 48 C REAL*8 PHI,X 49 C PHI=1.0D0 50 C RETURN 51 C END 52 C 53 C IF THESE ARE NOT SUFFICIENT ADD 54 C 55 C (1) WEIGHT=F(X) 56 C (2) IF(I.EQ.1)PHI=F(X) 57 C IF(I.EQ.2)PHI=F(X) 58 C ETC 59 C 60 C 62 C 63 SUBROUTINE TREND(N,IW,FX,DIST) IMPLICIT REAL*8 (A-H,O-Z) 64 REAL*8 NORM(50,50) 65 REAL*4 FX(71),ABS,DIST DIMENSION X(100),NAME(20),C(50),W(100),FPHI(50) 66 DIMENSION EX(100) COMMON SING 67 EXTERNAL WEIGHT,PHI 68 NN=71 SING=1.0D-25 69 NMAX=50 71 MMAX=100 72 STOP = 999.D0 C C CALCULATION OF K C K=# OF YEARS OF DATA IN THE SERIES-->DEL (INCLUDING GAPS) C K=0 DO88 IJ=1,NN IF(FX(IJ).EQ.0.0)K=K+1.0 IF(FX(IJ).NE.0.0)GOTO 92 88 CONTINUE 92 K=NN-K DO89 JK=1,K X(JK)=JK 89 CONTINUE M=NN-K DO91 I=1,K FX(I)=FX(I+M) 91 CONTINUE C READ NUMBER OF BASE FUNCTIONS AND WEIGHT PARAMETER 79 C 80 C 82 C READ X,F(X),W(X) 83 C 84 C 91 C MAKE CHECKS 92 C 93 IF(N.LT.1)GO TO 150 94 IF(N.GT.NMAX)GO TO 160 95 IF(K.GT.MMAX)GO TO 170 96 IF(N.GT.K)GO TO 180 97 C 98 C IF WEIGHT NOT ASSIGNED COMPUTE IT 99 C 100 IF(IW.EQ.0)GO TO 35 101 DO 30 I=1,K 102 W(I)=WEIGHT(X(I)) 103 IF(FX(I).EQ.0.0)W(I)=0.0 30 CONTINUE C 104 C FORM THE NORMAL EQUATIONS 105 C 106 35 DO 50 I=1,N 107 FPHI(I)=0.0D0 108 DO 50 J=I,N 109 NORM(I,J)=0.0D0 110 DO 40 L=1,K 111 WW=W(L) 112 XX=X(L) 113 PHII=PHI(I,XX) 114 PHIJ=PHI(J,XX) 115 NORM(I,J)=NORM(I,J)+WW*PHII*PHIJ 116 40 IF(I.EQ.J)FPHI(I)=FPHI(I)+WW*PHII*FX(L) 117 50 NORM(J,I)=NORM(I,J) 118 C 119 C INVERT THE NORMALS - ERROR EXIT FOR DETERMINANT LESS THAN SING 120 C 121 CALL CHOLD(NORM,NMAX,N,DETA,&190) 122 C 123 C COMPUTE THE COEFFICIENTS C(I) 124 C 125 DO 60 I=1,N 126 C(I)=0.0D0 127 DO 60 J=1,N 128 60 C(I)=C(I)+NORM(J,I)*FPHI(J) 129 C 130 C COMPUTE LEAST SQUARES DISTANCE,VARIANCE AND STANDARD DEVIATION 131 C 132 RHO=0.0D0 133 DO 80 I=1,K 134 WW=W(I) 135 XX=X(I) 136 PN=0.0D0 137 DO 70 J=1,N 138 70 PN=PN+C(J)*PHI(J,XX) 139 80 RHO=RHO+WW*(FX(I)-PN)**2 140 MN=K-N 141 IF(MN.EQ.0)GO TO 90 142 VAR=RHO/DFLOAT(MN) 143 STD=DSQRT(VAR) 144 TVAR = NORM(2,2) * VAR TSTD = DSQRT (TVAR) TSTD = TSTD * 3048.0 TRD = C(2) * 3048.0 C C CALCULATION OF DELTA H (FEET) C DELTA=TRD/30.48 C C CALCULATION OF WEIGHT C Q=(1.D0/(TSTD**2.D0))*0.0338*DIST INTER=0 145 GO TO 100 146 90 VAR=0.0D0 147 STD=0.0D0 148 INTER=1 149 C 150 C PRINT ALL RESULTS 151 C 152 100 CONTINUE 153 PRINT 1040 154 DO 110 I=1,K 155 110 PRINT 1050,X(I),FX(I),W(I) 156 PRINT 1060 157 DO 120 I=1,N 158 120 PRINT 1070,I,C(I) 159 IF(INTER.EQ.0)PRINT 1080,RHO,VAR,STD 160 PRINT 7009 PRINT 7010,TRD,TSTD WRITE(6,7020) DELTA,Q 7020 FORMAT(2(/),T30,'DELTA H (FT) =',F8.3,T55,'Q=',F8.3,//) IF(INTER.EQ.1)PRINT 1090 161 RETURN 900 CONTINUE 140 STOP 172 150 PRINT 1130 173 STOP 174 160 PRINT 1140,N,NMAX 175 STOP 176 170 PRINT 1150,K,MMAX 177 STOP 178 180 PRINT 1160,K,N 179 STOP 180 190 PRINT 1170 181 STOP 182 1040 FORMAT(45X,'INPUT CHECK'//22X,'X',21X,'F(X)',23X,'W(X)'//) 187 1050 FORMAT('0', 5X,F20.5, 5X,F20.5, 5X,F20.5) 188 1060 FORMAT('-',40X,'COMPUTED COEFFICIENTS') 189 1070 FORMAT('0',30X,'C(',I2,')',10X,F20.10) 190 1080 FORMAT('-',10X,'RHO=',F20.10,10X,'VAR=',F20.10,10X,'STD=',F20.10) 191 1090 FORMAT('-',10X,'K=N : RHO=VAR=STD= 0.0D0') 192 1130 FORMAT('0',5X,'ERROR - DEGREE OF FUNCTION < 1 CHECK DATA CARD 2') 197 1140 FORMAT('0',5X,'ERROR - N=',I4,5X,'NMAX=',I4,5X,'CHECK DATA CARD 2 198 1OR INCREASE NMAX AND DIMENSION') 199 1150 FORMAT('0',5X,'ERROR - K=',I4,5X,'MMAX=',I4,5X,'CHECK XSTOP OR INC 200 1REASE MMAX AND DIMENSION') 201 1160 FORMAT('0',5X,'ERROR - K=',I4,5X,'N=',I4,5X,'CHECK INPUT OR INCREA 202 1SE NO. OF DATA POINTS OR DECREASE NO. OF BASE FUNCTIONS') 203 1170 FORMAT('0',5X,'ERROR - DETERMINANT OF GRAMM''S MATRIX .LT.SING, DE 204 1CREASE SING OR CHECK BASE FUNCTIONS OR WEIGHT FUNCTION') 205 7009 FORMAT (////,T30,'TREND (CM/CENT)',10X,'STD') 7010 FORMAT(T30,F15.10,T50,F15.10,//) END 206 FUNCTION WEIGHT(X) REAL*8 WEIGHT,X WEIGHT = 1.0D0 RETURN END FUNCTION PHI (I,X) REAL*8 PHI,X C GO TO BASE FUNCTION (I) GO TO (10,20), I 10 PHI = 1.0D0 RETURN 20 PHI = X RETURN END SUBROUTINE CHOLD(A,IRDA,NA,DETA,*) 220 DOUBLE PRECISION A,DETA,SUM,SQRT,DSQRT,ABS,DABS,SING 221 DIMENSION A(IRDA,NA) 222 COMMON SING 223 SQRT(SUM) = DSQRT(SUM) 224 ABS(DETA) = DABS(DETA) 225 DETA = A(1,1) 226 A(1,1) = SQRT(A(1,1)) 227 IF(NA .EQ. 1) GO TO 6 228 DO 1 I = 2,NA 229 1 A(I,1) = A(I,1) / A(1,1) 230 DO 5 J = 2,NA 231 SUM = 0. 232 J1 = J - 1 233 DO 2 K = 1,J1 234 2 SUM = SUM + A(J,K) ** 2 235 DETA=DETA*(A(J,J)-SUM) 236 A(J,J) = SQRT(A(J,J) - SUM) 237 IF(J .EQ. NA) GO TO 5 238 J2 = J + 1 239 DO 4 I = J2,NA 240 SUM = 0. 241 DO 3 K = 1,J1 242 3 SUM = SUM + A(I,K) * A(J,K) 243 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 244 5 CONTINUE 245 6 IF(ABS(DETA).LT.SING)GO TO 16 246 DO 7 I = 1,NA 247 7 A(I,I) = 1. / A(I,I) 248 IF(NA .EQ. 1) GO TO 10 249 N1 = NA - 1 250 DO 9 J = 1,N1 251 J2 = J + 1 252 DO 9 I = J2,NA 253 SUM = 0. 254 I1 = I - 1 255 DO 8 K = J,I1 256 8 SUM = SUM + A(I,K) * A(K,J) 257 9 A(I,J) = - A(I,I) * SUM 258 10 DO 15 J = 1,NA 259 IF(J .EQ. 1) GO TO 12 260 J1 = J - 1 261 DO 11 I = 1,J1 262 11 A(I,J) = A(J,I) 263 12 DO 14 I = J,NA 264 SUM = 0. 265 DO 13 K = I,NA 266 13 SUM = SUM + A(K,I) * A(K,J) 267 14 A(I,J) = SUM 268 15 CONTINUE 269 RETURN 270 16 RETURN 1 271 END 272