DEFINE FILE 8(9458,130,E,IGRAV) C C ******************************************************************* C**** THIS MAIN PROGRAM READS GRAVITY DATA FROM DISK.THESE DATA ARE PO- C**** INT ANOMALIES STORED ON DISK ACCORDING TO THE NORMS OF THE IAGG. C**** THE POINT ANOMALIES ARE SEPARATED IN BLOCKS OF ONE DEGREE BY ONE C**** DEGREE. THE BLOCKS ARE IN CRESCENT ORDER OF LONGITUDE AND IN EACH C**** BLOCK THE DATA ARE IN CRESCENT ORDER OF LATITUDE. C**** AFTER ONE BLOCK IS RED, THE DATA ARE SEPARATED IN BLOCKS OF 20 MI C**** NUTES BY 20 MINUTES, A VALUE FOR THE MEAN ANOMALY IS EVALUATED FOR C**** EACH OF THESE BLOCKS AND FOR THE HOLE BLOCK OF ONE DEGREE BY ONE C**** DEGREE. C**** THE EVALUATION OF THE MEAN ANOMALY IS MADE BY BI-DIMENSIONAL POLI- C**** NOMIAL FITTING THECNIC WHEN THE DATA ARE ENOUGH OR BY A SIMPLE A- C**** RITHMETIC MEAN. C**** ALL THESE TASKS ARE ACCOMPLISHED USING THE SUBROUTINES THAT FOLLOW C**** THIS MAIN PROGRAM. C**** WRITTEN BY BLITZKOW - 1980 C ***************************************************************** C REAL *8 ANO(200),W(200),X(200),Y(200),XX(200),YY(200),ANOM(200, 11),WE(200),RM,RDL0,RDLA,AA,BB,ANM,VARAN,PCOE(121,1),SGM(121,121) 2,SNM,SNM1,SWE,SWE1 DIMENSION LAT(200),LON(200),LA(2),LO(2),LLO(200),LLA(200), 1NPP(50),XR(5),YR(5),XP(5),YP(5) C**** DEFINITION OF THE CONSTANTS AND INITIAL VALUES OF THE VARIABLES RM=6370900. NDP=2 IN=(NDP+1)*(NDP+1) IREC=4165 75 NT=0 ID1=0 I=1 J=1 C**** READING THE DATA 1 READ (8'IREC,10) LAT(I),LON(I),IERH,IERG,IFAA,NQUA 10 FORMAT (7X,I7,I8,10X,I1,16X,I1,37X,I5,28X,I6) IF (I.GT.1) GO TO 20 NQ=NQUA GO TO 21 20 IF (NQUA.NE.NQ) GO TO 40 21 IF (IERG-1) 16,17,17 16 SIGG=0.05 GO TO 22 17 SIGG=0.5 22 SIGH=0.5 C**** THE COLUMNS 27 AND 44 OF THE RECORD HAVE A CODE THAT REPRESENTS THE C**** GRAVITY VALUE AND ALTITUDE ERROR. THE CODE OF THE ALTITUDE IS NOT C**** USED AT PRESENT. AN ERROR OF 0.5M IS ASSUMED TO ALL POINTS. W(I)=1./(SIGG**2+(0.02969*0.5)**2) ANO(I)=FLOAT(IFAA)/100. I=I+1 IREC=IREC+1 GO TO 1 C**** SEPARATE THE DATA IN BLOCKS OF 20 MINUTES BY 20 MINUTES C**** ACCORDING TO THE LONGITUDE FOR EACH BLOCK OF 1 DEGREE BY 1 DEGREE. 40 IS=I-1 IREC=IREC-1 IF (IS.LT.3) GO TO 45 IF (IS.LT.15) ID1=1 CALL BLO20 (LAT,LON,ANO,W,IS,LLA,LLO,ANOM,WE,NPP) WRITE (6,2)((I,NPP(I)),I=1,9) 2 FORMAT(5X,5HBLOCK,I4,3X,17HNUMBER OF POINTS=,I3,/) C**** START THE COMPUTATION TO FIT THE POLYNOMIAL FOR BLOCKS OF 20 C**** MINUTES BY 20 MINUTES C**** DEFINITION OF THE INITIAL VALUES SWE=0 SNM=0 J=1 DC=20./60. DC1=DC DC2=-DC CX=(10.*4.*ATAN(1.)*RM)/(60.*180.) WRITE (6,3) 3 FORMAT (2(6X,8HLATITUDE,3X,9HLONGITUDE,3X,7HANOMALY),/) WRITE (6,4) (LLA(KK),LLO(KK),ANOM(KK,1),KK=1,IS) 4 FORMAT (2(7X,I7,4X,I8,3X,F7.2)) DO 60 I=1,9 SWE1=0 SNM1=0 ID=0 LL=NPP(I) IF(LL.LT.16.AND.LL.GT.2) ID=1 IF (LL.LT.3) ID=2 IF (LL.EQ.0) GO TO 60 DO 71 JJ=1,LL FLA=FLOAT(LLA(J))/10000. FLO=FLOAT(LLO(J))/10000. IF(J.GT.1) GO TO 51 ILA=FLA ILO=FLO FLO0=FLOAT(ILO)-0.5 IF(FLA) 52,53,53 52 FLA0=FLOAT(ILA)-0.5 GO TO 51 53 FLA0=FLOAT(ILA)+0.5 51 RDLA=4.*ATAN(1.)*FLA/180. IF (ID1.NE.1) GO TO 41 SWE=SWE+1./WE(J) SNM=SNM+ANOM(J,1)/WE(J) 41 X(J)=(FLA-FLA0)*4.*ATAN(1.)*RM/180000. Y(J)=(FLO-FLO0)*4.*ATAN(1.)*DCOS(RDLA)*RM/180000. GO TO (54,11),ID GO TO (61,62,63,64,65,66,67,68,69),I 62 DC2=0 GO TO 61 63 DC2=DC GO TO 61 64 DC1=0 DC2=-DC GO TO 61 65 DC1=0 DC2=0 GO TO 61 66 DC1=0 DC2=DC GO TO 61 67 DC1=-DC DC2=-DC GO TO 61 68 DC1=-DC DC2=0 GO TO 61 69 DC1=-DC DC2=DC 61 FA0=FLA0+DC1 FO0=FLO0+DC2 RDL0=FA0*4.*ATAN(1.)/180. CY=(10.*4.*ATAN(1.)*RM*DCOS(RDL0))/(60.*180.) XX(JJ)=(FLA-FA0)*4.*ATAN(1.)*RM/180. YY(JJ)=(FLO-FO0)*4.*ATAN(1.)*RM*DCOS(RDLA)/180. GO TO 55 54 SWE1=SWE1+1./WE(J) SNM1=SNM1+ANOM(J,1)/WE(J) 55 J=J+1 71 CONTINUE NTP=JJ-1 IF (ID.EQ.1) GO TO 56 CALL POLFIT (XX,YY,ANOM,NTP,WE,NDP,PCOE,SGM) CALL GMA (PCOE,CX,CY,NDP,SGM,ANM,VARAN) GO TO 57 56 ANM=SNM1/SWE1 VARAN=SWE1/FLOAT(LL-1) GO TO 73 57 WRITE (6,5) 5 FORMAT (/,10X,30HCOEFFICIENTS OF THE POLYNOMIAL) CALL MOUTD (PCOE,121,IN,1) 73 WRITE (6,6) ANM,VARAN,NTP,NQ,I 6 FORMAT (6X,13HMEAN ANOMALY=,F7.2,6X,9HVARIANCE=,F8.5,6X,17HNUMBER 1OF POINTS=,I3,6X,12HGRID NUMBER=,I6,1H-,I2,//) GO TO 60 11 J=J+NPP(I) 60 CONTINUE RDLA0=FLA0*4.*ATAN(1.)/180. AA=(0.5*4.*ATAN(1.)*RM)/180000. BB=(0.5*4.*ATAN(1.)*RM*COS(RDLA0))/180000. IF (ID1.EQ.1) GO TO 58 CALL POLFIT (X,Y,ANOM,IS,WE,NDP,PCOE,SGM) CALL GMA (PCOE,AA,BB,NDP,SGM,ANM,VARAN) GO TO 59 58 ANM=SNM/SWE VARAN=SWE/FLOAT(IS-1) GO TO 72 59 WRITE (6,7) 7 FORMAT (//,10X,30HCOEFFICIENTS OF THE POLYNOMIAL) CALL MOUTD (PCOE,121,IN,1) 72 WRITE (6,8) ANM,VARAN,IS,NQ 8 FORMAT (6X,13HMEAN ANOMALY=,F7.2,6X,9HVARIANCE=,F8.5,6X,17HNUMBER 1OF POINTS=,I3,6X,12HGRID NUMBER=,I6) XAREA=30.*COS(RDLA0) YAREA=30. RDLA1=FLOAT(ILA)*4.*ATAN(1.)/180. BB=AA CALL AREA (XAREA/2.54,YAREA/2.54) CALL SETPLT (-BB,-AA,BB,AA) CALL PRNTCH ('*') XR(1)=-BB XR(2)=BB XR(3)=BB XR(4)=-BB XR(5)=XR(1) YR(1)=AA YR(2)=AA YR(3)=-AA YR(4)=-AA YR(5)=AA CALL NOWPLT (0,XR(1),YR(1)) CY=CX DO 80 I=2,5 80 CALL NOWPLT (1,XR(I),YR(I)) DO 81 I=1,4 GO TO (82,83,84,85),I 82 XP(1)=BB-4.*CX/1000. YP(1)=AA XP(2)=XP(1) YP(2)=-AA GO TO 86 83 XP(1)=BB-2.*CX/1000. YP(1)=AA XP(2)=XP(1) YP(2)=-AA GO TO 86 84 XP(1)=-BB YP(1)=AA-2.*CY/1000. XP(2)=BB YP(2)=YP(1) GO TO 86 85 XP(1)=-BB YP(1)=AA-4.*CY/1000. XP(2)=BB YP(2)=YP(1) 86 CALL NOWPLT (0,XP(1),YP(1)) 81 CALL NOWPLT (1,XP(2),YP(2)) CALL PRNTCH ('+') DO 87 I=1,IS 87 CALL NOWPLT (0,Y(I),X(I)) CALL ENDPLT GO TO 88 45 WRITE (6,9) IS,NQ 9 FORMAT (/,10X,26HTHE BLOCK HAS A FEW POINTS,/,10X,17HNUMBER OF POI 1NTS=,I2,10X,12HGRID NUMBER=,I6) 88 IREC=IREC+1 IF (IREC.EQ.4233) GO TO 90 GO TO 75 90 STOP END SUBROUTINE BLO20 (LAT,LON,ANO,W,NP,LLA,LLO,AANO,WW,NNP) C**** PROGRAM TO SEPARATE THE DATA OF A BLOCK OF ONE DEGREE BY ONE DEGR C**** EE IN BLOCKS OF 20 MINUTES BY 20 MINUTES. C**** THE ORIGINAL DATA ARE SUPPOSED TO BE IN CRESCENT ORDER OF LATITU C**** DE. THE FINAL DATA WILL BE IN BLOCKS OF CRESCENT ORDER OF LONGITU C**** DE. C**** INPUT: C**** LAT,LON - GEODESIC COORDINATES OF THE POINTS. C**** ANO - GRAVITY ANOMALY. C**** W - VARIANCE OF THE ANOMALY. C**** NP - NUMBER OF POINTS IN THE BLOCK. C**** OUTPUT: C**** LLA,LLO -GEODESIC COORDINATES OF THE POINTS SEPARATED. C**** AANO - GRAVITY ANOMALY. C**** WW - VARIANCE OF THE ANOMALY. C**** NNP - ARRAY WITH DIMENSION 9 WITH THE NUMBER OF POINTS OF C**** IN EACH 20 MINUTES BY 20 MINUTES BLOCK. REAL *8 ANO(200),AANO(200,1),W(200),WW(200) DIMENSION LAT(200),LON(200),NNP(9),LA(2),LO(2),NP1(3),LLO(200), 1LLA(200) NT=0 SNP1=0 J=1 DO 10 I=1,NP IF (I.GT.1) GO TO 20 IK=1 FLA= FLOAT(LAT(I))/10000. FLO= FLOAT(LON(I))/10000. ILA=FLA ILO=FLO IF(FLA)11,12,12 11 RLA= (FLOAT(ILA)-0.3333)*10000. LA(1)=RLA LA(2)=LA(1)-3333 GO TO 13 12 RLA=(FLOAT(ILA)+0.3333)*10000. LA(2)=RLA LA(1)=LA(2)+3333 13 RLO=(FLOAT(ILO)-0.3333)*10000. LO(2)=RLO LO(1)=LO(2)-3333 C*** SEPARATE THE DATA BY STRIPS OF 20 MINUTES BY 20 MINUTES C*** ACCORDING TO LATITUDE 20 IF (IK.EQ.3) GO TO 10 IF (LAT(I).LT.LA(J)) GO TO 30 GO TO 10 30 IF (I-1) 14,14,15 14 NP1(IK)=0 J=J+1 IK=IK+1 GO TO 20 15 NP1(IK)=I-1-NT NT=I-1 IK=IK+1 IF(IK.EQ.3) GO TO 10 J=J+1 GO TO 20 10 CONTINUE IF (IK.LT.3) NP1(IK+1)=0 IF (IK.LT.2) NP1(IK+2)=0 NP1(IK)=I-1-NT NT=0 IJ=1 J=1 I1=1 C C**** SEPARATE THE DATA IN BLOCKS OF 20 MINUTES BY 20 MINUTES C**** ACCORDING TO LONGITUDE FOR EACH BLOCK OF 1 DEGREE BY 1 DEGREE. C DO 50 N=1,3 ID=1 IF (NP1(N).NE.0) GO TO 49 DO 42 JJ=1,3 IJ=JJ+3*(N-1) 42 NNP(IJ)=0 IJ=IJ+1 GO TO 50 49 NN=NP1(N) DO 48 II=1,NN IL=II+SNP1 C**** ID=1 THERE IS NO LATITUDE AND LONGITUDE IQUALS TO ZERO. IF (ID.EQ.1) GO TO 47 IF (LON(IL).EQ.0) GO TO 48 C**** 'J' IS THE NUMBER OF TIMES THAT THE STRIP IS PROCESSED. IN C**** THE LAST TIME IS NOT NECESSARY TO TEST FOR LO(J). IF (J.EQ.3) GO TO 52 47 IF (LON(IL)-LO(J)) 52,52,48 52 LLO(I1)=LON(IL) LLA(I1)=LAT(IL) AANO(I1,1)=ANO(IL) WW(I1)=W(IL) LON(IL )=0 LAT(IL)=0 ANO(IL)=0 W(IL)=0 I1=I1+1 48 CONTINUE NNP(IJ)=I1-1-NT NT=I1-1 IJ=IJ+1 ID=ID+1 J=J+1 IF(J.LT.4) GO TO 49 SNP1=SNP1+NP1(N) J=1 50 CONTINUE RETURN END SUBROUTINE POLFIT (XX,YY,OBV,NO,W,NP,CO,SIGM) C C FIT A BI-DIMENSIONAL POLYNOMIAL AND COMPUTE THE VARIANCE-COVARIAN- C CE MATRIX OF THE ADJUSTED COEFICIENTS. C INPUT: C XX,YY-PLANE COORDINATES OF THE POINTS IN METERS. C OBV-OBSERVED VALUES,E.G.,GRAVITY ANOMALY,GEOIDAL HEIGHT,ETC. C NO- NUMBER OF OBSERVATIONS. THE PROGRAM IS DIMENSIONED TO A C MAXIMUM OF 500 OBSERVATIONS. C W-VECTOR OF THE INVERSE VARIANCES OF THE OBSERVED VALUES C NP-ORDER OF THE POLYNOMIAL IN THE MERRY'S SENSE.THE PROGRAM C IS DIMENSIONED TO A MAXIMUM ORDER OF TEN. C OUTPUT: C CO-VECTOR OF THE ADJUSTED COEFFICIENTS OF THE POLYNOMIAL C SIGM- VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED COEFFICIENTS C C REAL *8 XX(200),YY(200),OBV(200,1),W(200),CO(121,1),SIGM(121, 1121),WE(200,200),A(200,121),EINV(121,121),C,OBP,V,VTWV,SIG2 2,TERM VTWV=0 NT=(NP+1)*(NP+1) DO 10 K=1,NO DO 10 J=1,NT 10 A(K,J)=TERM(XX(K),YY(K),NP,J) DO 22 J=1,NO DO 22 K=1,NO IF(J-K) 11,12,11 11 WE(J,K)=0. GO TO 22 12 WE(J,K)=W(K) 22 CONTINUE C C COMPUTE AND SOLVE THE NORMAL EQUATIONS C CALL NOREQ (A,NO,NT,WE,OBV,CO,EINV) C C COMPUTATION OF THE GRAVITY ANOMALY OF EACH POINT USING THE ADJUSTED C POLYNOMIAL AND COMPUTATION OF THE RESIDUALS C DO 20 K=1,NO OBP=0. DO 30 J=1,NT C=TERM(XX(K),YY(K),NP,J) 30 OBP=OBP+CO(J,1)*C V=OBV(K,1)-OBP 20 VTWV=VTWV+V*V*W(K) DGF=FLOAT(NO-NT) SIG2=VTWV/DGF WRITE (6,6) SIG2 6 FORMAT (/,10X,29HVARIANCE FACTOR A POSTERIORI=,F10.4,) C C COMPUTATION OF THE VARIANCE-COVARIANCE MATRIX C DO 40 K=1,NT DO 40 J=1,NT 40 SIGM(K,J)=SIG2*EINV(K,J) RETURN END 0393 ¶bbbBBB ý- 0 1ý ¶BBBBJ0 ËËËËËËËË FUNCTION TERM(X,Y,N,L) 1 C COMPUTE EACH TERM X**J*Y**K OF THE BIDIMENSIONAL POLYNOMIAL 1 C INPUT: 1 C X,Y-COORDINATES OF THE POINT 1 C N-ORDER OF THE POLYNOMIAL IN THE MERRY'S SENSE 1 C L-ORDER OF THE TERM IN THE POLYNOMIAL 1 C IN A POLYNOMIAL OF THE SECOND ORDER, FOR EXAMPLE,THE ORDER OF THE 1 C COEFFICIENTS IS:C00,C01,C02,C10,C11,C12,C20,C21,C22 1 C L=1 FOR C00, L=2 FOR C01, L=3 FOR C02, AND SO ON 1 C 1 REAL *8 X,Y,TERM 1 IN=N+1 1 J=(L-1)/IN 1 K=L-J*IN-1 1 IF(J.EQ.0.AND.K.EQ.0) GO TO 20 1 IF(J.EQ.0) GO TO 21 1 IF(K.EQ.0) GO TO 22 1 TERM= X**J*Y**K 1 GO TO 10 1 20 TERM=1. 1 GO TO 10 1 21 TERM=Y**K 1 GO TO 10 1 22 TERM=X**J 1 10 RETURN 1 END 1 SUBROUTINE NOREQ (A,NR,NC,W,O,X,NN) 1 C COMPUTE THE NORMAL EQUATIONS FROM THE MATHEMATICAL MODEL IN THE 1 C PARAMETRIC FORM 1 C INPUT: 1 C A-IS THE MATRIX OF THE COEFFICIENTS IN THE LINEAR MODEL 1 C NR,NC-NUMBER OF ROWS AND COLUMNS RESPECTIVELY OF THE MATRIX A 1 C W-IS THE WEIGHT MATRIX 1 C O-IS THE VECTOR OF OBSERVATIONS 1 C OUTPUT: 1 C X-IS THE OUTPUT VECTOR OF COMPUTED UNKNOWNS 1 C N-IS THE INVERSE MATRIX OF THE COEFFICIENT OF THE NORMAL EQUATIONS. 1 C THIS MATRIX WILL BE NECESSARY TO COMPUTE THE VARIANCE-COVARIANCE 1 C MATRIX SIGM=SIG2*N**(-1) 1 REAL *8 NN(121,121),D 1 REAL *8 A(200,121),W(200,200),O(200,1),X(121,1),AT(121,200),B( 1 1121,200),U(121,1),AN(121,121) 1 CALL TRNSD (AT,121,A,200,NR,NC) 1 CALL MMULD (B,121,AT,121,W,200,NC,NR,NR) 1 CALL MMULD (AN,121,B,121,A,200,NC,NR,NC) 1 DO 10 I=1,NC 1 DO 10 J=1,NC 1 10 NN(I,J)=AN(I,J) 1 CALL MMULD (U,121,B,121,O,200,NC,NR,1) 1 CALL SPIN (NN,NC,121,D,IEX) 1 CALL MMULD (X,121,NN,121,U,121,NC,NC,1) 1 RETURN 1 END 1 SUBROUTINE GMA (C,A,B,NP,SIGM,GMAN,SIGAN) C C COMPUTE THE MEAN ANOMALY OF A RECTANGLE OF DIMENSIONS 2*A BY 2*B C AFTER A BI-DIMENSIONAL POLYNOMIAL OF ORDER NP HAS BEEN INTERPOLA- C TED TO GIVEN POINT ANOLALIES. THE MEAN ANOMALY IS THE INTEGRAL OF C THE POLYNOMIAL DIVIDED BY THE AREA OF THE RECTANGLE. C INPUT: C C- POLYNOMIAL COEFFICIENTS COMPUTED PREVIOUSLY. C A,B- LINEAR DIMENSIONS OF THE SEMI-SIDES OF THE RECTANGLE. C NP- ORDER OF THE POLYNOMIAL IN THE MERRY'S SENSE. THE PROGRAM C IS DIMENSIONED TO A MAXIMUM ORDER OF TEN. C SIGM- THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED COEFFI- C CIENTS OF THE POLYNOMIAL. C OUTPUT: C GMAN- GRAVITY MEAN ANOMALY OF THE RECTANGLE. C SIGAN- VARIANCE OF THE COMPUTED MEAN ANOMALY. C REAL *8 C(121,1),SIGM(121,121),IN(100),CC(100,1),CCT(1,100), 1SIGMA(100,100),GMAN,SIGAN,CT(1,100) C C DEFINITION OF THE CONSTANTS USED IN THE SUBROUTINE C SIGAN=0. GMAN=0. IJ=NP/2*2+2 C VERIFY IF THE ORDER OF THE POLYNOMIAL IS ODD OR EVEN TPAR= FLOAT(NP)/2.-FLOAT(NP/2) IP=NP+3 IF(TPAR.EQ.0) IP=NP+2 ID=0 IND=1 C C COMPUTATION OF THE MEAN ANOMALY C IF(NP.EQ.1) GO TO 10 DO 20 I=2,IJ,2 IF(I-2) 15,15,16 16 IND=IND-2+IP 15 II=I-2 DO 20 K=2,IJ,2 ID=ID+1 IN(ID)=IND KK=K-2 FR=FLOAT((II+1)*(KK+1)) IF(IND.EQ.1) GO TO 21 COE=C(IND,1)/FR IF(II) 11,11,12 11 GMAN=GMAN+COE*B**KK CC(ID,1)=(B**KK)/FR GO TO 23 12 IF(KK) 13,13,14 13 GMAN=GMAN+COE*A**II CC(ID,1)=(A**II)/FR GO TO 23 14 GMAN=GMAN+COE*A**II*B**KK CC(ID,1)=(A**II*B**KK)/FR GO TO 23 21 GMAN=GMAN+C(IND,1) CC(IND,1)=1. 23 IND=IND+2 20 CONTINUE GO TO 30 10 GMAN=C(IND,1) SIGAN=SIGM(IND,IND) GO TO 24 C C COMPUTATION OF THE VARIANCE OF THE MEAN ANOMALY C 30 DO 40 I=1,ID I1=IN(I) DO 40 J=1,ID I2=IN(J) 40 SIGMA(I,J)=SIGM(I1,I2) CALL TRNSD (CCT,1,CC,100,ID,1) CALL MMULD (CT,1,CCT,1,SIGMA,100,1,ID,ID) DO 50 I=1,ID 50 SIGAN= SIGAN+CT(1,I)*CC(I,1) 24 RETURN END