C THIS PROGRAM IS DESIGNED FOR THE CALCULATION OF A GEOID UNDULATION C AT A POINT WHOSE LATITUDE AND LONGITUDE IS SPECIFIED. THE PROGRAM C IS DESIGNED TO USE THE POTENTIAL COEFFICIENT MODEL EGM96 AND A C SET OF SPHERICAL HARMONIC COEFFICIENTS OF A CORRECTION TERM. C THE CORRECTION TERM IS COMPOSED OF SEVERAL DIFFERENT COMPONENTS C THE PRIMARY ONE BEING THE CONVERSION OF A HEIGHT ANOMALY TO A GEOID C UNDULATION. THE PRINCIPLES OF THIS PROCEDURE WERE INITIALLY C DESCRIBED IN THE PAPER: USE OF POTENTIAL COEFFICIENT MODELS FOR GEOID C UNDULATION DETERMINATION USING A SPHERICAL HARMONIC REPRESENTATION C OF THE HEIGHT ANOMALY/GEOID UNDULATION DIFFERENCE BY R.H. RAPP, C JOURNAL OF GEODESY, 1996. C THIS PROGRAM IS DESIGNED TO BE USED WITH THE CONSTANTS OF EGM96 C AND THOSE OF THE WGS84(G873) SYSTEM. THE UNDULATION WILL REFER TO C THE WGS84 ELLIPSOID. C SPECIFIC DETAILS ON THE UNDULATION COMPUTATION WILL BE FOUND IN THE C JOINT PROJECT REPORT DESCRIBING THE DEVELOPMENT OF EGM96. C THIS PROGRAM IS A MODIFICATION OF THE PROGRAM DESCRIBED IN THE C FOLLOWING REPORT: C A FORTRAN PROGRAM FOR THE COMPUTATION OF GRAVIMETRIC QUANTITIES FROM C HIGH DEGREE SPHERICAL HARMONIC EXPANSIONS, RICHARD H. RAPP, C REPORT 334, DEPARTMENT OF GEODETIC SCIENCE AND SURVEYING, THE OHIO C STATE UNIVERSITY, COLUMBUS, 1982 C THIS PROGRAM WAS PUT IN THIS FORM IN DEC 1996. C RHRAPP.F477.NONLY C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C The input files consist of: C C correction coefficient set ("CORRCOEF") => UNIT = 1 C potential coefficient set ("EGM96") => UNIT = 12 C points at which to compute (INPUT.DAT") => UNIT = 14 C C The output file is: C C computed geoid heights ("OUTF477") => UNIT = 20 C C FILE ASSIGNMENT REVISIONS AT NIMA, DECEMBER 1996. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DIMENSIONS OF P,Q,HC,HS MUST BE AT LEAST ((MAXN+1)*(MAXN+2))/2, C DIMENSIONS OF SINML,COSML,SCRAP MUST BE AT LEAST MAXN, C WHERE MAXN IS MAXIMUM ORDER OF COMPUTATION C THE CURRENT DIMENSIONS ARE SET FOR A MAXIMUM DEGREE OF 360 IMPLICIT REAL*8 (A-H,O-Z) REAL*8 P(65341),SCRAP(361),RLEG(361),DLEG(361),RLNN(361), *SINML(361),COSML(361) REAL*8 HC(65341),HS(65341),CC(65341),CS(65341) DATA RAD/57.29577951308232D0/,HT/0.0D0/ C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCC CCCCCCCCCCC CCCCCCCCCC INPUT/OUPUT FILE NAMES CCCCCCCCCCC CCCCCCCCCC CCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CORRECTION COEFFICIENT FILE open(1,file='CORRCOEF',status='old') C C INPUT DATA FILE open(14,file='INPUT.DAT',status='old') C C POTENTIAL COEFFICIENT FILE open(12,file='EGM96',status='old') C C OUPUT FILE open(20,file='OUTF477.DAT',status='unknown') C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cccc CALL ERRSET(208,256,-1,2,IX,0) cccc010 READ(5,900) NMAX NMAX = 360 900 FORMAT(I3) WRITE(6,910) NMAX 910 FORMAT(1H1,//' MAXIMUM DEGREE =',I4//) L=65341 DO 71 I=1,L CC(I)=0.0D0 71 CS(I)=0.0D0 C THE CORRECTION COEFFICIENTS ARE NOW READ IN 72 READ(1,*,END=79)N,M,T1,T2 IG=(N*(N+1))/2+M+1 CC(IG)=T1 CS(IG)=T2 GO TO 72 79 CONTINUE C THE POTENTIAL COEFFICIENTS ARE NOW READ IN AND THE REFERENCE C EVEN DEGREE ZONAL HARMONIC COEFFICIENTS REMOVED TO DEGREE 6 CALL DHCSIN(NMAX,F,RJ2,RJ4,RJ6,HC,HS) C SETTING IFLAG=1 PREVENTS LEGENDRE FUNCTION DERIVATIVES BEING TAKEN C IN SUBROUTINE LEGFDN IFLAG=1 IR=0 K=NMAX+1 FLATL=91.0D0 C READ GEODETIC LATITUDE,LONGITUDE AT POINT UNDULATION IS WANTED cccc 030 READ(5,100,END=090) FLAT,FLON 030 READ(14,*,END=090) FLAT,FLON 100 FORMAT(6X,2F10.5) C COMPUTE THE GEOCENTRIC LATITUDE,GEOCENTRIC RADIUS,NORMAL GRAVITY CALL RADGRA(FLAT,FLON,HT,RLAT,GR,RE) IF(FLATL.EQ.FLAT) GO TO 040 RLAT1=RLAT RLAT=1.5707963267948966D0-RLAT FLATL=FLAT DO 25 J=1,K M=J-1 CALL LEGFDN(M,RLAT,RLEG,DLEG,NMAX,IR,RLNN,IFLAG) DO 26 I=J,K N=I-1 LOC=(N*(N+1))/2+M+1 26 P(LOC)=RLEG(I) 25 CONTINUE 040 RLON=FLON/RAD CALL DSCML (RLON,NMAX,SINML,COSML) CALL HUNDU(U,NMAX,P,HC,HS,SINML,COSML, *GR,RE,RLAT1,CC,CS,HACO) C U IS THE GEOID UNDULATION FROM THE EGM96 POTENTIAL COEFFICIENT MODEL C INCLUDING THE HEIGHT ANOMALY TO GEOID UNDULATION CORRECTION TERM C AND A CORRECTION TERM TO HAVE THE UNDULATIONS REFER TO THE C WGS84 ELLIPSOID. THE GEOID UNDULATION UNIT IS METERS. cccc WRITE(6,101)FLAT,FLON,U WRITE(20,101) FLAT,FLON,U 101 FORMAT(2F14.7,F10.3) GO TO 30 830 CONTINUE 90 STOP END SUBROUTINE HUNDU(UNDU,NMAX,P,HC,HS, *SINML,COSML,GR,RE,ANG,CC,CS,HACO) IMPLICIT REAL*8 (A-H,O-Z) cccc DIMENSION SINML(1),COSML(1),P(1) cccc REAL*8 HC(1),HS(1),CC(1),CS(1) DIMENSION SINML(361),COSML(361),P(65341) REAL*8 HC(65341),HS(65341),CC(65341),CS(65341) C CONSTANTS FOR WGS84(G873) DATA GM/.3986004418D15/,AE/6378137.0D0/ C GM IN UNITS OF M**3/S**2 AR=AE/RE ARN=AR AC=0.0 A=0.0 B=0.0 K=3 DO 030 N=2,NMAX ARN=ARN*AR K=K+1 SUM=P(K)*HC(K) SUMC=P(K)*CC(K) SUM2=0.0 DO 020 M=1,N K=K+1 TEMPC=CC(K)*COSML(M)+CS(K)*SINML(M) TEMP=HC(K)*COSML(M)+HS(K)*SINML(M) SUMC=SUMC+P(K)*TEMPC 020 SUM=SUM+ P(K) *TEMP AC=AC+SUMC 30 A=A+SUM*ARN AC=AC+CC(1)+P(2)*CC(2)+P(3)*(CC(3)*COSML(1)+CS(3)*SINML(1)) HACO=AC/100.D0 UNDU=A*GM/(GR*RE) C ADD HACO TO CONVERT HEIGHT ANOMALY ON THE ELLIPSOID TO THE UNDULATION C ADD -0.53M TO MAKE UNDULATION REFER TO THE WGS84 ELLIPSOID. UNDU=UNDU+HACO-0.53D0 RETURN END SUBROUTINE DSCML (RLON,NMAX,SINML,COSML) IMPLICIT REAL*8 (A-H,O-Z) cccc DIMENSION SINML(1),COSML(1) DIMENSION SINML(361),COSML(361) A=DSIN(RLON) B=DCOS(RLON) SINML(1)=A COSML(1)=B SINML(2)=2.0*B*A COSML(2)=2.0*B*B-1.0 DO 010 M=3,NMAX SINML(M)=2.0*B*SINML(M-1)-SINML(M-2) 010 COSML(M)=2.0*B*COSML(M-1)-COSML(M-2) RETURN END SUBROUTINE DHCSIN (NMAX,F,J2,J4,J6,HC,HS) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER IC REAL*8 J2,J4,J6,J8,J10,K REAL*8 HC(65341),HS(65341) C CONSTANTS FROM WGS84(G873) DATA GM,OMEGA,AE/.3986004418D20,7.292115D-5,6378137.0D0/ DATA RF/298.257223563D0/ C THE EVEN DEGREE ZONAL COEFFICIENTS GIVEN BELOW WERE COMPUTED FOR THE C WGS84(G873) SYSTEM OF CONSTANTS AND ARE IDENTICAL TO THOSE VALUES C USED IN THE NIMA GRIDDING PROCEDURE. COMPUTED USING SUBROUTINE C GRS WRITTEN BY N.K. PAVLIS J2=0.108262982131D-2 J4=-.237091120053D-05 J6=0.608346498882D-8 J8=-0.142681087920D-10 J10=0.121439275882D-13 M=((NMAX+1)*(NMAX+2))/2 DO 001 N=1,M HC(N)=0.0 001 HS(N)=0.0 912 FORMAT(2I4,2D30.20) cccc 02 READ(12,912,END=3)N,M,C,S 02 READ(12,*,END=3)N,M,C,S,EC,ES IF(N.GT.NMAX) GO TO 002 N=(N*(N+1))/2+M+1 HC(N)=C HS(N)=S GO TO 002 3 HC(4)=HC(4)+J2/DSQRT(5.D0) HC(11)=HC(11)+J4/3.0D0 HC(22)=HC(22)+J6/DSQRT(13.D0) HC(37)=HC(37)+J8/DSQRT(17.D0) HC(56)=HC(56)+J10/DSQRT(21.D0) RETURN END SUBROUTINE LEGFDN(M,THETA,RLEG,DLEG,NMX,IR,RLNN,IFLAG) C C THIS SUBROUTINE COMPUTES ALL NORMALIZED LEGENDRE FUNCTION C IN "RLEG" AND THEIR DERIVATIVES IN "DLEG". ORDER IS ALWAYS C M , AND COLATITUDE IS ALWAYS THETA (RADIANS). MAXIMUM DEG C IS NMX . ALL CALCULATIONS IN DOUBLE PRECISION. C IR MUST BE SET TO ZERO BEFORE THE FIRST CALL TO THIS SUB. C THE DIMENSIONS OF ARRAYS RLEG, DLEG, AND RLNN MUST BE C AT LEAST EQUAL TO NMX+1 . C C THIS PROGRAM DOES NOT COMPUTE DERIVATIVES AT THE POLES . C C IF IFLAG = 1 , ONLY THE LEGENDRE FUNCTIONS ARE C COMPUTED. C C ORIGINAL PROGRAMMER :OSCAR L. COLOMBO, DEPT. OF GEODETIC SCIENCE C THE OHIO STATE UNIVERSITY, AUGUST 1980 . ****************** C C IMPLICIT REAL*8 (A-H,O-Z) cccc DIMENSION RLEG(1),DLEG(1),RLNN(1) DIMENSION RLEG(361),DLEG(361),RLNN(361) 2, DRTS(1300),DIRT(1300) NMX1 = NMX+1 NMX2P = 2*NMX+1 M1 = M+1 M2 = M+2 M3 = M+3 IF(IR.EQ.1) GO TO 10 IR = 1 DO 5 N = 1,NMX2P DRTS(N) = DSQRT(N*1.D0) 5 DIRT(N) = 1.D0/DRTS(N) 10 COTHET = DCOS(THETA) SITHET = DSIN(THETA) IF(IFLAG.NE.1.AND.THETA.NE.0.D0)SITHI = 1.D0/SITHET C C COMPUTE THE LEGENDRE FUNCTIONS . C RLNN(1) = 1.D0 RLNN(2) = SITHET*DRTS(3) DO 15 N1 = 3,M1 N = N1-1 N2 = 2*N 15 RLNN(N1) = DRTS(N2+1)*DIRT(N2)*SITHET*RLNN(N1-1) IF(M.GT.1) GO TO 20 IF(M.EQ.0) GO TO 16 RLEG(2) = RLNN(2) RLEG(3) = DRTS(5)*COTHET*RLEG(2) GO TO 20 16 RLEG(1) = 1.D0 RLEG(2) = COTHET*DRTS(3) 20 CONTINUE RLEG(M1) = RLNN(M1) IF(M2.GT.NMX1) GO TO 35 RLEG(M2) = DRTS(M1*2+1)*COTHET*RLEG(M1) IF(M3.GT.NMX1) GO TO 35 DO 30 N1 = M3,NMX1 N = N1-1 IF(M.EQ.0.AND.N.LT.2.OR.M.EQ.1.AND.N.LT.3) GO TO 30 N2 = 2*N RLEG(N1) = DRTS(N2+1)*DIRT(N+M)*DIRT(N-M)*(DRTS(N2-1)*COTHET* 2 RLEG(N1-1)-DRTS(N+M-1)*DRTS(N-M-1)*DIRT(N2-3)*RLEG(N1-2)) GO TO 30 30 CONTINUE 35 IF(IFLAG.EQ.1) RETURN IF(SITHET.EQ.0.D0) WRITE(6,99) 99 FORMAT(//' *** LEGFDN DOES NOT COMPUTE DERIVATIVES AT THE POLES 2 *****************'//) IF(SITHET.EQ.0.D0) RETURN C C COMPUTE ALL THE DERIVATIVES OF THE LEGENDRE FUNCTIONS. C RLNN(1) = 0.D0 RLN = RLNN(2) RLNN(2) = DRTS(3)*COTHET DO 40 N1 = 3, M1 N = N1-1 N2 = 2*N RLN1 = RLNN(N1) RLNN(N1) = DRTS(N2+1)*DIRT(N2)*(SITHET*RLNN(N)+COTHET*RLN) RLN = RLN1 40 CONTINUE DLEG(M1) = RLNN(M1) IF(M2.GT.NMX1) RETURN DO 60 N1 = M2,NMX1 N = N1-1 N2 = N*2 DLEG(N1) = SITHI*( N *RLEG(N1)*COTHET-DRTS(N-M)*DRTS(N+M)* 2 DRTS(N2+1)*DIRT(N2-1)*RLEG(N)) 60 CONTINUE RETURN END SUBROUTINE RADGRA(FLAT,FLON,HT,RLAT,GR,RE) IMPLICIT REAL*8(A-H,O-Z) REAL*8 N,K C THIS SUBROUTINE COMPUTES GEOCENTRIC DISTANCE TO THE POINT, C THE GEOCENTRIC LATITUDE,AND C AN APPROXIMATE VALUE OF NORMAL GRAVITY AT THE POINT BASED C THE CONSTANTS OF THE WGS84(G873) SYSTEM ARE USED DATA AE/6378137.D0/,E2/.00669437999013D0/,RAD/57.29577951308232D0/ DATA GEQT,K/9.7803253359D0,.00193185265246D0/ FLATR=FLAT/RAD FLONR=FLON/RAD T1=DSIN(FLATR)**2 N=AE/DSQRT(1.-E2*T1) T2=(N+HT)*DCOS(FLATR) X=T2*DCOS(FLONR) Y=T2*DSIN(FLONR) Z=(N*(1.-E2)+HT)*DSIN(FLATR) N=AE/DSQRT(1.-E2*T1) C COMPUTE THE GEOCENTRIC RADIUS RE=DSQRT(X**2+Y**2+Z**2) C COMPUTE THE GEOCENTRIC LATITUDE RLAT=DATAN(Z/DSQRT(X**2+Y**2)) C COMPUTE NORMAL GRAVITY:UNITS ARE M/SEC**2 GR=GEQT*(1.+K*T1)/DSQRT(1.-E2*T1) RETURN END