C TEST OF COVARIANCE FUNCTION SUBROUTINES, PROGRAMMED BY C.C.TSCHERNING, C DEP.GEODETIC SCIENCE, OSU AND GEODAETISK INSTITUT,KOEBENHAVN JUNE 75. C REFERENCES: C (A) TSCHERNING,C.C.: COVARIANCE EXPRESSIONS FOR SECOND AND LOWER ORDER C DERIVATIVES OF THE ANOMALOUS POTENTIAL, REPORTS OF THE DEP. OF C GEODETIC SCIENCE NO. 225,1975. C IMPLICIT REAL *8(A-H,O-Z),LOGICAL (L) COMMON /CMCOV/CI(12),CR(51),SIGMA0(300),SIGMA(300),HMAX,KI(25),N1, *N2,LOCAL,LSUM C THE COMMON AREA IS USED FOR THE TRANSFER OF DATA TO AND FROM THE SUB- C ROUTINE COVAX. DIMENSION COV(181,7),KP(7),KQ(7),HP(7),HQ(7) *,KX(3),KY(3),IP(3),IQ(3),LA(3),SM(2001) DATA GM,RE/3.98D14,6371.0D3/ *,D0,D1,D2/0.0D0,1.0D0,2.0D0/,PI/3.1415926535D0/ C RE IS THE MEAN RADIUS OF THE EARTH AND GM IS THE PRODUCT OF THE GRAVI- C TATIONAL CONSTANT AND THE MASS OF THE EARTH. C WRITE(6,10) 10 FORMAT('0TEST OF COVARIANCE FUNCTION SUBROUTINES, VERS. JUNE 75.', *//,' COVARIANCES BETWEEN QUANTITIES OF KIND KP,KQ ARE COMPUTED. ', */,' THE KINDS AND CORRESPONDING UNITS ARE AS FOLLOWS: (E=EOTVOS' *,'):',/, *' (1) THE HEIGHT ANOMALY (METERS), (2) THE NEGATIVE RADIAL DER-',/ *,' IVATIVE DIVIDED BY THE RADIAL DISTANCE (E), (3) THE GRAVITY',/, *' ANOMALY (MGAL), (4) THE RADIAL DERIVATIVE OF (3) (E), (5) THE',/ *' SECOND ORDER RADIAL DERIVATIVE (E), (6),(7) THE LATITUDE AND',/, *' THE LONGITUDE COMPONENTS OF THE DEFLECTIONS OF THE VERTICAL',/, *' (ARCSECONDS), (8),(9) THE DERIVATIVES OF (3) IN NORTHERN AND',/, *' EASTERN DIRECTION, RESPECTIVELY (E), (10),(11) THE DERIVATIVE',/ *,' OF (2) IN THE SAME DIRECTIONS (E), (12)-(14) THE SECOND ORDER', */,' DERIVATIVES IN NORTHERN, (NORTHERN,EASTERN) AND EASTERN',/, *' DIRECTIONS, RESPECTIVELY (E).',//) C C INPUT OF THE VALUE OF A LOGICAL VARIABLE, LTEST, TRUE WHEN TEST OUT- C PUT IS NEEDED. READ(5,11)LTEST 11 FORMAT(L2) C INPUT OF QUANTITIES SPECIFYING THE DEGREE-VARIANCE MODEL TO BE USED IN C THE FOLLOWING SEQUENCE: THE RATIO BETWEEN AN ADOPTED BJERHAMMAR-SPHERE C RADIUS (RB) AND RE, SQUARED, THE QUANTITY A(I) IN REF(A),EQ.(17), DI- C VIDED BY RB**2 IN UNITS OF MGAL**2, THE INTEGERS K(2),K(3) OF EQ. C (17), WHEN APPLICABLE, OTHERWISE A ZERO, THE VALUE OF A LOGICAL VARI- C ABLE, LOCAL, TRUE,WHEN THE DEGREE-VARIANCES UP TO-AND INCLUSIVE DEGREE C N ARE ZERO AND FALSE WHEN EMPIRICAL ANOMALY DEGREE-VARIANCES UP TO C ORDER N WILL BE INPUT, THE INTEGER N, THE INTEGER KT EQUAL TO THE C DEGREE-VARIANCE MODEL NUMBER (1,2 OR 3), THE VALUE OF THE LOGICAL C VARIABLE LSUM, WHICH IS TRUE WHEN A FINITE LEGENDRE SERIES (MAXIMAL C DEGREE 2000) MUST BE USED FOR THE EVALUATION OF COVARIANCES IN ALTI- C TUDES GREATHER THAN HMAX AND OTHERWISE EQUAL TO FALSE, THE VALUE OF C N2 (WHICH MUST BE LESS THAN OR EQUAL TO 2000), THE VALUE OF THE HEIGHT C HMAX IN METERS, AND FINALLY THE VALUE OF LAST1, TRUE, WHEN THE LIST C OF INPUT PARAMETERS IS THE LAST ONE. C 100 READ(5,9)S,A,KI(3),KI(4),LOCAL,N,KT,LSUM,N2,HMAX,LAST1 9 FORMAT(2D14.7,2I5,L2,2I5,L2,I5,D14.7,L2) IF (N2.LT.2) N2 = 2 IF (N2.GT.2001) N2 = 2000 N2 = N2+1 KI(5) = KT WRITE(6,12)S,A,KI(3),KI(4),N,KT 12 FORMAT('0PARAMETERS SPECIFYING THE MODEL DEGREE-VARIANCES:',/, *' S,A =',2D14.7,/,' K0-K3,N,KT= -2 -1',4I5/) IF (LSUM) WRITE(6,6)HMAX,N2 6 FORMAT(' WHEN THE HEIGHT OF ONE OF THE POINTS OF EVALUATION IS ABO *VE ',D14.7,' METERS',/,' WILL THE COVARIANCES BE EVALUATED BY MEAN *S OF A LEGENDRE SERIES HAVING ',I5,' TERMS.') RB2 = RE*RE*S LMODEL = N.EQ.0 C CI(8) = A*RB2*1.0D-10 CI(10) = S C CONVERTING THE CONSTANT A INTO UNITS OF (M/SEC)**4 AND SUBSEQUENT C STORING OF A AND S IN THE ARRAY CI ACCORDING TO THE SPECIFICATIONS C GIVEN IN COVAX. C IF (N.NE.0) GO TO 101 C N EQUAL TO ZERO IMPLIES, THAT ALL DEGREE-VARIANCES ARE EQUAL TO THE C MODEL DEGREE-VARIANCES. THIS AGAIN IS EQUIVALENT TO HAVING THE DEGREE- C VARIANCES OF DEGREE 0,1,2 EQUAL TO ZERO, I.E. THE COVARIANCE FUNCTION C USED IS A 2'-ORDER LOCAL COVARIANCE FUNCTION. LOCAL = .TRUE. N = 2 C 101 N1 = N+1 C INPUT OF EMPIRICAL ANOMALY DEGREE-VARIANCES IN UNITS OF MGAL**2. IF(.NOT.LOCAL)READ(5,13)(SIGMA0(I), I = 1, N1) 13 FORMAT(12F6.2) IF (.NOT.LOCAL)WRITE(6,7)(SIGMA0(I), I = 1, N1) 7 FORMAT(' EMPIRICAL ANOMALY DEGREE-VARIANCES IN UNITS OF MGAL**2:', */,25(12F6.2/)) C N2 = 2001 CALL COVAX(SM) C THE ARRAY SM OCCURRING IN THE CALL OF COVAX IS USED TO STORE THE C DEGREE-VARIANCES WHEN LSUM IS TRUE. THE DIMENSION OF THE ARRAY IS C TRANSFERRED TO THE SUBROUTINE BY MEANS OF THE VARIABLE N2 OCCURRING C IN THE COMMON AREA /CMCOV/. C C CONVERSION OF ANOMALY DEGREE-VARIANCES TO POTENTIAL-DEGREE-VARIANCES. WRITE(6,8)(SIGMA0(I), I = 1, N1) 8 FORMAT('0MODEL DEGREE-VARIANCE CORRECTIONS:',/,50(6(1X,D11.4),/)) C KPP=0 KQQ=0 C INPUT OF QUANTITIES SPECIFYING THE OUTPUT TABLE. THE TABLE WILL CON- C TAIN MT COLUMNS OF COVARIANCES OF KINDS KP, KQ (TO BE INPUT SUBSEQUE- C NTLY) COMPUTED IN POINTS P AND Q. EACH COLUMN WILL CONTAIN NCOV VALUES C CORRESPONDING TO Q MOVING IN AN AZIMUTH (ALFA) IN STEPS OF LENGTH DT C (MINUTES). THE SPECIFICATION CONSIST OF NCOV, DT, ALFA IN DEG.,MIN., C SEC AND A LOGICAL, LAST2, WHICH IS TRUE WHEN THIS IS THE LAST SPECI- C FICATION FOR THE DEGREE-VARIANCE MODEL UNDER CONSIDERATION. 102 READ(5,14)NCOV,DT,IDEG,MIN,SEC,LAST2 14 FORMAT(I5,F7.2,I5,I3,F6.2,L2) CALL RAD(ALFA,IDEG,MIN,SEC,1) RT = 60.0*DT/206264.806 LPOLE = DABS(ALFA).LT.1.0D-6.OR.DABS(ALFA-PI).LT.1.0D-6 IF (LPOLE) GO TO 103 C CA = DCOS(ALFA) SA = DSIN(ALFA) 103 IF (LTEST) WRITE(6,15)NCOV,DT,IDEG,MIN,SEC,ALFA 15 FORMAT('0NCOV,DT,DEG,MIN,SEC,ALFA=',I4,F6.2,I5,I3,F6.2,D14.6) C MV = 0 MT = 1 C IF (NCOV.LE.181) GO TO 104 NCOV = 181 WRITE(6,37) 37 FORMAT(' NCOV TOO BIG, FIXED TO 181.') 104 MT = MT+1 IF (MV.NE.0.AND.MV.NE.3) GO TO 112 C INPUT OF INTEGERS KP AND KQ SIGNIFYING THE KIND OF QUANTITIES BETWEEN C WHICH WE WANT TO COMPUTE THE COVARIANCES. THE VALUES OF KP,KQ MUST BE C EQUAL TO THE EQUATION NUMBERS OF REF(A), WHICH DEFINES THE QUANTITIES C (1) - (14). C ON THE SAME PUNCH CARD INPUT OF THE HEIGHTS OF P AND Q AS WELL,(IN C METERS),AND OF A LOGICAL VARIABLE LAST3, TRUE WHEN THIS KIND OF COVA- C RIANCES ARE THE LAST ONES TO BE COMPUTED WITH THE CHOOSED FORM OF THE C TABLE. THREE SETS OF VALUES MAY BE PUNCHED ON ONE CARD, CF. FORMAT C STATEMENT 15. READ(5,16)(KX(K),KY(K),IP(K),IQ(K),LA(K),K = 1, 3) 16 FORMAT(3(2I3,2I8,L2)) MV = 0 112 MV = MV+1 KP(MT) = KX(MV) KQ(MT) = KY(MV) HP(MT) = IP(MV) HQ(MT) = IQ(MV) LAST3 = LA(MV) C KI(6) = KP(MT) KI(7) = KQ(MT) LNEW = KP(MT).NE.KPP.OR.KQ(MT).NE.KQQ C COMPUTATION OF CONSTANTS NEEDED FOR THE COVARIANCE COMPUTATION, WHICH C ARE INDEPENDENT OF T AND THE HEIGHTS BY THE CALL OF COVBX. IF (LNEW) CALL COVBX KPP = KP(MT) KQQ = KQ(MT) C IF (LTEST.AND.LNEW) *WRITE(6,17)(CI(K),K=1,7),(KI(K),K=6,25),(SIGMA(K),K=1,N1) 17 FORMAT('0CI:',7D11.4,/,' KI:',20I3,/,' SI:',5D11.4,/,59(4X,5D11.4/ *)) C DO 120 M = 1, NCOV IF (MT.EQ.2) COV(M,1) = (M-1)*DT RV = (M-1)*RT T = DCOS(RV) C RV IS EQUAL TO THE SPHERICAL DISTANCE BETWEEN P AND Q IN UNITS OF RAD- C IANS. U = DSIN(RV) IF (LPOLE) GO TO 105 SQ = U*CA CQ = DSQRT(D1-SQ *SQ) SD = U*SA/CQ CD = T/CQ GO TO 106 C 105 SD = D0 CD = D1 IF (RV.GT.PI/D2) CD = -D1 CQ = T I = 1 IF (ALFA.LT.D0) I = -1 SQ = U*I C C TRANSFER OF COORDINATE INFORMATION TO THE SUBROUTINE ACCORDING TO C THE SPECIFICATIONS GIVEN IN THE SUBROUTINE. 106 CR(1) = T CR(2) = HP(MT) CR(3) = HQ(MT) CR(4) = D0 CR(5) = SQ CR(6) = D1 CR(7) = CQ CR(8) = SD CR(9) = CD CR(10) = GM/(RE+HP(MT))**2 CR(11) = GM/(RE+HQ(MT))**2 IF (LTEST)WRITE(6,18)SQ,CQ,SD,CD 18 FORMAT('0SQ,CQ,SD,CD=',4D12.5) CALL COVCX(COV(M,MT)) C IF (.NOT.LTEST) GO TO 120 KK = KI(8)+1 WRITE(6,19)((CR(I*8+K+3),K=1,8),I=1,KK) 19 FORMAT(' CR:',8D11.4,/,4(4X,8D11.4,/)) 120 CONTINUE C IF (.NOT.(LAST3.OR.MT.EQ.7)) GO TO 104 C C OUTPUT OF A TABLE OF COVARIANCES. WRITE(7,30) WRITE(6,30) 30 FORMAT(' ') WRITE(6,20)IDEG,MIN,SEC WRITE(7,20)IDEG,MIN,SEC 20 FORMAT( ' TABLE OF COVARIANCES:',/, *' BETWEEN QUANTITIES OF KIND KP AND KQ, EVALUATED IN P,Q,',/, *' HAVING SPHERICAL DISTANCE PSI, HEIGHTS HP, HQ',/, *' AND AN AZIMUTH OF',I5,' D',I3,' M',F6.2,' SEC FROM P TO Q.') WRITE(6,30) WRITE(7,30) WRITE(6,21)(KP(I),I=2,MT) WRITE(7,21)(KP(I),I=2,MT) 21 FORMAT(' KP= ',6(I6,5X)) WRITE(6,22)(KQ(I),I=2,MT) WRITE(7,22)(KQ(I),I=2,MT) 22 FORMAT(' KQ= ',6(I6,5X)) WRITE(6,23)(HP(I),I=2,MT) WRITE(7,23)(HP(I),I=2,MT) 23 FORMAT(' HP= ',6(1X,F10.1)) WRITE(6,24)(HQ(I),I=2,MT) WRITE(7,24)(HQ(I),I=2,MT) 24 FORMAT(' HQ= ',6(1X,F10.1)) WRITE(6,26) WRITE(7,26) 26 FORMAT(' PSI') C DO 113 K = 1, NCOV RM = COV(K,1) IE = IDINT(RM/60.0D0) RM = RM-60.0D0*IE WRITE(7,25)IE,RM,(COV(K,I), I = 2, MT) 113 WRITE(6,25)IE,RM,(COV(K,I), I = 2, MT) 25 FORMAT(I4,F6.2,6F11.5) C MT = 1 IF (.NOT.LAST3) GO TO 104 C IF (.NOT.LAST2) GO TO 102 C IF (.NOT.LAST1) GO TO 100 STOP END SUBROUTINE RAD(RA,IDEG,MIN,SEC,MODE) IMPLICIT REAL*8(A-H,O-Z) DATA RS,PI/206264.806D0,3.1415926535D0/ IG = 1 IF (IDEG.LE.0 .AND. MODE.EQ.1) IG = -1 IDEG = IDEG*IG IF (MODE.NE.1 .OR. MIN .GE. 0) GO TO 20 IG = -1 MIN = MIN*IG C 20 GO TO (30,40,50,60),MODE 30 RA =IDEG*3600.0D0+MIN*60.0D0+SEC GO TO 70 40 RA =IDEG*3600.0D0+SEC*60.0D0 GO TO 70 50 RA = SEC*3600.0D0 GO TO 70 60 RA = SEC*3240.0D0 70 RA = RA/RS 80 IF (DABS(RA).LE.PI) GO TO 90 RA = RA-2.0D0*PI*DSIGN(RA,1.0D0) GO TO 80 90 RETURN END SUBROUTINE COVAX(SM) C THE SUBROUTINE COMPUTES THE COVARIANCE BETWEEN TWO QUANTITIES OF A C KIND SPECIPIED THROUGH THE VALUE OF TWO INTEGER VARIABLES (STORED IN C KI(6) AND KI(7), SEE BELOW). THE QUANTITIES ARE EVALUATED IN TWO C POINTS, P AND Q, THE COORDINATES OF WHICH ARE GIVEN IMPLICITLY BY THE C VALUES OF CR(1) - CR(9). C C THE COVARIANCE FUNCTION USED IS DEFINED ACCORDING TO A DEGREE-VARIANCE C MODEL AND A SET OF EMPIRICAL (POTENTIAL) DEGREE-VARIANCES. THE DEGREE- C VARIANCE MODEL IS SPECIFIED THROUGH THE VALUES OF KI(1)-KI(5),CI(8)- C CI(10) AND THE PARAMETERS N1 AND LOCAL OCCURRING IN THE COMMON BLOCK C /CMCOV/. EMPIRICAL ANOMALY DEGREE-VARIANCES WILL HAVE TO BE STORED IN C SIGMA0 WHEN LOCAL IS FALSE, AND ARE USED FOR THE COMPUTATION OF RESI- C DUAL POTENTIAL DEGREE-VARIANCES, (SEE REF(A), EQ.(16)). C C THE SUBROUTINE HAS THREE ENTRIES, COVAX, COVBX AND COVCX, WHICH HAVE C TO BE CALLED IN THIS SEQUENCE. C C BY THE CALL OF COVAX, THE KIND OF COVARIANCE FUNCTION TO BE USED IS C DETERMINED. THE VALUE OF KI(5) WILL DETERMINE THE DEGREE-VARI- C ANCE MODEL (1,2 OR 3, CF.REF(A),EQ.(17)) THAT WILL BE USED. THE QUAN- C TITIES K(2),K(3) MUST BE STORED IN KI(3),KI(4), AND BE EQUAL TO ZERO C WHEN NOT USED (EG.,KI(3),KI(4) BOTH ZERO WHEN KI(5)=1). THE QUANTITY C A(I) MUST BE STORED IN CI(8) IN UNITS OF (M/SEC)**4, AND THE SQUARE OF C THE RATIO BETWEEN THE RADIUS OF THE BJERHAMMAR-SPHERE (RB) AND THE C MEAN RADIUS OF THE EARTH (RE) MUST BE STORED IN CI(10). C C THERE ARE THEN THREE POSSIBILITIES: C (1) ONE OF THE DEGREE-VARIANCE MODELS IS USED WITHOUT MODIFICATIONS. C THE SUMMATION LIMIT P OF REF.(A),EQ.(20) IS THEN FIXED TO 3. C BECAUSE THIS IS EQUIVALENT TO REQUIRING THE FIRST 3 DEGREE-VARIAN- C AREA /CMCOV/ MUST BE EQUAL TO 3 AND .TRUE., RESPECTIVELY. C CES TO BE ZERO, THE VARIABLES N1 AND LOCAL STORED IN THE COMMON C (2) A NUMBER (N1) OF THE ANOMALY DEGREE-VARIANCES (DEGREE ZERO TO C N1-1) ARE PUT EQUAL TO EMPIRICAL DETERMINED QUANTITIES. THE ANO- C MALY DEGREE-VARIANCE OF DEGREE K WILL HAVE TO BE STORED IN C SIGMA0(K+1) IN UNITS OF MGAL**2 WHEN CALLING COVAX. LOCAL MUST BE C EQUAL TO .FALSE.. COVAX WILL CONVERT THE ANOMALY DEGRE5-VARIANCES C INTO POTENTIAL DEGREE-VARIANCES. C (3) THE N1 FIRST DEGREE-VARIANCES (DEGREE 0 - N1-1) ARE EQUAL TO ZERO. C THIS MEANS, THAT THE VALUES OF A (N1-1)-ORDER LOCAL COVARIANCE C FUNCTION WILL BE COMPUTED. LOCAL MUST HAVE THE VALUE .TRUE.. C IN ALL CASES N1 MUST BE LESS THAN 300. C C THE COVARIANCES WILL GENERALLY BE COMPUTED BY CLOSED EXPRESSIONS, BUT C THEY MAY IN CERTAIN CASES BE USELESS IN BIG ALTITUDES OF NUMERICAL C REASONS, CF. REF(A), SECTION 4. IN THEESE CASES MUST THE LOGICAL VARI- C ABLE LSUM BE TRUE AND THE VARIABLE HMAX MUST HAVE ASSIGNED A VALUE C EQUAL TO THE CRITICAL ALTITUDE. WHEN LSUM IS TRUE AND THE HEIGHT OF C P OR Q IS GREATHER THAN HMAX, WILL THE SERIES REF(A), EQ.(16), ABBRE- C VIATED TO DEGREE N2-1 BE USED FOR THE COMPUTATION OF THE COVARIANCES. C THE VALUES OF LSUM, N2 AND HMAX WILL (IN THE SAME WAY AS FOR THE PARA- C METERS SPECIFYING THE DEGREE-VARIANCE MODEL) BE TRANSFERRED TO COVAX C THROUGH THE COMMON AREA /CMCOV/, BUT AN ARRAY SM IS TRANSFERRED AS A C PARAMETER IN THE CALL IN ORDER TO ENABLE VARIABLE DIMENSIONING (SPECI- C FIED BY THE VARIABLE N2 IN /CMCOV/). C C THE CALL OF COVAX WILL ALSO INITIALIZE CERTAIN VARIABLES USED IN C SUBSEQUENT COMPUTATIONS. C C THE CALL OF COVBX WILL FIX CERTAIN CONSTANTS USED FOR THE COMPUTA- C TIONS, WHICH ARE INDEPENDENT OF THE POINTS P AND Q. WHEN COVBX IS CAL- C LED, THE KIND OF QUANTITIES BETWEEN WHICH THE COVARIANCE IS TO BE C COMPUTED MUST BE SPECIFIED. THIS IS DONE BY STORING IN KI(6) AND C KI(7) INTEGERS EQUAL TO THE EQUATION NUMBERS OF REF.A, EQ.(1) - (14) C DEFINING THE QUANTITIES. C C THE CALL OF COVCX WILL RESULT IN THE COMPUTATION OF THE COVARIANCE , C WHICH IS TRANSFERRED TO THE CALLING PROGRAM THROUGH THE VARIABLE COV. C INFORMATION RELATED TO THE COORDINATES OF P AND Q MUST BE STORED IN C THE ARRAY CR WHEN COVCX IS CALLED, SEE BELOW. C REFERENCES: C (A) TSCHERNING,C.C.: COVARIANCE EXPRESSIONS FOR SECOND AND LOWER ORDER C DERIVATIVES OF THE ANOMALOUS POTENTIAL, REPORTS OF THE DEP. OF C GEODETIC SCIENCE NO. 225,1975. C (B) TSCHERNING,C.C. AND R.H.RAPP: CLOSED COVARIANCE EXPRESSIONS C FOR GRAVITY ANOMALIES, GEOID UNDULATIONS, AND DEFLECTIONS OF C THE VERTICAL IMPLIED BY ANOMALY DEGREE-VARIANCE MODELS. DEP- C ARTMENT OF GEODETIC SCIENCE, THE OHIO STATE UNIVERSITY, C REPORT NO. 208, 1974. C IMPLICIT REAL *8(A-H,O-Z), LOGICAL (L) C COMMON /CMCOV/CI(12),CR(51),SIGMA0(300),SIGMA(300),HMAX,KI(25),N1, *N2,LOCAL,LSUM C THE COMMON BLOCK CONTAINS THE VALUES OF PARAMETERS USED FOR THE COM- C PUTATIONS AND RETURN VALUES OF FUNCTIONS AND CONSTANTS, WHICH HAVE C BEEN USED IN THE COMPUTATIONS. C PARAMETERS USED FOR THE COMPUTATIONS: C CI(8) = THE CONSTANT A(I) OF REF.(A), EQ.(17) IN UNITS OF (M/SEC)**4 C CI(10) THE SQUARE OF THE RATIO BETWEEN THE BJERHAMMAR-SPHERE RADIUS C (RB) AND THE MEAN RADIUS OF THE EARTH (RE). C CR(1) COSINE OF THE SPHERICAL DISTANCE BETWEEN P AND Q, C CR(2),CR(3) THE HEIGHT OF P, Q, RESPECTIVELY, (UNITS METERS), C CR(4),CR(5) SINE OF THE LATITUDE OF P, Q, RESPECTIVELY, C CR(6),CR(7) COSINE OF THE LATITUDE OF P, Q, RESPECTIVELY, C CR(8),CR(9) SINE AND COSINE OF THE LONGITUDE DIFFERENCE, C CR(10),CR(11) THE REFERENCE GRAVITY IN P, Q, RESPECTIVELY (WHEN C USED, OTHERWISE STORE 1.0D0), (UNITS M/SEC**2). C SIGMA0(1)-SIGMA0(N1) MUST CONTAIN THE EMPIRICAL ANOMALY DEGREE- C VARIANCES IN UNITS OF MGAL**2. C KI(3) = K(2) OF DEG.VAR. MODEL 2 OR 3, C KI(4) = K(3) OF DEG.VAR. MODEL 3, CF. REF.(A), EQ.(17). C KI(5) = THE DEG.VAR. MODEL NUMBER, (EQUAL TO 1, 2 OR 3), C KI(6),KI(7) THE INTEGER SPECIFYING THE KIND OF QUANTITY WHICH IS C ASSOCIATED WITH P, Q, RESPECTIVELY, C N1 = THE NUMBER OF EMPIRICAL DEGREE-VARIANCES USED (LOCAL =.FALSE.) C OR (ORDER+1) OF THE LOCAL COVARIANCE FUNCTION USED (LOCAL=.TRUE.). C HMAX, N2, LSUM. HMAX IS THE HEIGHT ABOVE WHICH THE LEGENDRE SERIES C OF MAXIMAL DEGREE N2-1 WILL BE USED FOR THE COMPUTATION OF THE CO- C VARIANCES WHEN LSUM IS TRUE. N2 MUST BE GREATHER THAN 2 AS WELL AS C GREATHER THAN N1. C RETURN VALUES: C CI(1)-CI(7), THE QUANTITIES C(J,Q) OF REF.(A), EQ.(47), WITH C CI(1) - CI(KI(5)+1) = C(J,Q), CI(5) = C(KI(5)+2,Q), C CI(6) = C(KI(5)+3,Q), CI(7) = C(KI(5)+4,Q), C CI(9) = RB**2, CI(11),CI(12) QUANTITIES USED TO GIVE THE COMPUTED C COVARIANCES THE PROPER UNITS. C CR(ND*8+12), THE VALUES OF THE ND'TH DERIVATIVE OF THE SUM OF THE C FINITE LEGENDRE-SERIES, CF.REF.(A), EQ.(20),(48) AND (52). C CR(ND*8+13) - CR(ND*8+19), THE VALUES OF THE ND'TH DERIVATIVES OF C THE FUNCTIONS F(-2), F(-1), F(KI(3)), F(KI(4)), S0, S1, S2, CF. REF. C (A), EQ. (42), (41), (39), (39), (30), (34) AND (35). C SIGMA0(1) - SIGMA0(N1) THE POTENTIAL DEGREE-VARIANCE CORRECTIONS, C CF. REF.(A), EQ.(16), (AFTER THE CALL OF COVAX). C SIGMA(4) - SIGMA(N1), THE POTENTIAL DEGREE-VARIANCES MULTIPLIED BY C THE FACTORS GIVEN IN REF.(A), TABLE 1. C SIGMA(1) - SIGMA(3), THE DEGREE-VARIANCES OF DEGREE 0,1,2 MINUS C TERMS OF THE SAME DEGREES ACQUIRED FROM REF.(A), EQ.(34),(35),(41) C AND (42). C KI(8),KI(9) THE NUMBER OF DIFFERENTIATIONS IN RADIAL DIRECTION AND C WITH RESPECT TO T = COS(SPHERICAL DIST.) TO BE PERFORMED. C KI(10) - KI(15) THE CONSTANTS I,K,J,M,J1,M1 OF REF.(A), SECTION 2. C KI(16) - KI(19) THE QUANTITIES M(1) - M(4) OF REF.(A), EQ.(26)-(29). C KI(20),KI(21) THE EXPONENT OF THE REFERENCE GRAVITY, C KI(22),KI(23) THE EXPONENT OF THE RADIAL DISTANCE AND C KI(24),KI(25) THE EXPONENT OF COSINE OF THE LATITUDE, OF P, Q RES- C PECTIVELY WITH WHICH THEESE QUANTITIES ARE USED IN THE COVARIANCE C COMPUTATIONS. C DIMENSION K7(15),K9(15),K11(15),K13(15),K15(15),K17(15),K19(15), *K21(15),K23(15),C11(15),K8(15),L(7),LN(7),CX(6,8),SM(N2) *,C(6),V(6),U(6),G(6),P(6),R(6),SS1(4),D(36),RM(6),Q(6) C THE ARRAY SM IS USED TO STORE THE DEGREE-VARIANCES WHEN THE LOGICAL C VARIABLE LSUM IS TRUE. IN CASE THE SUBSCRIPT LIMIT IS CHANGED IS IT C NECESSARY TO CHANGE THE VALUE OF THE VARIABLE N2 ACCORDINGLY. C EQUIVALENCE (CX(1,1),C(1)),(CX(1,2),V(1)),(CX(1,3),U(1)), *(CX(1,4),G(1)),(CX(1,5),P(1)),(CX(1,6),R(1)),(CX(1,7),SS1(1)), *(CX(2,8),SS2) C DATA D0,D1,D2,D3,RE/0.0D0,1.0D0,2.0D0,3.0D0,6371.0D3/ *,K7/5*0,6*1,3*2,0/,K9/5*1,2,3,2,3,2,3,2,2,3,1/,K11/11*0,2,3,3,0/, *K13/11*1,2,3,3,1/,K15/0,1,-1,-1,1,0,0,-1,-1,1,1,0,0,0,0/, *K17/3*0,2,2,10*0/,K19/1,4*0,1,1,8*0/,K21/0,2,1,2,2,1,1,7*2,0/, *K23/6*0,1,0,1,0,1,0,1,2,0/,K8/0,1,1,2,2,0,0,4*1,4*0/, *C11/1.0D0,1.0D9,1.0D5,2*1.0D9,2*-206264.806D0,7*1.0D9,1.0D0/ C THE ARRAYS K7 - K23 CONTAINS TABLES OF QUANTITIES RELATED TO THE KIND C OF COVARIANCES (1 - 14) WHICH MAY BE COMPUTED. THE ELEMENT WITH SUBS- C CRIPT 15 IS DUMMY (RESERVED FOR PROGRAM EXTENSIONS). THEIR ACTUAL VA- C LUES WILL AFTER CALL OF COVBX BE STORED IN THE ELEMENTS OF THE ARRAY C KI HAVING SUBSCRIPTS 8 - 25. C K7 CONTAINS THE ORDER OF DIFFERENTIATION WITH RESPECT TO T,K8 THE C ORDER OF DIFFERENTIATION WITH RESPECT TO THE RADIUS, GF.REF(A),TABLE C 1. K9,K11,K13 THE KIND OF DIFFERENTIATIONS TO BE COMPUTED WITH RESPECT C TO THE LATITUDE (2) AND THE LONGITUDE (3), CF.REF(A),SECTION 3. K15 C AND K17 CONTAINS AN INTEGER, WHICH WILL BE ADDED TO THE DEGREE. THE C SUM WILL THEN BE MULTIPLIED WITH THE DEGREE-VARIANCE OF THE CORRESPON- C DING DEGREE WHEN A FIRST AND/OR SECOND DIFFERENTIATION WITH RESPECT C TO THE RADIAL DISTANCE HAS TAKEN PLACE. C C11 CONTAIN QUANTITIES USED TO GIVE THE COVARIANCES THE PROPER UNITS. C KT = KI(5) KT1 = KT+1 IF (KT.LT.3) GO TO 15 DO 16 K = KT, 2 16 KI(K+2) = D0 15 KI(1) = -2 KI(2) = -1 C IF ((KT.LT.3).OR.(KT.EQ.3.AND.KI(4).GT.KI(3))) GO TO 17 C ASSURING, THAT KI(4).GT.KI(3), BECAUSE THIS FACT IS USED IN SUB- C SEQUENT COMPUTATIONS. K = KI(3) KI(3) = KI(4) KI(4) = K 17 II = KI(3) JJ = KI(4) SM(1) = D0 SM(2) = D0 N3 = N1 A = CI(8) S = CI(10) RB2 = S*(RE**2) CI(9) = RB2 RB2 = RB2*1.0D-10 T = D0 Q(1) = D0 RM(1) = D0 C SIGMA0(1) = D0 SIGMA0(2) = D0 IF (LOCAL) SIGMA0(3) = D0 IF (.NOT.LOCAL) SIGMA0(3) = SIGMA0(3)*RB2/S**4 IF (N1.LT.4) GO TO 14 DO 13 K = 4, N1 IF (.NOT.LOCAL) T = SIGMA0(K)*S**(-K-1)*RB2 GO TO (10,11,12),KT 10 KK = 1 GO TO 13 11 KK = K+II-1 GO TO 13 12 KK = (K+II-1)*(K+JJ-1) 13 SIGMA0(K) = (T-A*(K-2)/((K-3)*KK))/(K-2)**2 14 RETURN C C ENTRY COVBX C BY THE CALL OF COVBX ALL QUANTITIES NECESSARY FOR THE COMPUTATION OF C THE COVARIANCE, BUT INDEPENDENT OF THE POSITION OF THE POINTS P AND Q, C ARE COMPUTED. RB2 = CI(9) CI(11) = D1 IF (KI(6).EQ.15.OR.KI(7).EQ.15) GO TO 19 C DO 20 M = 1, 2 K = KI(M+5) C FOR M = 1, K IS EQUAL TO THE KIND EVALUATED IN P AND FOR M = 2 EQUAL C TO THE KIND EVALUATED IN Q. KI(M+9) = K9(K) KI(M+11) = K11(K) KI(M+13) = K13(K) KI(M+15) = K15(K) KI(M+17) = K17(K) KI(M+19) = K19(K) KI(M+21) = K21(K) KI(M+23) = K23(K) C 20 CI(11) = CI(11)*C11(K) C KQ = K KP = KI(6) KI(8) = K7(KP)+K7(KQ) KI(9) = K8(KP)+K8(KQ) 19 ND = KI(8) NR = KI(9) C ND AND NR ARE THE NUMBER OF DIFFERENTIATIONS WITH RESPECT TO T AND THE C RADIAL DISTANCES, RESPECTIVELY. C C UPDATING THE DEGREE-VARIANCES, CF. REF(A), TABLE 1. SIGMA(1) = D0 SIGMA(2) = D0 IF (LSUM) N1 = N2 DO 21 M = 3, N1 B = D1 DO 22 I = 1, 4 22 IF (KI(I+15).NE.0) B = B*(M+KI(I+15)-1) IF (M.LE.N3) SIGMA(M) = SIGMA0(M)*B IF (.NOT.LSUM.OR.M.EQ.3) GO TO 21 DO 48 K = 1, KT1 48 B = B/(M+KI(K)-1) C STORING THE MODIFIED DEGREE-VARIANCES OF DEGREE M-1 IN SM(M) AND AD- C DIND THE DEGREE-VARIANCE CORRECTIONS FOR M .LE. N3. SM(M) = B*A IF (M.LE.N3) SM(M) = SM(M)+SIGMA(M) 21 CONTINUE IF (N1.GT.2) SM(3) = SIGMA(3) IF (LSUM) N1 = N3 C C EVALUATION OF THE QUANTITIES C(J,NR), CF.REF(A), TABLE 2. DO 23 K = 1, 7 23 CI(K) = D0 C DO 25 K = 1, KT1 CI(K) = D1 DO 25 KQ = 1, KT1 25 IF (K.NE.KQ) CI(K) = CI(K)/(KI(KQ)-KI(K)) C CF.REF(A),EQ.(19). WE WILL THEN COMPUTE THE QUANTITIES GIVEN IN REF(A) C TABLE 2. IF (NR.LT.2) GO TO 29 KP = KI(16)+KI(17)+KI(18)+KI(19) IF (NR.EQ.4) M = KI(16)*(KI(17)+KI(18)+KI(19))+KI(17)*(KI(18)+ *KI(19))+KI(18)*KI(19) C GO TO (26,27,28),KT 26 CI(NR+3) = D1 IF (NR.GT.2) CI(NR+2) = KP+3 IF (NR.EQ.4) CI(NR+1) = M+3*KP+7 GO TO 29 27 IF (NR.GT.2) CI(NR+2) = D1 IF (NR.EQ.4) CI(NR+1) =-KI(3)+3+KP GO TO 29 28 IF (NR.EQ.4) CI(NR+1) = D1 29 IF (NR.EQ.0) GO TO 31 C DO 30 KP = 1, 4 DO 30 K = 1, KT1 30 IF (KI(KP+15).NE.0) CI(K) = CI(K)*(KI(KP+15)-KI(K)) C C THE LOGICAL ARRAYS L AND LN REGISTER WHICH TERMS THAT WILL HAVE TO C BE EVALUATED , RESPECTIVELY NOT EVALUATED IN REF.(A), EQ. (47). 31 DO 38 K = 1, 7 L(K) = DABS(CI(K)).GT.1.0D-15 38 LN(K) = .NOT.(L(K)) C DO 32 K = 3, 7 DO 32 M = 1, 3 IF (M.EQ.1.AND.K.GT.5.OR.(M+KI(K)-1).EQ.0.AND.K.LT.5.OR.LN(K)) *GO TO 32 GO TO (34,34,35,35,34,36,37),K 34 B = D1 GO TO 33 35 B = D1/(M+KI(K) -1) GO TO 33 36 B = (M-1) GO TO 33 37 B = (M-1)*(M-1) 33 SIGMA(M) = SIGMA(M)-A*CI(K)*B 32 CONTINUE SIGMA(3) = SIGMA(3)-A*CI(2) C ND1 = ND+1 ND2 = ND+2 RETURN C ENTRY COVCX(COV) C COMPUTATION OF THE COVARIANCE IN A SPECIFIC PAIR OF POINTS. THE VALUE C IS RETURNED THROUGH THE PARAMETER COV. C THE COVARIANCES COMPUTED WILL BE IN UNITS CORRESPONDING TO THE KIND OF C QUANTITIES, I.E. FOR KIND (1) METERS, (2) EOTVOS (E), (3) MGAL, C (4),(5) E, (6),(7) ARCSECONDS, (8) - (14) E. C THE FOLLOWING QUANTITIES MUST BE STORED IN THE ELEMENTS OF THE ARRAY C CR WHEN COVCX IS CALLED: (1) COSINE TO THE SPHERICAL DISTANCE BETWEEN C P AND Q, (2),(3) THE HEIGHT OF P, Q RESPECTIVELY, (4),(5) SINE OF THE C THE LATITUDE OF P, Q, RESPECTIVELY, (6),(7) COSINE OF THE LATITUDE OF C P, Q, RESPECTIVELY, (8),(9) SINE AND COSINE OF THE LONGITUDE DIFFER- C ENCE. THE REFERENCE GRAVITY WILL HAVE TO BE STORED IN CR(10),CR(11) C FOR P, Q RESPECTIVELY (WHEN USED, OTHERWISE STORE 1.0). T = CR(1) HP = CR(2) HQ = CR(3) SP = CR(4) SQ = CR(5) CP = CR(6) CQ = CR(7) SD = CR(8) CD = CR(9) RP = RE+HP RQ = RE+HQ C IN HEIGH ALTITUDES AND WHEN LSUM IS TRUE WILL THE COVARIANCE BE COM- C PUTED BY A SUMMATION OF THE LEGENDRE-SERIES ABBREVIATED TO DEGREE C N2-1. LSUMC = LSUM .AND. (HP.GT.HMAX .OR. HQ.GT.HMAX) C COMPUTATION OF THE CONSTANT USED TO CONVERT THE COVARIANCE INTO C PROPER UNITS. CI(12) = CI(11)/(CP**KI(24)*CQ**KI(25)*RP**KI(22)*RQ**KI(23) **CR(11)**KI(21)*CR(10)**KI(20)) C S = RB2/(RP*RQ) S2 = S*S ST = S*T T2 = T*T P2 = (D3*T2-D1)/D2 P3 = (D3*ST+D1)/D2 C C INITIALIZING ARRAY ELEMENTS. NOTE THE USE OF THE EQUIVALENCING. DO 50 K = 1, 8 DO 50 M = 1, ND2 50 CX(M,K) = D0 DO 51 K = 1, ND2 C(K) = D0 51 D(K) = D0 DO 52 K = 1, 40 52 CR(K+11) = D0 C C SUMMATION AND DIFFERENTIATION OF THE LEGENDRE SERIES, CF.REF(A),EQ. C (49) AND (51). IF (LSUMC) N1 = N2 K1 = N1 K2 = N1+1 K = N1-1 DO 54 M = 1, N1 GI = (D2*K+D1)*S/K1 GJ = -K1*S2/K2 K2 = K1 K1 = K K = K-1 IF (.NOT.LSUMC) SI = SIGMA(K2) IF (LSUMC) SI = SM(K2) I2 = 0 I1 = 1 DO 53 I = 2, ND2 B = D(I) D(I) = C(I) C(I) = GI*(D(I)*T+I2*D(I1))+GJ*B+SI SI = D0 I2 = I1 53 I1 = I 54 CONTINUE IF (LSUMC) N1 = N3 C C COMPUTATION OF THE QUANTITIES D(1)-D(36),CF.REF(A),SECTION 3. IF (ND.EQ.0) GO TO 55 C D(1) = D1 CS = CP*SQ SC = SP*CQ SCC = SC*CD CC = CP*CQ CCS = CC*SD CSC = CS*CD D(2) = CS-SCC D(3) = CCS D(7) = SC-CSC D(13) = -CCS IF (ND.EQ.1) GO TO 55 C D(4) = -T SCS = SC*SD CCC = CC*CD SS = SP*SQ SSC = SS*CD CSS = CS*SD D(5) = -SCS D(6) = -CCC D(8) = CC+SSC D(9) = -CSS D(14) = SCS D(15) = CCC D(19) = -T D(25) = CSS D(31) = -CCC IF (ND.EQ.2) GO TO 55 C SSS = SS*SD D(10) = -SC+CSC D(11) = SSS D(12) = CSC D(16) = CCS D(17) = -SCC D(18) = CCS D(20) = -CS+SCC D(21) = -CCS D(26) = -SSS D(27) = -CSC D(32) = SCC D(33) = -CCS IF (ND.EQ.3) GO TO 55 C D(22) = T D(23) = SCS D(24) = CCC D(28) = -CSS D(29) = SSC D(30) = -CSS D(34) = CCC D(35) = SCS D(36) = CCC 55 IF (LSUMC) GO TO 75 C C COMPUTATION OF THE FUNCTIONS L=R(1), N=1/RN, M=RM(2), F0=P(2), CF. C REF.(A), EQ. (31)-(33),(40) AND (77A). RL2 = D1-D2*ST+S2 RL = DSQRT(RL2) R(1) = RL RL1 = D1/RL RN = D1/(D1+RL-ST) RL2 = D1/RL2 RNL = RN*RL1 RM(2) = D1-RL-ST P(2) = S*DLOG(D2*RN) RL3 = RL2*RL1 RL5 = RL3*RL2 S3 = S2*S R(2) = -S*RL1 IF (ND.EQ.0) GO TO 56 C C COMPUTATION OF THE DERIVATIVES WITH RESPECT TO T. C CF. REF.(A), EQ. (77B),(69A),(57). R(3) = -S2*RL3 RM(3) = -R(2)-S P(3) = S2*(RNL+RN) IF (ND.EQ.1) GO TO 56 C C CF. REF.(A), EQ. (77C),(69B),(58). R(4) = -D3*S3*RL5 RM(4) = -R(3) P(4) = S3*(RL3+(D1+(D2+RL1)*RL1)*RN)*RN IF (ND.EQ.2) GO TO 56 C C CF. REF.(A), EQ. (77D),(69C),(59). RL4 = RL2*RL2 RL7 = RL5*RL2 S4 = S2*S2 R(5) = -15.0D0*S4*RL7 RM(5) = -R(4) P(5) = S4*(D3*RL5+((D3+D3*RL1)*RL3+D2*(D1+(D3+(D3+RL1)*RL1)*RL1) **RN)*RN)*RN IF (ND.EQ.3) GO TO 56 C C CF. REF.(A), EQ. (69D),(60). S5 = S4*S RL6 = RL4*RL2 RM(6) = -R(5) P(6) = D3*S5*((5.0D0*RL7+((4.0D0+5.0D0*RL1)*RL5+((4.0D0+(8.0D0 *+4.0D0*RL1)*RL1)*RL3+(2.0D0+(8.0D0+(12.0D0+(8.0D0+D2*RL1)*RL1) **RL1)*RL1)*RN)*RN)*RN)*RN) C 56 IF (LN(2)) GO TO 58 C COMPUTATION OF THE FUNCTIONS L=R(1), N=1/RN, M=RM(2), F0=P(2), CF. C REF.(A), EQ. (31)-(33),(40) AND (77A). RL2 = D1-D2*ST+S2 RL = DSQRT(RL2) R(1) = RL RL1 = D1/RL RN = D1/(D1+RL-ST) RL2 = D1/RL2 RNL = RN*RL1 RM(2) = D1-RL-ST P(2) = S*DLOG(D2*RN) RL3 = RL2*RL1 RL5 = RL3*RL2 S3 = S2*S R(2) = -S*RL1 IF (ND.EQ.0) GO TO 56 C C COMPUTATION OF THE DERIVATIVES WITH RESPECT TO T. C CF. REF.(A), EQ. (77B),(69A),(57). R(3) = -S2*RL3 RM(3) = -R(2)-S P(3) = S2*(RNL+RN) IF (ND.EQ.1) GO TO 56 C C CF. REF.(A), EQ. (77C),(69B),(58). R(4) = -D3*S3*RL5 RM(4) = -R(3) P(4) = S3*(RL3+(D1+(D2+RL1)*RL1)*RN)*RN IF (ND.EQ.2) GO TO 56 C C CF. REF.(A), EQ. (77D),(69C),(59). RL4 = RL2*RL2 RL7 = RL5*RL2 S4 = S2*S2 R(5) = -15.0D0*S4*RL7 RM(5) = -R(4) P(5) = S4*(D3*RL5+((D3+D3*RL1)*RL3+D2*(D1+(D3+(D3+RL1)*RL1)*RL1) **RN)*RN)*RN IF (ND.EQ.3) GO TO 56 C C CF. REF.(A), EQ. (69D),(60). S5 = S4*S RL6 = RL4*RL2 RM(6) = -R(5) P(6) = D3*S5*((5.0D0*RL7+((4.0D0+5.0D0*RL1)*RL5+((4.0D0+(8.0D0 *+4.0D0*RL1)*RL1)*RL3+(2.0D0+(8.0D0+(12.0D0+(8.0D0+D2*RL1)*RL1) **RL1)*RL1)*RN)*RN)*RN)*RN) C 56 IF (LN(2)) GO TO 58 C COMPUTATION OF THE FUNCTION F-1 AND ITS DERIVATIVES, CF. REF.(A), C EQ. (41) AND (61) - (65). U(2) = S*(RM(2)+T*P(2)) IF (ND2.LT.3) GO TO 58 DO 57 K = 3, ND2 57 U(K) = S*(RM(K)+T*P(K)+(K-2)*P(K-1)) C 58 IF (LN(1)) GO TO 60 C COMPUTATION OF THE FUNCTION F-2 AND ITS DERIVATIVES, CF. REF.(A) EQ. C (42), AND (65)- (68). DO 59 K = 2, ND2 GO TO (61,61,62,63,64,65),K 61 CY = S*(D1-T2)/4.0D0 GO TO 59 62 CY = -ST/D2 GO TO 59 63 CY = D3*P(2)-S/D2 GO TO 59 64 CY = 9.0D0*P(3) GO TO 59 65 CY = 18.0D0*P(4) 59 V(K) = S*(RM(K)*P3+S*((K-2)*D3*RM(K-1)/D2+P2*P(K)+D3*T*P(K-1)* *(K-2)+CY)) C 60 IF (LN(3)) GO TO 73 C COMPUTATION OF THE FUNCTION F1 AND ITS DERIVATIVES, CF. REF.(A) EQ. C (36), REF.(B), EQ.(101) AND REF.(A), EQ.(70),(71). Q(2) = DLOG(D1+D2*S/(D1-S+RL)) IF (ND.EQ.0) GO TO 66 Q(3) = S2*RNL IF (ND.EQ.1) GO TO 66 Q(4) = S3*((RL1+D1)*RN+RL2)*RNL IF (ND.EQ.2) GO TO 66 Q(5) = S4*(D3*RL4+((D2+D3*RL1)*RL2+(D2 +(4.0D0+D2*RL1)*RL1)*RN) **RN)*RNL IF (ND.EQ.3) GO TO 66 Q(6) = D3*S5*(5.0D0*RL6+((D3+5.0D0*RL1)*RL4+((D2+(6.0D0+4.0D0*RL1) **RL1)*RL2+(6.0D0+(6.0D0+D2*RL1)*RL1)*RNL)*RN)*RN)*RNL C C COMPUTATION OF THE FUNCTION F2 AND ITS DERIVATIVES, CF. REF.(A), EQ. C (3),(72)-(75). 66 P(2) = (RL-D1+T*Q(2))/S IF (ND.EQ.0) GO TO 68 DO 67 K = 3, ND2 67 P(K) = (R(K-1)+T*Q(K)+(K-2)*Q(K-1))/S 68 I1 = II-1 K1 = 1 J1 = I1 IF (I1.GE.2) GO TO 49 DO 49 M = 2, ND2 IF (I1.EQ.0) G(M) = Q(M) IF (I1.EQ.1) G(M) = P(M) 49 CONTINUE IF (L(4)) J1 = JJ-1 IF (J1.LE.1) GO TO 71 C C CF. REF.(A), EQ. (38),(76). DO 71 K = 2, J1 DO 69 M = 2, ND2 B = Q(M) Q(M) = P(M) 69 P(M) = (R(M-1)+(2*K-1)*((M-2)*Q(M-1)+T*Q(M))-K1/S*B)/(K*S) IF (K.NE.I1) GO TO 71 DO 70 M = 2, ND2 70 G(M) = P(M) 71 K1 = K C 73 IF (LN(6)) GO TO 72 C CF. REF.(A), EQ. (34),(55). SS1(2) = S2*(T-S)*RL3 IF (ND.GT.0) SS1(3) = S2*(RL3+D3*(T-S)*S*RL5) C C CF. REF.(A), EQ. (35). 72 IF (L(7)) SS2= S2*((T+S)*RL3+D3*S*(T2-D1)*RL5) C C ADDING THE DIFFERENT TERMS, CF. REF.(A), EQ. (22),(47). C TIPLIED BY RB**2 IN UNITS OF MGAL**2, THE INTEGERS K(2),K(3) OF EQ. 75 DO 78 M = 2, ND2 C CF. REF.(A), EQ. (50),(52). C(M) = S*C(M) CR(M*8 -4) = C(M) DO 78 K = 1, 7 IF (LN(K)) GO TO 78 C STORING THE TERMS FOR TRANSFER TO THE CALLING PROGRAM USING THE COMMON C AREA /CMCOV/. CR(M*8+K -4) = A*CX(M,K+1)*CI(K) IF (K.EQ.5) CR(M*8+K-4) = -CR(M*8+K-4) C(M) = C(M)+CR(M*8+K -4) 78 CONTINUE C C INTEGERS SPECIFYING THE KINDS OF DIFFERENTIATION WITH RESPECT TO THE C LATITUDES AND/OR THE LONGITUDES, CF. REF.(A), SECTION 3. I = KI(10) J = KI(12) K = KI(11) M = KI(13) J1 = KI(14) M1 = KI(15) C C COMPUTATION OF THE DERIVATIVES OF ORDER ND WITH RESPECT TO THE LATI- C TUDES AND THE LONGITUDES, CF. REF.(A), EQ. (43) - (46). GO TO (80,81,82,83,84),ND1 80 COV = C(2) GO TO 85 81 COV = C(3)*D(I+6*(K-1)) GO TO 85 82 COV = D(I)*D(J1)*D(6*(K-1)+1)*D(6*(M1-1)+1)*C(4)+D(I+J+6*(K+M-1)) **C(3) GO TO 85 83 COV = D(I+J+6*(K+M-1))*C(3)+(D(I+J)*D(6*(K+M-1)+1)+D(I+6*(K-1)) **D(J1+6*(M1-1))+D(I+6*(M1-1))*D(J1+6*(K-1)))*C(4) *+D(I)*D(J1)*D(6*(K-1)+1)*D(6*(M1-1)+1)*C(5) GO TO 85 84 COV = D(I+J+6*(K+M-1))*C(3)+(D(I+J+6*(K-1))*D(6*(M-1)+1) *+D(I+6*(M-1+K))*D(J)+D(J+6*(K+M-1))*D(I)+D(I+J+6*(M-1)) **D((K-1)*6+1)+D(I+J)*D(6*(K+M-1)+1)+D(I+6*(K-1))*D(J+6*(M-1)) *+D(I+6*(M-1))*D(J+6*(K-1)))*C(4)+(D(I+J)*D(6*(K-1)+1)*D(6*(M-1)+1) *+D(I+6*(K-1))*D(J)*D(6*(M-1)+1)+D(I+6*(M-1))*D(J)*D(6*(K-1)+1) *+D(J+6*(K-1))*D(I)*D(6*(M-1)+1)+D(J+6*(M-1))*D(I)*D(6*(K-1)+1) *+D(6*(K+M-1)+1)*D(I)*D(J))*C(5)+D(I)*D(J)*D(6*(K-1)+1)*D(6*(M-1) *+1)*C(6) C C GIVING THE COVARIANCE THE PROPER UNITS. 85 COV = COV*CI(12) C RETURN END