C MICROSOFT FORTRAN SOURCE, "CLENQT.FOR" (NIMA/GIMG, 10 OCT 1996) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C*********THE PARAMETER VALUES ARE AS FOLLOWS: C C*******NMAX DENOTES MAXIMUM DEGREE & ORDER ATTAINED. C*******ND = NUMBER OF C OR S GEOPOTENTIAL COEFFICIENTS. C*******NDD = SIZE OF A & B CLENSHAW ARRAYS. C PARAMETER (NMAX=180,NP1=NMAX+1,NP2=NMAX+2) PARAMETER (ND=NP1*NP2/2) PARAMETER (NP3=NMAX+3,NP4=NMAX+4,NDD=NP3*NP4/2) C C********THE ARRAYS ARE AS FOLLOWS: C C********CNM & SNM STORE THE GEOPOTENTIAL COEFFICIENTS. C********IV IS A LOCATING ARRAY. C********ALL OTHER ARRAYS STORE CLENSHAW COEFFICIENTS NEEDED FOR C********THE SELECTED GRAVIMETRIC QUANTITIES THAT ARE COMPUTED. C DIMENSION T11(0:NP2),TG1(0:NP2),TP1(0:NP2) DIMENSION T12(0:NP2),TG2(0:NP2),TP2(0:NP2) DIMENSION S11(0:NP2),S12(0:NP2),SP1(0:NP2),SP2(0:NP2) DIMENSION SV15(0:NP2),SV1(0:NP2),SV2(0:NP2),RNN(0:NMAX) DIMENSION A(NDD),B(NDD),CNM(ND),SNM(ND),IV(NP3) DIMENSION AS(0:NP2),CR(0:NMAX),SR(0:NMAX),S2(0:NP2) DIMENSION SG1(0:NP2),SG2(0:NP2),SG(0:NP2),SHT(0:NP2) DIMENSION C2N(5) C C**********WGS84 MODEL VALUES ARE USED. C DATA C2/108262.9989050D-8/,RKM/3.9860015D14/ DATA AE/6378137.D0/,ESQ/.00669437999013D0/,F/298.257D0/ DATA GRAVA/9.7803267714D0/,STAR/.001931851386D0/ C C**********WGS72 MODEL VALUES ARE USED. C C DATA C2/108263.0D-8/,RKM/3.986005D14/ C DATA AE/6378135.D0/,ESQ/.006694317778D0/,F/298.26D0/ C DATA GRAVA/9.7803327D0/,STAR/.005278994D0/ C OPEN (12, ERR=350, FILE='EGM180.NOR', STATUS='OLD') TEMPP = 400.D0 TEMPL = 400.D0 HGHT = .0000001D0 GAMMAA = RKM/AE**2 GA = RKM/AE F = 1.D0/F T1 = 1.D0 T2 = 2.D0 T3 = 3.D0 PI = DACOS(-1.D0) PI2 = PI/2.D0 RPD = PI/180.D0 RPS = RPD/3600.D0 RHO = 1.D0/RPS C C********READ IN POTENTIAL COEFFICIENTS. C PRINT*,'READING COEFFICIENT FILE' 400 CONTINUE READ(12,FMT=431,END=440) N, M, CBAR, SBAR 431 FORMAT(2I5,2D15.8) IF(N .GT. NMAX)GO TO 440 LL = N*(N+1)/2 + M + 1 CNM(LL) = CBAR SNM(LL) = SBAR C GO TO 400 440 CONTINUE PRINT*,'COMPLETED READING COEFFICIENTS' PRINT*,' ' C SR(0) = 0.D0 CR(0) = 1.D0 C C********BUILD LOCATING ARRAY. C DO 1 I = 1,NP3 IV(I) = I*(I-1)/2 1 CONTINUE DO 21 I = 0,NMAX RNN(I) = I - 1 21 CONTINUE C C*******MODIFY CNM EVEN ZONAL COEFFICIENTS C DO 22 I=2,5 C2N(I) = (-1.D0)**(I+1) * (3.D0*(ESQ**I)/((2.D0*I+1.D0)* & (2.D0*I+3.D0))) * (1.D0 - I + 5.D0*I*C2/ESQ) 22 CONTINUE C C*******WGS84 VALUES ARE USED. CB2 = -1.D0 * C2/DSQRT(5.D0) CB4 = -1.D0 * C2N(2)/3.D0 CB6 = -1.D0 * C2N(3)/DSQRT(13.D0) CB8 = -1.D0 * C2N(4)/DSQRT(17.D0) CB10 = -1.D0 * C2N(5)/DSQRT(21.D0) C C*******WGS72 VALUES ARE USED. C CB2 = -4.841732D-04 C CB4 = 7.8305D-07 C CB6 = 0.0 C CB8 = 0.0 C CB10 = 0.0 C CNM(4) = CNM(4) - CB2 CNM(11) = CNM(11) - CB4 CNM(22) = CNM(22) - CB6 IF( NMAX .GT. 6 )CNM(37) = CNM(37) - CB8 IF( NMAX .GT. 9 )CNM(56) = CNM(56) - CB10 C C********BUILD ALL CLENSHAW COEFFICIENT ARRAYS. C DO 2 I = 1,NMAX RI = I AS(I) = - DSQRT((T2*RI+T1)/(T2*RI)) 2 CONTINUE C DO 6 M = 0,NMAX MP1 = M + 1 RM = M DO 5 N = MP1,NMAX RN = N LL = IV(N+1)+M+1 A(LL) = -DSQRT(((T2*RN+T1)*(T2*RN-T1))/((RN-RM)* & (RN+RM))) B(LL) = DSQRT(((T2*RN+T1)*(RN+RM-T1)*(RN-RM-T1))/ & ((T2*RN-T3)*(RN+RM)*(RN-RM))) 5 CONTINUE 6 CONTINUE C C********PROMPT USER FOR COMPUTATION POINT DATA. C 66 CONTINUE WRITE(6,333) 333 FORMAT(1X,'ENTER LAT,LONG IN DEGREES FOLLOWED BY GEODETIC HT.') WRITE(6,334) 334 FORMAT(1X,'IN METERS. PUT ALL 3 ON SAME LINE. HAVE LAT. BETWEEN', . /,1X,'-90 & + 90 AND LONGITUDE BETWEEN 0 & 360 DEGREES',//) READ (5,*) PHID, RLAMD, HT PHI = PHID*RPD C C********PERFORM CLENSHAW SUMMATION FOR EACH PT. C*********COMPUTE GEOCENTRIC DISTANCE,COLATITUDE AND NORMAL GRAVITY. C IF(TEMPP .EQ. PHID .AND. HGHT .EQ. HT)GO TO 67 C RN = AE/DSQRT(1.D0-ESQ*DSIN(PHI)**2) T22 = (RN+HT)*DCOS(PHI) X2Y2 = T22**2 Z1 = (RN*(1.D0-ESQ)+HT)*DSIN(PHI) R = DSQRT(X2Y2+Z1**2) PSI = DATAN(Z1/DSQRT(X2Y2)) C C*******WGS84 VALUES ARE USED. GRAVN = GRAVA *(1.D0 + STAR * DSIN(PHI)**2) & / DSQRT(1.D0-ESQ*DSIN(PHI)**2) C C*******WGS72 VALUES ARE USED. C GRAVN = GRAVA *(1.D0 + STAR * DSIN(PHI)**2 C & + 0.000023461 * DSIN(PHI)**4) C GRAVN = GRAVN - HT*.3086D-5 DC = DCOS(PSI) TH = PI2-PSI Y = DSIN(TH) T = DCOS(TH) XX = T/Y F1 = AE/R F2 = (AE/R)**2 P11 = DSQRT(T3)*Y*F1**3 P11HT = DSQRT(T3)*Y*F2 D11 = - DSQRT(T3)*XX*F1**3 C 67 CONTINUE IF(TEMPL .EQ. RLAMD)GO TO 68 C RLAM = RLAMD * RPD SR(1) = DSIN(RLAM) CR(1) = DCOS(RLAM) DO 40 J = 2,NMAX SR(J)=T2*CR(1)*SR(J-1)-SR(J-2) CR(J)=T2*CR(1)*CR(J-1)-CR(J-2) 40 CONTINUE C 68 CONTINUE C DO 10 M = NMAX,0,-1 IF(TEMPP .EQ. PHID .AND. HGHT .EQ. HT)GO TO 11 C DO 9 N = NMAX,M,-1 LL = IV(N+1)+M+1 LL2= IV(N+2)+M+1 LL3= IV(N+3)+M+1 S11(N) = -A(LL2)*T*F1*S11(N+1)-B(LL3)*F2*S11(N+2)+CNM(LL) S12(N) = -A(LL2)*T*F1*S12(N+1)-B(LL3)*F2*S12(N+2)+SNM(LL) SG1(N) = -A(LL2)*T*F1*SG1(N+1)-B(LL3)*F2*SG1(N+2)+CNM(LL)*RNN(N) SG2(N) = -A(LL2)*T*F1*SG2(N+1)-B(LL3)*F2*SG2(N+2)+SNM(LL)*RNN(N) SP1(N) = -A(LL2)*T*F1*SP1(N+1)-B(LL3)*F2*SP1(N+2)-A(LL2)*F1* & S11(N+1) SP2(N) = -A(LL2)*T*F1*SP2(N+1)-B(LL3)*F2*SP2(N+2)-A(LL2)*F1* & S12(N+1) 9 CONTINUE C T11(M) = S11(M) T12(M) = S12(M) TG1(M) = SG1(M) TG2(M) = SG2(M) TP1(M) = SP1(M) TP2(M) = SP2(M) 11 CONTINUE C SV15(M)=-AS(M+1)*Y*F1*SV15(M+1)+T11(M)*CR(M)+T12(M)*SR(M) SV1(M)=-AS(M+1)*Y*F1* SV1(M+1)+TP1(M)*CR(M)+TP2(M)*SR(M) SV2(M)=-AS(M+1)*Y*F1* SV2(M+1)+ AS(M+1)*XX*F1*SV15(M+1) S2(M)=-AS(M+1)*Y*F1*S2(M+1)-T11(M)*M*SR(M)+T12(M)*M*CR(M) SG(M)=-AS(M+1)*Y*F1*SG(M+1)+TG1(M)*CR(M)+TG2(M)*SR(M) SHT(M)=-AS(M+1)*Y*F1*SHT(M+1)+T11(M)*CR(M)+T12(M)*SR(M) 10 CONTINUE C SUMX = (SV1(1)+SV2(1))*P11 + SV15(1)*D11 +(TP1(0)+TP2(0))*F2 SUMX = (-Y*SUMX*GAMMAA/GRAVN)*RHO SUME = S2(1)*P11 SUME = -(GAMMAA*SUME*RHO)/(GRAVN*DC) SUMG = (TG1(0)+TG2(0))*F2 + SG(1)*P11 SUMG = SUMG*GAMMAA*1.D5 SUMHT = (T11(0)+T12(0))*F1 + SHT(1)*P11HT SUMHT = SUMHT*GA/GRAVN DR = -(SUMG + 2.D0*SUMHT*GRAVN/R*1.D5) DNS = - GRAVN*SUMX/RHO*1.D5 DEW = - GRAVN*SUME/RHO*1.D5 DEFL = DSQRT(SUMX**2+SUME**2) TEMPP = PHID TEMPL = RLAMD HGHT = HT C WRITE (6,FMT=102) 102 FORMAT(4X,'PHI',4X,'LAMDA',5X,'HT',6X,'N',3X,'DELTA G', & 2X,'DIST R',2X,'DIST NS') WRITE (6,FMT=100) PHID,RLAMD,HT,SUMHT,SUMG,DR,DNS 100 FORMAT(1X,7(F7.2,1X)) WRITE (6,FMT=103) 103 FORMAT(26X,'DIST EW',4X,'XI',5X,'ETA',4X,'THETA') WRITE (6,FMT=101) DEW,SUMX,SUME,DEFL 101 FORMAT(25X,4(F7.2,1X)) WRITE (6,FMT=104) 104 FORMAT(//) WRITE(6,FMT=335) 335 FORMAT(1X,'ENTER 1 IF YOU HAVE ANOTHER COMPUTATION POINT') WRITE(6,FMT=336) 336 FORMAT(1X,'ENTER 0 IF YOU WANT TO STOP') READ (5,FMT=*) IQUIT IF(IQUIT.NE.0) GO TO 66 350 STOP END