DIMENSION XN(2000),XI(2000),ETA(2000),DG(2000),XLT(120),XLN(120), *XIGC(120),ETAGC(120),UNDGC(120),DGGC(120),GEOH(120,120),WRK(1257) C C PROGRAM TO COMPUTE GEOIDAL HEIGHTS, MERIDIAN AND PRIME VERTICAL C COMPONENTS OF THE DEVIATION OF THE VERTICAL AND FREE-AIR ANOMALIES C USING SPHERICAL HARMONIC COEFFICIENTS. C THIS PROGRAM HAS BEEN DEVELOPED BY G.LACHAPELLE AND MODIFIED BY C D. BLITZKOW - 1980. C REAL *8 MODEL COMMON/B/DEGLAT,DLAT,NLAT,DEGLON,DLON,NLON DATA MODEL/'GEM 10'/ C C DEFINITION OF THE CONSTANTS. C NDO - HIGHER ORDER OF THE COEFFICIENTS. C IPRINT - IDENTIFICATION THAT COULD BE NEGATIVE, ZERO OR POSITIVE: C NEGATIVE - COMPUTES EVERYTHING RELATED TO THE ELLIPSOID C TO WHICH THE GEM 10 MODEL IS REFERED. C ZERO - COMPUTES EVERYTHING FOR THE ELLIPSOID MENTIONED C ABOVE AND FOR THE LOCAL ELLIPSOID. C POSITIVE - COMPUTES EVERYTHING REFERED ONLY TO THE LOCAL C ELLIPSOID. C C DEGLAT,DEGLTB,DEGLON,DEGLNB - LATITUDES AND LONGITUDES LIMITS OF C THE AREA WHERE THE VALUES WILL BE C COMPUTED. C DLAT,DLON - INTERVAL OF THE COMPUTATION POINTS FOR LATITUDE AND C LONGITUDE RESPECTIVELY. C NDO = 30 IPRINT = 0 DEGLAT=1. DEGLTB=-25. DLAT=1.0 DEGLON=228. DEGLNB=320. DLON=1.0 C C WHEN THE VALUE OF LOWEST IS 1, VALUES OF THE LONGITUDES MEASURED C BY WEST ARE ALSO PRINTED; OTHERWISE IT PRINTS ONLY THE VALUES MEA- C SURED BY EAST. C LOWEST=1 NLAT=(DEGLAT-DEGLTB)/DLAT+1 NLON=(DEGLNB-DEGLON)/DLON+1 SCRD=3.24E05/ASIN(1.0) CALL NDEV77 (0,XLT,XLN,DGGC,XIGC,ETAGC,UNDGC,0,NDO) C C PARAMETERS OF NORMAL REFERENCE FIELD RELATED TO GEM 10 C A=6378140. F=1./298.257 C C PARAMETERS OF THE LOCAL REFERENCE FIELD C ABF=6378160. FBF=0.335292371299E-02 DABF=A-ABF DFBF=F-FBF C C TRANSLATION COMPONENTS OF THE ELLIPSOID AND GEOIDAL HEIGHT OF ZE- C RO ORDER. C TX=80.80 TY=14.81 TZ=44.01 XNO=0. J1=0 XLAT=DEGLAT*3600./SCRD DLATRD=DLAT*3600./SCRD DLONRD=DLON*3600./SCRD DO 100 ILAT=1,NLAT K1=NLAT+1-ILAT XLON=DEGLON*3600./SCRD DO 40 ILON=1,NLON XLT(ILON)=XLAT XLN(ILON)=XLON 40 XLON=XLON+DLONRD CALL NDEV77(1,XLT,XLN,DGGC,XIGC,ETAGC,UNDGC,NLON,NDO) NNA=NLON*(ILAT-1)+1 NNB=NNA+NLON-1 JJ=0 SINLAT=SIN(XLAT) DO 50 J=NNA,NNB K2=J-J1 JJ=JJ+1 XN(J)=UNDGC(JJ)+XNO+DABF-A*DFBF*SINLAT*SINLAT GEOH(K1,K2)=XN(J) C C THE GEOIDAL HEIGHTS HAVE BEEN PUT IN THE MATRIX 'GEOH' TO ALLOW C THE USE OF THE PROGRAM 'CONMAP' TO PLOT THE VALUES. C XI(J)=XIGC(JJ) ETA(J)=ETAGC(JJ) 50 DG(J)=DGGC(JJ) J1=J1+NLON 100 XLAT=XLAT-DLATRD C C CALL CONMAP C C PRINT RESULTS C DX=0.0 DY=0.0 DZ=0.0 IF(IPRINT)110,110,150 110 CONTINUE C CALL TABLE (XN,ABF,FBF,XNO,DX,DY,DZ,MODEL,1,LOWEST) CALL TABLE (XN,ABF,FBF,XNO,DX,DY,DZ,MODEL,1,LOWEST) C CALL TABLE (XI,ABF,FBF,XNO,DX,DY,DZ,MODEL,2,LOWEST) CALL TABLE (XI,ABF,FBF,XNO,DX,DY,DZ,MODEL,2,LOWEST) C CALL TABLE(ETA,ABF,FBF,XNO,DX,DY,DZ,MODEL,3,LOWEST) CALL TABLE(ETA,ABF,FBF,XNO,DX,DY,DZ,MODEL,3,LOWEST) C CALL TABLE (DG,ABF,FBF,XNO,DX,DY,DZ,MODEL,4,LOWEST) CALL TABLE (DG,ABF,FBF,XNO,DX,DY,DZ,MODEL,4,LOWEST) IF (IPRINT)200,150,150 150 CONTINUE DX=TX DY=TY DZ=TZ J1=0 XLAT=DEGLAT*3600./SCRD DO 180 ILAT=1,NLAT K1=NLAT+1-ILAT SINLAT=SIN(XLAT) COSLAT=COS(XLAT) XLON=DEGLON*3600./SCRD NNA=NLON*(ILAT-1)+1 NNB=NNA+NLON-1 DO 170 J=NNA,NNB K2=J-J1 SINLON=SIN(XLON) COSLON=COS(XLON) XN(J)=XN(J)-(-COSLAT*COSLON*DX-COSLAT*SINLON*DY-SINLAT*DZ-DA+A*DF* *SINLAT*SINLAT) C C THE GEOIDAL HEIGHTS HAVE BEEN PUT IN THE MATRIX 'GEOH' TO ALLOW C THE USE OF THE PROGRAM 'CONMAP' TO PLOT THE VALUES. C GEOH(K1,K2)=XN(J) XI(J)=XI(J)+(SINLAT*COSLON*DX+SINLAT*SINLON*DY-COSLAT*DZ+2.*A* *SINLAT*COSLAT*DF)*SCRD/A ETA(J)=ETA(J)+(SINLON*DX-COSLON*DY)*SCRD/A 170 XLON=XLON+DLONRD J1=J1+NLON 180 XLAT=XLAT-DLATRD CALL TABLE (XN,ABF,FBF,XNO,DX,DY,DZ,MODEL,1,LOWEST) C CALL TABLE (XI,ABF,FBF,XNO,DX,DY,DZ,MODEL,2,LOWEST) C CALL TABLE(ETA,ABF,FBF,XNO,DX,DY,DZ,MODEL,3,LOWEST) C C SPECIAL WIDTH AND HEIGHT FOR THE MAP ARE SELECTED DEFINING DY AND C DX. C C CALL CONMAP (GEOH,120,NLAT,NLON,1,1,2.,0,0,1.,1.,WRK,1257) C CALL ENDPLT 200 CONTINUE STOP END SUBROUTINE NDEV77(KKK,XLT,XLN,DGGC,XIGC,ETAGC,UNDGC,NRP,NDO) C C OCTOBER 77 - MODIFICATION OF SUBROUTINE NDEVG. C READ SPHERICAL HARMONIC COEFFICIENTS AND COMPUTE GEOIDAL HEIGHTS, C MERIDIAN AND PRIME VERTICAL COMPONENTS OF THE DEVIATION OF THE C VERTICAL AND FREE-AIR ANOMALIES. C REAL *8 MODEL DATA MODEL/'GEM 10'/ DIMENSION XJ(703),XK(703),P(37,37),PS(37) DIMENSION XLT(1),XLN(1),XIGC(1),ETAGC(1),UNDGC(1),DGGC(1) IF(KKK)1,1,3 C PARAMETERS OF GRS67. 1 A=6378160. F=0.335292371299E-02 EE=0.669460532856E-02 XJ2=0.10827E-02 RMEAN=6.371E06 GGG=9.798E05 SCRD=3.24E05/ASIN(1.0) C C READ THE SPHERICAL HARMONIC COEFFICIENTS FROM CARDS. C NCOEFF=(((NDO+1)*(NDO+2))/2)-3 DO 9 I=1,NCOEFF XJ(I)=0.0 9 XK(I)=XJ(I) C WRITE (6,6) MODEL C 6 FORMAT (/,9X,34HSPHERICAL HARMONIC COEFFICIENTS - ,A6,1X,5HMODEL, C 1,/,17X,20H(UNITS ARE IN 10**6),//) 8 FORMAT(2I3,1X,2F10.5) 5 READ(5,8,END=4)I,J,C,S IJK=(I*(I+1))/2+J-2 XJ(IJK)=C*10.E-07 XK(IJK)=S*10.E-07 C WRITE (6,7) I,J,C,S C 7 FORMAT (10X,2I4,2F20.5) GO TO 5 C C EVEN ZONAL COEFF. OF NORMAL FIELD (UP TO N=6 C 4 XJN2=-XJ2/SQRT(5.) XJN4=(EE*EE/35.)*(10.*XJ2/EE-1.) XJN6=-(EE*EE*EE/21.)*(15.*XJ2/EE-2.)/SQRT(13.) XJ(1)=XJ(1)-XJN2 XJ(8)=XJ(8)-XJN4 XJ(19)=XJ(19)-XJN6 RETURN C C CALCUALTE THE VALUES FOR NRP STATIONS. C 3 DO 100 K=1,NRP IF (K.EQ.1) GO TO 10 IF ((XLT(K)-XLT(K-1)).LT.1.0E-05) GO TO 15 10 XLAT=XLT(K) C XLATG GEOCENTRIC LATITUDE XLATG= ATAN((1.-EE)*TAN(XLAT)) SINLTG=SIN(XLATG) COSLTG=COS(XLATG) CALL SPHFT3(NDO,SINLTG,COSLTG,P,PS) C CALCULATE DSIN(XLATG)/D(XLAT) COSLAT=COS(XLAT) SINLAT=SIN(XLAT) DXLATG=(1.-EE)*COSLTG/(COSLAT*COSLAT+(1.-EE)*(1.-EE)*SINLAT 1*SINLAT) 15 XLONG=XLN(K) XISOM=0.0 ETASOM=0.0 UNDSOM=0.0 DGSOM=0.0 IJK=0 NNN=NDO+1 DO 50 N=3,NNN NA=N-2 DO 50 MA=1,N M=MA-1 IJK=IJK+1 XLONGM=XLONG*M NM=N-MA IF(NM)20,25,20 20 DP=P(MA,N) GO TO 30 25 DP=PS(N) 30 COSM=COS(XLONGM) SINM=SIN(XLONGM) XISOM=XISOM+(XJ(IJK)*COSM+XK(IJK)*SINM)*DP*DXLATG ETASOM=ETASOM+(-XJ(IJK)*M*SINM+XK(IJK)*M*COSM)*P(N,MA) UNDSOM=UNDSOM+(XJ(IJK)*COSM+XK(IJK)*SINM)*P(N,MA) 50 DGSOM=DGSOM+NA*(XJ(IJK)*COSM+XK(IJK)*SINM)*P(N,MA) XIGC(K)=-XISOM*SCRD ETAGC(K)=-ETASOM*SCRD/COSLAT UNDGC(K)=UNDSOM*RMEAN 100 DGGC(K)=DGSOM*GGG RETURN END SUBROUTINE SPHFT3(N,SINLTG,COSLTG,P,PS) C C THIS SUBROUTINE WAS PREPARED BY GERARD LACHAPELLE. C TECHNISCHE HOCHSCHUDE IN GRAZ, JUNE-JULY 1974 C THIS SUBROUTINE CALCULATE THE ASSOCIATED LEGENDRE FUNCTIONS AND THEI C FIRST DERIVATIVES UP TO AND INCLUDING DEGREE AND ORDER N TO BE SPE C CIFIED. THE DIMENSION OF P AND PS ARE P(N+1,N+1) AND PS(N+1) C SINLTG AND COSLTG ARE THE SINUS AND COSINUS OF THE GEOCENTRIC LATITU C THE ASSOCIATED LEGENDRE FUNCTIONS ARE STORED IN THE LOWER TRIANGULAR C PAR OF THE MATRIX P. THE FIRST DERIVATIVES OF THE ASSOCIATED LEGE C FUNCTIONS ARE STORED IN THE UPPER TRIANGULAR PART OF P (ZONAL AND C TESSERAL) AND IN THE VECTOR PS (SECTORIAL). C THE ASSOCIATED LEGENDRE FUNCTION OF DEGREE A AND ORDER B IS C THEN STORED IN P(A+1,B+1). ITS FIRST DEVIRATIVE IS STORED IN P (B+ C IF A IS DIFFERENT FROM B AND IN PS(A+1) IF A=B. C REAL P(37,37),PS(37) N1 = N+1 DO 20 I=1,N1 DO 20 J=1,N1 20 P(I,J)=0.0 P(1,1) = 1.0 P(1,2) = COSLTG*1.732050808 P(2,1)=SINLTG*1.732050808 P(2,2) = P(1,2) PS(1) = 0.0 PS(2) = -P(2,1) C C COMPUTATION OF SECTORIAL FUNCTIONS (N=M) AND THEIR FIRST DERIVATIVES C P2 = 2. COSP = COSLTG COSP1 = 1. DO 50 I = 3,N1 N2 = I - 1 P2 = P2*2. COSP = COSP*COSLTG T = (4.*N2+2.)*FACT(2*N2,N2-1)/FACT(N2-1,0) P(I,I) = SQRT(T/(N2*N2))*COSP/P2 PS(I) = -SQRT(T)*COSP1*SINLTG/P2 50 COSP1 = COSP1*COSLTG C C COMPUTATION OF ZONAL AND TESSERAL FUNCTIONS C DO 70 I = 3,N1 N2 = I - 1 DO 70 J = 1,N2 M = J - 1 70 P(I,J) = SQRT((2.*N2-1.)*(2*N2+1.)/((N2-M)*(N2+M)))*SINLTG*P(I-1,J 1)-SQRT((2.*N2+1.)*(N2+M-1.)*(N2-M-1)/((2.*N2-3)*(N2-M)*(N2+M)))*P( 2I-2,J) C C COMPUTATION OF THE FIRST DERIVATIVES OF THE ZONAL AND TESSERAL FUNCT C DO 90 I = 3,N1 N2 = I - 1 DO 90 J = 1,N2 M = J - 1 90 P(J,I) = (SQRT((N2-M)*(2.*N2+1.)/((N2+M)*(2.*N2-1.)))*(N2+M)*P(I- 11,J)-N2*SINLTG*P(I,J))/(COSLTG*COSLTG) RETURN END SUBROUTINE TABLE(X,A,F,XNO,DX,DY,DZ,MODEL,NDOVG,LOWEST) C C PRINT GEOIDAL HEIGHTS, MERIDIAN AND PRIME VERTICAL COMPONENTS AND C FREE-AIR ANOMALIES. C REAL X(1),LG(20),LGX(20) REAL *8 MODEL COMMON/B/DEGLAT,DLAT,NLAT,DEGLON,DLON,NLON 1 FORMAT(1H1,// ,60X,'GEOID UNDULATIONS',/) 2 FORMAT(1H1,// ,45X,'MERIDIAN COMPONENTS OF THE DEVIATION OF * THE VERTICAL',/) 3 FORMAT(1H1,// ,45X,'PRIME VERTICAL COMPONENTS OF THE DEVIAT *ION OF THE VERTICAL',/) 4 FORMAT(1H1,// ,55X,'FREE-AIR GRAVITY ANOMALIES',/) 10 FORMAT (20X,'A =',F10.1,5X,'F =',E12.5, 15X,'XNO =',F6.2,5X,'DX =',F6.1,5X,'DY =',F6.1,5X, *'DZ =',F6.1,5X,A6,//,8X,15F8.2) 11 FORMAT(/) 12 FORMAT(8X,15F8.2) 17 FORMAT(2X,F6.2,15F8.1) NCOL=15 NLINE=50 NTLON=(NLON+NCOL-1)/NCOL NTLAT=(NLAT+NLINE-1)/NLINE NLONZ=0 DO 60 I=1,NTLON DEGLOX=DEGLON+DLON*(I-1)*NCOL K=1 20 NLONZ=NLONZ+1 IF(DEGLOX.GT.360.)DEGLOX=DEGLOX-360. LG(K)=DEGLOX LGX(K)=360.-DEGLOX IF(NLONZ.EQ.NLON) GO TO 25 DEGLOX=DEGLOX+DLON IF(K.EQ.NCOL) GO TO 25 K=K+1 GO TO 20 25 DEGLAX=DEGLAT NLATX=1 30 M=1 J=NLATX/NLINE GO TO(31,32,33,34),NDOVG 31 WRITE(6,1) GO TO 36 32 WRITE(6,2) GO TO 36 33 WRITE(6,3) GO TO 36 34 WRITE(6,4) 36 WRITE(6,10)A,F,XNO,DX,DY,DZ,MODEL,(LG(L),L=1,K) IF(LOWEST.EQ.1)WRITE(6,12)(LGX(L),L=1,K) WRITE(6,11) 35 KKA=NLON*(J *NLINE+M-1)+(I-1)*NCOL+1 KKB=KKA+K-1 WRITE(6,17)DEGLAX,(X(MM),MM=KKA,KKB) DEGLAX=DEGLAX-DLAT IF(NLATX.EQ.NLAT) GO TO 60 NLATX=NLATX+1 IF(M.EQ.NLINE) GO TO 30 M=M+1 GO TO 35 60 CONTINUE RETURN END FUNCTION FACT (N,M) C CALCULATE FACT N/FACT M C RESTRICTION - N MUST BE GREATER OR EQUAL TO M IF (N.EQ.0.OR.N.EQ.1) GO TO 20 IF (N.EQ.M) GO TO 20 FACT=N N1=N K=N-M-1 DO 10 I=1,K N1=N1-1 10 FACT=FACT*N1 RETURN 20 FACT=1. RETURN END