C C C C MEAN C C TO COMPUTE MEAN GRAVITY ANOMALIES FROM POINT DATA C VERSION 1.0 ... OCTOBER, 1973 ... C. MERRY. C C THE METHOD USES A SECOND-ORDER POLYNOMIAL TO FIT POINT DATA IN C EACH BLOCK. THIS IS INTEGRATED ANALYTICALLY TO OBTAIN A MEAN VALUE C FOR THE BLOCK. IF THERE ARE NOT SUFFICIENT DATA POINTS AVAILABLE, C A WEIGHTED ARITHMETIC MEAN IS USED INSTEAD. A FULL DESCRIPTION OF C THE METHOD IS IN :- MERRY AND VANICEK (1974) A METHOD FOR C ASTROGRAVIMETRIC GEOID DETERMINATION. TECH. REP. NO. 27, DEPT. OF C SURV. ENGG. , UNB, FREDERICTON, CANADA. C C C INPUT DATA C C UNIT 5 (CARD READER): C 1) INITIAL LAT. , LONG., STORAGE BLOCK SIZE OF POINT DATA. C (ALL IN DEGREES) - 2F8.3,F8.6 C 2) FINAL LAT. LONG., INVERSE OF MEAN BLOCK SIZE. C (ALL IN DEGREES) - 2F8.2,F8.4 C C UNIT 2 (TAPE) : C CONTAINS GRAVITY DATA SUPPLIED BY EARTH PHYSICS BRANCH, EMR. C C NOTES : C 1) FINAL AND INITIAL LAT., LONG. REFER TO LIMITS OF C RECTANGULAR AREA FOR WHICH MEANS ARE TO BE COMPUTED. FOR THE EMR C GRAVITY TAPE, LONG. IS POSITIVE WEST OF GREENWICH, AND STORAGE C BLOCK SIZE IS 1.0 DEGREES. C 2) THE EMR TAPE CONTAINS EXTRA DATA NOT NEEDED HERE. ONLY C THE FOLLOWING ARE NEEDED FOR EACH GRAVITY POINT: LAT., LONG., C ANOMALY (FREE-AIR OR BOUGUER ARE AVAILABLE), HEIGHT ACCURACY, C YEAR OF OBSERVATION. C C C MAIN PROGRAM SELECTS DATA FROM TAPE, IN STORAGE BLOCK SEGMENTS C AND TRANSFERS TO PORT IMPLICIT REAL*8(A-H,O-Z) DIMENSION FB(2000),AC(2000),P(2000),A(2000) COMMON /BLOCK1/ P,A,FB,AC,BPD,BLD DATA K /'-'/ DATA IZ /1H0/ 100 FORMAT(I2,I3,I6,F6.2,F7.2,I2,F7.1,I2,I5,I1,I1,F4.1,I2,F8.1,I2,F6.1 *,I2,F6.1,I2) 103 FORMAT('1',16X,'LAT.',8X,'LONG.',4X,'MEAN ANOM.',6X,'NO.',4X,'STD. * DEV.',//) 105 FORMAT('0',10X,F6.2,3X,F6.2,3X,F6.2) REWIND 2 C C READ AREA LIMITS AND BLOCKSIZE C READ(5,*) BP,BL,DIN READ(5,*) CX,CY,SIZE WRITE(6,105) BP,BL,DIN WRITE(6,105) CX,CY,SIZE WRITE(6,103) BLL=BL-DIN BPX=BP-DIN BPD=BP+DIN 4 CONTINUE BLD=BL+DIN C C READ DATA FROM TAPE C 2 BACKSPACE 2 J=1 10 READ(2,100,END=3)IY,IP,NUM,PDE,ALDE,IC,AHT,IHC,IEP,IDEP,IDC,TC, *ITC,GRAV,IGC,FA,IFC,BA,IBA PHI=PDE ALAM=ALDE C C TEST IF WITHIN AREA LIMITS C IF(PHI.LT.BPX) GOTO10 IF(ALAM.GT.CY) GOTO10 IF(ALAM.LT.BLL) GOTO10 C P(J)=PHI A(J)=ALAM FB(J)=FA C C CALCULATE ACCURACY OF ANOMALY, FROM HEIGHT ACCURACY C IF(IFC.EQ.0) ACC=12.0 IF(IFC.EQ.1) ACC=0.01 IF(IFC.EQ.2) ACC=0.05 IF(IFC.EQ.3) ACC=0.1 IF(IFC.EQ.5) ACC=1.0 IF(IFC.EQ.6) ACC=3.0 IF(IFC.EQ.7) ACC=5.0 IF(IFC.EQ.8) ACC=10.0 31 CONTINUE AC(J)=ACC IF(PHI.GE.BPD) GOTO1 IF(ALAM.GE.BLD) GOTO5 J=J+1 GOTO10 5 CONTINUE J=J-1 C C FOR EACH STORAGE BLOCK COMPLETELY READ, TRANSFER DATA TO PORT C CALL PORT(SIZE,J) BLD=BLD+DIN GOTO2 1 CONTINUE J=J-1 C C FOR EACH STORAGE BLOCK COMPLETELY READ, TRANSFER DATA TO PORT C CALL PORT(SIZE,J) IF(BPD.GT.CX) RETURN BPD=BPD+DIN GOTO4 3 CONTINUE STOP END SUBROUTINE PORT(SIZE,J) C C PORT SELECTS DATA FOR EACH MEAN BLOCK, TESTS FOR QUANTITY AND C DISTRIBUTION OF DATA IN BLOCK. IF INADEQUATE, WEIGHTED MEAN OF DATA C IS TAKEN. IF ADEQUATE, BEST-FITTING POLYNOMIAL IS COMPUTED AND C INTEGRATED, VIA ORTHO AND INTEG. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION D(9),SUMD(9,9) DIMENSION IQ(4) DIMENSION FB(2000),AC(2000),P(2000),A(2000) DIMENSION FZ(600),CA(600),AP(600),AL(600) COMMON /BLOCK1/ P,A,FB,AC,BPD,BLD COMMON /BLOCK2/ AP,AL COMMON /BLOCK3/ NK COMMON /BLOCK4/ D 101 FORMAT(3F12.3,I10,F12.3) 102 FORMAT(10X,3F12.3,I10,F12.3) IF(J.EQ.0) RETURN C C SET MEAN BLOCK SIZE AND CALCULATE NUMBER OF MEAN BLOCKS PER C STORAGE BLOCK C IDN=IFIX(SNGL(SIZE)) DN=1.D0/SIZE RHO=206264.80625D0/3600.D0 R=6371031.D0 C C SET ORDER OF POLYNOMIAL C NK=2 NAK=(NK+1)*(NK+1) BPS=BPD-1.D0 DO 11 M=1,IDN BPS=BPS+DN BPQ=BPS-DN BLS=BLD-1.D0 DO 10 KK=1,IDN BLS=BLS+DN BLQ=BLS-DN LI=0 C C SELECT DATA FOR EACH MEAN BLOCK C DO 9 L=1,J IF(P(L).LT.BPQ.OR.P(L).GT.BPS) GOTO4 IF(A(L).LT.BLQ.OR.A(L).GT.BLS) GOTO4 LI=LI+1 FZ(LI)=FB(L) CA(LI)=1.D0/AC(L)**2 AP(LI)=P(L) AL(LI)=A(L) IF(LI.GT.599) GOTO7 4 CONTINUE 9 CONTINUE 7 CONTINUE IF(LI.EQ.0) GOTO5 AI=LI+1.0D-40 PS=BPS-0.5D0*DN PL=BLS-0.5D0*DN DO 13 IE=1,4 13 IQ(IE)=0 C C TEST FOR DATA IN EACH QUADRANT AND FOR SUFFICIENCY (MORE THAN C 15 DATA POINTS) C DO 12 IE=1,LI AP(IE)=(AP(IE)-PS)/RHO*R AL(IE)=(AL(IE)-PL)/RHO*R*DCOS(PS/RHO) II=2.05+DSIGN(1.D0,AP(IE))+(1+DSIGN(1.D0,AL(IE)))/2 IQ(II)=1 12 CONTINUE ICH=IQ(1)*IQ(2)*IQ(3)*IQ(4) C C IF TEST FAILED, WEIGHTED ARITHMETIC MEAN IS COMPUTED IF(LI.GT.15.AND.ICH.NE.0) GOTO3 SUM=0.D0 ASP=0.D0 DO 6 K=1,LI SUM=SUM+FZ(K)*CA(K) ASP=ASP+CA(K) 6 CONTINUE AIJ=ASP+1.0D-50 AM=SUM/AIJ C C POLYNOMIAL TO COMPUTE ACCURACY OF MEAN ANOMALY C SNM=11.5575-0.317464*AI+3.6509D-03*AI**2-1.9788D-05*AI**3+4.9341D- *08*AI**4 SNM=0.3333333D0 IF(AI.GT.200) SNM=0.3D0 C C RESULTS WRITTEN ON PRINTER AND DISC C WRITE(6,102) PS,PL,AM,LI,SNM WRITE(11,101) PS,PL,AM,LI,SNM GOTO2 3 CONTINUE C C IF TEST PASSED , POLYNOMIAL IS FITTED TO DATA AND INTEGRATED C CALL ORTHO(LI,NAK,CA,FZ,9,D,SUMD) CALL INTEG(DN,PS,PL,SUMD,SMU,SK) SKK=DSQRT(SK) C C RESULTS WRITTEN ON PRINTER AND DISC C WRITE(6,102) PS,PL,SMU,LI,SKK WRITE(11,101) PS,PL,SMU,LI,SKK 2 CONTINUE 5 CONTINUE 10 CONTINUE 11 CONTINUE RETURN END SUBROUTINE INTEG(DN,PS,PL,SUMD,SMU,SK) C C INTEG COMPUTES THE MEAN GRAVITY ANOMALY USING THE ANALYTICAL C EVALUATION OF THE INTEGRATION OF APPROXIMATING POLYNOMIAL. C (SEE APPENDIX 1 OF TECH. REP. 27 ) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION SUMD(9,9) DIMENSION D(9),VD(9),VE(9) COMMON /BLOCK3/ NK COMMON /BLOCK4/ D RHO=206264.80625D0/3600.D0 R=6371031.D0 NAK=(NK+1)*(NK+1) DO 1 J=1,NAK 1 VD(J)=0.D0 DNN=0.5D0*DN A=DNN/RHO*R B=DNN/RHO*R*DCOS(PS/RHO) VD(1)=1.D0 VD(3)=B**2/3.D0 VD(7)=A**2/3.D0 VD(9)=VD(3)*VD(7) SMU=D(1)+D(3)*VD(3)+D(7)*VD(7)+D(9)*VD(9) DO 2 I=1,NAK SUM=0.D0 DO 3 J=1,NAK 3 SUM=SUM+SUMD(I,J)*VD(J) 2 VE(I)=SUM SK=0.D0 DO 4 I=1,NAK 4 SK=SK+VE(I)*VD(I) RETURN END SUBROUTINE ORTHO(N,M,W,F,MRD,D,SUMD) C THIS SUBROUTINE ORTHOGONALIZES THE MATRIX PHI USING THE GRAM-SCHMIDT 0275 C METHOD, COMPUTES THE FOURIER COEFFICIENTS OF THE ORTHOGONALIZED MATRIX 0276 C DERIVES THE COEFFICIENTS OF PHI,COMPUTES THE VARIANCES OF THE FOURIER 0277 C COEFFICIENTS AND THE VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 0278 C INPUTS : 0279 C 1. PHI(OPTIONAL - COULD BE FUNCTION SUBPROGRAM INSTEAD) - AN N BY M 0280 C CONTAINING THE BASE FUNCTIONS EVALUATED FOR EACH OBSERVATION 0281 C 2. N - THE NUMBER OF OBSERVATIONS 0282 C 3. M - THE NUMBER OF BASE FUNCTIONS (>2) 0283 C 4. W - WEIGHT FUNCTIONS C 5. F - FUNCTIONAL VALUES 0285 C 6. NRD - THE MAXIMUM ROW DIMENSION OF PHI C 7. MRD - THE MAXIMUM COLUMN DIMENSION OF PHI 0288 C OUTPUTS : 0295 C 1. ALPHA - AN MRD BY M MATRIX CONTAINING THE ALPHA'S USED IN COMPUTI 0296 C THE ORTHOGONALIZED MATRIX AND IN COMPUTING THE COEFFICIENTS OF PH 0297 C 2. C - THE M FOURIER COEFFICIENTS OF THE ORTHOGONALIZED MATRIX 0298 C 3. D - THE M COEFFICIENTS OF THE INPUT MATRIX PHI 0299 C 4. SUMD - THE VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 0301 C 5. SC2 - THE SQUARES OF THE NORMS OF THE ORTHOGONALIZED MATRIX C 6. V - THE N RESIDUALS 0304 IMPLICIT REAL*8(A-H,O-Z) 0305 DIMENSION W(N),F(N),D(M),SUMD(MRD,M) DIMENSION ALPHA(9,9),C(9),SC2(9) C TEST FOR NEGATIVE DEGREES OF FREEDOM 0310 IF (N.LT.M) GO TO 100 0311 K=1 0312 ALPHA(M,M)=1.D0 0313 C DETERMINE THE ALPHA'S FOR COMPUTATION OF ORTHOGONALIZED MATRIX 0314 10 DO 3 J=K,M 0315 IF(J.NE.K) GO TO 6 0316 ALPHA(K,K)=1.D0 0317 GO TO 3 0318 6 SC1=0.D0 0319 SC2(K)=0.D0 0320 SC3=0.D0 0321 DO 2 I=1,N 0322 P=PHI(I,K) 0323 IF(K.EQ.1) GO TO 4 0324 K1=K-1 0325 DO 5 J1=1,K1 0326 5 P=P+ALPHA(J1,K)*PHI(I,J1) 0327 4 SC1=SC1+W(I)*PHI(I,J)*P 0328 SC3=SC3+F(I)*W(I)*P 0329 2 SC2(K)=SC2(K)+W(I)*P**2 0330 ALPHA(J,K)=-SC1/SC2(K) 0331 ALPHA(K,J)=ALPHA(J,K) 0332 3 CONTINUE 0333 C DETERMINE THE FOURIER COEFFICIENTS FOR THE ORTHOGONALIZED MATRIX 0334 C(K)=SC3/SC2(K) 0335 K=K+1 0336 IF(K.LT.3) GO TO 10 0337 C DETERMINE THE ALPHA'S USED IN COMPUTING THE COEFFICIENTS OF PHI 0338 JK=K-1 0339 9 JL=K 0340 JK=JK-1 0341 JJ=K-JK-1 0342 DO 8 LM=1,JJ 0343 JL=JL-1 0344 8 ALPHA(JK,K)=ALPHA(JK,K)+ALPHA(JK,JL)*ALPHA(K,JL) 0345 IF(JK.NE.1) GO TO 9 0346 IF(K.LT.M) GO TO 10 0347 C DETERMINE THE LAST FOURIER COEFFICIENT 0348 SC2(K)=0.D0 0349 SC3=0.D0 0350 DO 7 I=1,N 0351 P=PHI(I,K) 0352 K1=K-1 0353 DO 1 J=1,K1 0354 1 P=P+ALPHA(J,K)*PHI(I,J) 0355 SC2(K)=SC2(K)+W(I)*P**2 0356 7 SC3=SC3+F(I)*W(I)*P 0357 C(K)=SC3/SC2(K) 0358 C DETERMINE THE COEFFICIENTS OF PHI 0359 DO 13 I=1,M 0363 D(I)=C(I) 0364 IF(I.EQ.M) GO TO 13 0365 II=I+1 0366 DO 14 J=II,M 0367 14 D(I)=D(I)+ALPHA(I,J)*C(J) 0368 13 CONTINUE 0369 C COMPUTE THE VARIANCE OF THE FOURIER COEFFICIENTS AND THE VARIANCE-COVA 0370 C MATRIX OF THE COEFFICIENTS 0371 DO 15 I=1,M 0372 DO 15 J=1,M 0373 15 SUMD(I,J)=0.D0 0374 SC4=0.D0 0375 DO 22 I=1,N 0376 PN=0.D0 0377 DO 21 J=1,M 0378 21 PN=PN+D(J)*PHI(I,J) 0379 V2=(F(I)-PN)**2 22 SC4=SC4+V2*W(I) 0382 SIGMA2=SC4/(N-M) DO 23 I=1,M 0389 SUMC=SIGMA2/SC2(I) DO 23 J=1,I 0390 DO 23 K=J,I 0391 23 SUMD(J,K)=SUMD(J,K)+ALPHA(J,I)*ALPHA(K,I)*SUMC DO 24 I=1,M 0393 IT=I+1 0394 IF(IT.GT.M) GO TO 30 0395 DO 24 J=IT,M 0396 24 SUMD(J,I)=SUMD(I,J) 0397 30 RETURN 100 WRITE(6,101) 0414 101 FORMAT('0','*ERROR* NEGATIVE DEGREES OF FREEDOM') 0415 RETURN END FUNCTION PHI(L,K) C C PHI COMPUTES THE ELEMENTS OF THE DESIGN MATRIX FOR ORTHO C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(600),Y(600) COMMON /BLOCK2/ X,Y COMMON /BLOCK3/ NK NN=NK+1 I=(K-1)/NN J=K-I-(NN-1)*I-1 IF(X(L).EQ.0.0.AND.I.EQ.0) X(L)=1.0 IF(Y(L).EQ.0.0.AND.J.EQ.0) Y(L)=1.0 PHI=X(L)**I*Y(L)**J RETURN END