C GENERATING A "GEOID" FOR A SMALL AREA ON THE SURFACE OF THE EARTH. C AREA,A,IS A SPHERICAL CAP OF RADIUS (PSI0) AND POINT,O,THE CENTRE. C THE GEOIDAL UNDULATION,IS EXPRESSED AS A SERIES OF SPHERICAL CAP C FUNCTIONS R AND S ,ORTHOGONAL IN THE AREA: C NM NM C L N C N=ZIG ZIG(A R +B S ). WHERE C N=0 M=0 NM NM NM NM C _ C R =COS(M*AZ)*P (COS(K*PSI)), C NM NM C _ C S =SIN(M*AZ)*P (COS(K*PSI)), C NM NM C WHERE K IS THE SCALE K=PI/PSI0. C C R AND S,ARE ORTHOGONAL WITH A WEIGHT FUNCTION ,W ,DEFINED AS C C W(PSI0,PSI)=K*SIN(K*PSI)/SIN(PSI). C C AZ,PSI- ARE POLAR COORDINATES OF A POINT P IN THE AREA. C C THE POLAR SYSTEM IS A LOCAL SYSTEM CENTERD AT POINT O (CENTRE OF C THE AREA) WHOSE MERIDIAN MEASURES THE ZERO AZIMUTH (AZ),PSI IS THE C SPHERICAL DISTANCE BETWEEN THE POINT AND THE CENTER. C C A ,B - ARE THE COEFFICIENT TO BE DETEMINED. THEY ARE COMPUTED C NM NM C BY SURFACE INTEGRALS RESEMBLING THE STOKES' INTEGRAL WITH THE C FOLLOWING KERNELS OF INTEGRATION: C C C K =COS(M*AZ)*F (PSI0,PSI), C NM NM C S C K =SIN(M*AZ)*F (PSI0,PSI). WHERE C NM NM C C F =K*INTERAL OF (F1*F2*F3) C NM C WHERE C _ INF _ _ C F1=SIN(K*PSI) ,F2=P (K*PSIS) ,F3=ZIG(P (PSI)*P (PSIS)/(I-1)/C). C NM I=M IM IM C WHERE C=2 FOR M=0 AND C=4 FOR M>0. C NOTE: FOR FURTHER INFORMATIONS,SEE 'LOCAL EVALUATION OF THE GEOID C FROM GRAVITY DATA' BY MEHDI NAJAFI,A THESIS SUBMITTED IN C PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF C MASTER OF SCIENCE IN ENGINEERING IN THE DEPARTEMENT OF C SURVEYING ENGINEERING IN UNB. C THE MAIN PROGRAM CONSISTS OF THE FOLLOWING PARTS: C PART I- COMPUTES THE KERNELS F (PSI0,PSI) FOR N,M=0...NMAX, C NM C AND FOR PSI FROM ZERO UP TO "PSIMAX" IN STEPS C OF 10 MINUTES AND 1 DEGREE. C THEY ARE STORED IN AN ARRAY CALLED QNM. C PART II- COMPUTES THE GEOIDAL HEIGHT OF POINT O USING THE STOKES' C INTEGRAL AND THE COEFFICIENTS USING THE SURFACE INTEGRALS. C ALL THE INTEGRATIONS ARE EXTENDED OVER A RECTANGLE OF C DIMENSION "LATOUT".GRAVITY DATA ARE SUPLIED BY TWO FILES; C FILE #1(5X5 MINUTE) TO BE USED IN INNER ZONE "LATINN" AND C FILE #2(1X1 DEGREE) TO BE USED IN OUTER ZONE. C PART III- COMPUTES THE GEOID FOR A GIVEN AREA ,AND C PRINTS OUT THE COEFFICIENTS AND THE GEOIDAL HEIGHTS. C C INPUT DATA- C PART I- NMAX- IS THE MAXIMUM VALUE FOR N. C NCAP- MAXIMUM DEGREES UP TO WHICH THE ASSOCIATED C LEGENDRE FUNCTIONS ARE COMPUTED. C PSI0- SPHERICAL RADIUS OF THE CAP(AREA). C INFINT- THE MAXIMUM LIMIT OF THE SUMMATION IN F3. C MPOINT- THE NUMBER OF BASE POINTS USED BY THE C GAUSS LEGENDRE QUADRATURE.IT CAN TAKE ON 2,3,4,5,6 C 10 OR 15. NOTE: THE INPUT VARIABLES ARE ASSIGNED C IN THE PROGRAM. C PART II- LATA- LATITUDE OF POINT O IN DEGREES. C LONA- LONGITUDE OF POINT O IN DEGREES.BOTH LATA AND C LONA SHOULD BE GIVEN AS INTEGER DEGREES. C THESE ARE READ FROM UNFORMATTED CARD NO 1. C T=R/4*PI*G- A CONSTANT IN METERS PER MILLIGALS C WHERE R IS THE MEAN RADIUS OF THE EARTH AND G IS THE C MEAN GRAVITY OVER THE EARTH. C HLAT20- DIMENSION OF INNERMOST ZONE (MINUTE) C LATINN- DIMENSION OF INNER ZONE (IN DEGREES). C LATOUT- DIMENSION OF OUTER ZONE (IN DEGREES). C FILE #1- A TAPE OF (5X5) MINUTES FREE-AIR GRAVITY C ANOMALIES GIVEN IN MILLIGALS,LOGICAL UNIT NO 55. C FILE #2- A TAPE OF (1X1) DEGREE FREE-AIR GRAVITY C ANOMALIES GIVEN IN MILLIGALS,LOGICAL UNIT NO 12. C PARTIII- SMLF,BIGF,SMLL,BIGL- ARE TWO PARALLELS AND TWO C MERIDIANS ENCLOSING THE AREA OF INTEREST. C DF,DL- ARE INCREMENTS IN LATITUDE AND LONGITUDE. C THES ARE READ FROM UNFORMATED CARD NO 2. C C OUTPUT DATA- C PARTIII- THE COEFFICIENTS AND THE GEOID. C C M. NAJAFI 9/1980 IMPLICIT REAL*8(A-H,O-Z) COMMON/CF/ COF COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CQ/ QNM(121,45),QSINT1,QSINT2,PSI00,IQ1,IQ2,IQNM,JQNM COMMON/CZ/ ZPP(20,16),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M COMMON/SKS/ STKA(4500),STKB(2000),STKC(900),PSIA,PSIB,PSIC, &HLAT20,HLON20,DG,DGIN,DGINA, &HLTA,HLNA,HLTB,HLNB,HLTC,HLNC,INTSTA,INTSTB,INTSTC, &IDGXYA,IDGXYB,IDGXYC,IDGXYD,NINP,NINPA REAL*8 PN(16,16),FNM(20,20) EXTERNAL F C REAL*8 ZINMST(100)/100*0.D0/,ZINNER(100)/100*0.D0/,ZOUTER(100)/ &100*0.D0/ REAL*8 AVKLS(100),ABNM(100) IPN=16 IPNM=541 PI=DARCOS(-1.D0) DR=PI/180.D0 C INPUT DATA **************(PART I)************************* C NCAP=10 C INFINT=180 C C NMAX IS THE MAXIMUM ORDER UP TO WHICH THE KERNELS ARE COMPUTED. C IT CAN TAKE ON ANY VALUE FROM 0 TO NCAP-1. C NMAX=5 C MPOINT=15 C PSI0=10.D0 C C THE VARIABLE (COF) MUST BE KEPT EQUAL TO UNITY UNLESS WE WANT C TO COMPUTE THE GEOID FOR A DIFFERENT ,THAN THE SELECTED CAP, C SIZE CAP. C COF=DSQRT(0.5D0) COF=1.D0 C NMP1=NMAX+1 C C INPUT DATA **********( PART II)**************************** C C GEODETIC COORDINATES OF POINT O (THE CENTRE),AND THE FOUR C BOUNDARIES(TWO PARALLELS AND TWO MERIDIANS) OF THE AREA. C READ(5,*) LATA,LONA READ(5,*) SMLF,SMLL,BIGF,BIGL,DF,DL C C DIMENSIONS OF THE INNER AND OUTER RECTANGLE. C LATINN=8 LATOUT=20 C C DIMENSION OF THE INNERMOST ZONE. C HLAT20=10.D0/60.D0*DR HLAT20=HLAT20/2 HLON20=HLAT20 C WRITE(6,64) WRITE(6,65) LATA WRITE(6,68) LONA WRITE(6,66) WRITE(6,67) BIGF,BIGL WRITE(6,69) SMLF,SMLL WRITE(6,63) DF,DL 64 FORMAT('1',31X,'DEGREES') 63 FORMAT(/////,6X,'INTERVAL IN LATITUDES',5X,F5.2,2X,'DEG.S'/ ' & INTERVAL IN LONGITUDS',5X,F5.2,2X,'DEG.S') 65 FORMAT('0',5X,'LATITUDE OF THE CENTRE',5X,I5) 68 FORMAT('0',5X,'LONGITUD OF THE CENTRE',5X,I5) 66 FORMAT(/////,26X,'LATITUDE (DEG.S)',10X,'LONGITUDE (DEG.S)') 67 FORMAT('0',' NORTH EAST CORNER',12X,F6.2,16X,F7.2) 69 FORMAT('0',' SOUTH WEST CORNER',12X,F6.2,16X,F7.2) C SCRD=206264.80625 RDDG=3600/SCRD TT=6371000/1.2312529D+7 20 FORMAT(F11.0,2F12.7,10F10.1) 22 FORMAT(12X,I3,1X,I4,11X,I5,15X,F6.1,14X,F6.1,11X,F6.1) C 25 FORMAT(I5,1X,I5,2X,I4,1X,I2) XLONA=LONA*DR XLATA=LATA*DR SINLA=DSIN(XLATA) COSLA=DCOS(XLATA) C C (A PATTERN TO DIVIDE THE AREA AROUND COMPUTATION POINT,THE CENTRE) C C THE THREE ARRAYS:STKA(NSTA),STKB(NSTB) AND STKC(NSTC) CONTAIN THE C STOKES FUNCTION VALUES;THE STOKES FUNCTION IS EVALUATED IN THREE C INTERVALS:INTSTA,INTSTB AND INTSTC SECONDS. INTSTA=3 INTSTB=30 INTSTC=120 C NSTA=4500 NSTB=2000 NSTC=900 C PSIA=NSTA*INTSTA PSIB=NSTB*INTSTB PSIC=NSTC*INTSTC C C THE THREE NESTED RECTANGLES OF DIMENSIONS:NLATA,NLATB AND NLATC C IN DEGREES,INTO WHICH THE DATA BLOCKS ARE SUBDIVIDED INTO SMALLER C BLOCKS OF DIMENSION:IDGXYA,IDGXYB AND IDGXYC IN SECONDS.DATA C DATA BLOCKS BEYOND THE NLATC RECTANGLE ARE SUBDIVIDED INTO IDGXYD C SIZE BLOCKS. NLATA=4 NLATB=8 NLATC=10 C IDGXYA=60 IDGXYB=300 IDGXYC=1800 IDGXYD=3600 C C EXPANDING THE DIMENSIONS OF THE RECTANGLES,REGARDING THE LATITUDE C OF THE CENTRE. C IF(LATA.GE.50) GO TO 80 NLONA=NLATA NLONB=NLATB NLONC=NLATC NLOND=LATINN NLONE=LATOUT GO TO 83 80 CONTINUE IF(LATA.GE.70) GO TO 81 NLONA=NLATA*2 NLONB=NLATB*2 NLONC=NLATC*2 NLOND=LATINN*2 NLONE=LATOUT*2 HLON20=HLAT20*2 GO TO 83 81 CONTINUE IF(LATA.GE.75) GO TO 82 NLONA=NLATA*3 NLONB=NLATB*3 NLONC=NLATC*3 NLOND=LATINN*3 NLONE=LATOUT*3 HLON20=HLAT20*3 GO TO 83 82 CONTINUE NLONA=NLATA*4-1 NLONB=NLATB*4-1 NLONC=NLATC*4-1 NLOND=LATINN*4-1 NLONE=LATOUT*4-1 83 CONTINUE C I=NLATA/2.D0 HLTA=I*1.D0-0.5D0 I=NLONA/2.D0 HLNA=I*1.D0-0.5D0 I=NLATB/2.D0 HLTB=I*1.D0-0.5D0 I=NLONB/2.D0 HLNB=I*1.D0-0.5D0 I=NLATC/2.D0 HLTC=I*1.D0-0.5D0 I=NLONC/2.D0 HLNC=I*1.D0-0.5D0 C C SIDES OF THE INNER AND OUTER ZONES. C I=LATA-LATINN/2.0 J=LATA+LATINN/2.0 K=LONA-NLOND/2.0 L=LONA+NLOND/2.0 A1=I+0.5 B1=J-0.5 C1=K+0.5 D1=L-0.5 I=LATA-LATOUT/2.0 J=LATA+LATOUT/2.0 K=LONA-NLONE/2.0 L=LONA+NLONE/2.0 A2=I+0.5 B2=J-0.5 C2=K+0.5 D2=L-0.5 C***********************(PART I )******************************* C C CONSTRUCTING A TABLE OF NORMALIZED ASSOCIATED LEGENDRE FUNCTIONS. C NCAP1=NCAP+1 SDIS=0.D0 PSINT=10.D0/60 DO 10 IR=1,IPNM PHI=90.D0-SDIS CALL LEGDRE(1,PHI,PN,IPN,NCAP1) K=0 DO 11 I=1,NCAP1 DO 11 J=1,I K=K+1 11 PNM(IR,K)=PN(I,J) 10 SDIS=SDIS+PSINT JPNM=K C C COMPUTING THE KERNELS FOR A PSI FROM 0.D0 TO PSIMAX DEGREES. C PSIMAX=DARCOS(DSIN(A2*DR)*SINLA+DCOS(A2*DR)*COSLA*DCOS(XLONA- &C2*DR)) PSIMAX=PSIMAX/DR+1.D0 PSI=0.D0 QSINT1=10.D0/60 QSINT2=1.D0 IQ1=PSI0/QSINT1+1.5D0 PSI00=(IQ1-1)*QSINT1 IQ2=(PSIMAX-PSI00)/QSINT2+0.5D0 IQNM=IQ1+IQ2 DO 30 I=1,IQNM CALL ZGMAPP(MPOINT,NMAX,0.D0,PSI0*DR) AA=0.D0 BB=PSI0*DR CALL ISOKNL(AA,BB,NMAX,MPOINT,FNM,II,JJ) J=0 DO 32 IR=1,NMP1 DO 32 IC=1,IR J=J+1 32 QNM(I,J)=FNM(IR,IC) QSINT=QSINT1 IF(I.GE.IQ1) QSINT=QSINT2 30 PSI=PSI+QSINT JQNM=J C C*****************************( PART II)**************************** C C A PRIORI COMPUTATION OF THE STOKES FUNCTIONS. C DPSI=INTSTA/SCRD PSI=0.D0 DO 100 I=1,NSTA PSI=PSI+DPSI 100 STKA(I)=SF(PSI) DPSI=INTSTB/SCRD PSI=0.D0 DO 110 I=1,NSTB PSI=PSI+DPSI 110 STKB(I)=SF(PSI) DPSI=INTSTC/SCRD PSI=0.D0 DO 120 I=1,NSTC PSI=PSI+DPSI 120 STKC(I)=SF(PSI) C C FILE #1 - MEAN GRAVITY ANOMALIES OF (5X5) MINUTE. C IN1=0 DGIN=0.D0 DGINA=0.D0 NINP=0 NINPA=0 MLAT=5 MLON=5 SOMA=0.D0 REWIND 55 1 READ(55,20,END=5) R1,XF,XL,DG,R2,R3,R4,R5,R6,R7,R8,R9,R10 IF(XF.LT.A1.OR.XF.GT.B1.OR.XL.LT.C1.OR.XL.GT.D1) GO TO 1 C C COMPUTING THE INNER ZONE (A1,B1,C1,D1). C IN1=IN1+1 MLAT=5 MLON=5 IF(XF.GE.50.D0) MLON=10 IF(XF.GE.70.D0) MLON=30 XSIGMA=(MLAT*60.0/SCRD)*(MLON*60.0/SCRD) DSIGMA=XSIGMA*DCOS(XF*DR) C CALL STKSM(SINLA,COSLA,XLATA,XLONA,XF,XL,MLAT,MLON,S) SOMA1=S*DG*DSIGMA SOMA=SOMA+SOMA1 C CALL CSKNM(SINLA,COSLA,XLATA,XLONA,XF,XL,MLAT,MLON,NMAX, & AVKLS,NK) DO 40 KK=1,NK 40 ZINNER(KK)=ZINNER(KK)+AVKLS(KK)*DG*DSIGMA GO TO 1 5 CONTINUE C C CONTRIBUTION OF INNER ZONE. C HINNER=SOMA*TT DO 41 I=1,NK 41 ZINNER(I)=ZINNER(I)*TT C C CONTRIBUTION OF INNER MOST ZONE. C DGIN=DGIN/NINP HINMST=8*PI*TT*DGIN*DSQRT(HLAT20*HLON20*COSLA/PI) C DGINA=DGINA/NINPA SR2=4*HLAT20*HLON20*COSLA/PI K=1 ZINMST(K)=Q(0,0,0.D0)*PI*SR2*DGINA*TT IF(NMAX.EQ.0) GO TO 200 DO 31 N=1,NMAX K=K+1 ZINMST(K)=Q(N,0,0.D0)*PI*SR2*DGINA*TT DO 33 M=1,N K=K+1 ZINMST(K)=0.D0 K=K+1 33 ZINMST(K)=0.D0 31 CONTINUE 200 CONTINUE IN1=IN1-NINP C C FILE #2 - MEAN GRAVITY ANOMALIES (1X1) DEGREE FOR CANADA. C REWIND 12 IN2=0 MLAT=60 MLON=60 SOMA=0.D0 XSIGMA=(MLAT/60.D0*DR)*(MLON/60.D0*DR) 2 READ(12,22,END=6) LATT,LONN,IREC,FMIN,FMAX,DG XF=LATT XL=360.D0+LONN IF(XF.GT.A1.AND.XF.LT.B1.AND.XL.GT.C1.AND.XL.LT.D1) GO TO 2 IF(XF.LT.A2.OR.XF.GT.B2.OR.XL.LT.C2.OR.XL.GT.D2) GO TO 2 C C COMPUTING OUTER ZONE (A2-A1,B2-B1,C2-C1,D2-D1). C IN2=IN2+1 DSIGMA=XSIGMA*DCOS(XF*DR) C CALL STKSM(SINLA,COSLA,XLATA,XLONA,XF,XL,MLAT,MLON,S) SOMA1=S*DG*DSIGMA SOMA=SOMA+SOMA1 C CALL CSKNM(SINLA,COSLA,XLATA,XLONA,XF,XL,MLAT,MLON,NMAX, & AVKLS,NK) DO 60 KK=1,NK 60 ZOUTER(KK)=ZOUTER(KK)+AVKLS(KK)*DG*DSIGMA GO TO 2 6 CONTINUE C C CONTRIBUTION OF THE OUTER ZONE. C HOUTER=SOMA*TT DO 61 I=1,NK 61 ZOUTER(I)=ZOUTER(I)*TT C C FINAL COMPUTATION OF THE GEOID FOR THE CENTRE AND THE COEFFICIENTS. C TOTALN=HINMST+HINNER+HOUTER DO 90 I=1,NK 90 ABNM(I)=ZINMST(I)+ZINNER(I)+ZOUTER(I) C WRITE(6,112) HINMST,NINP WRITE(6,113) HINNER,IN1 WRITE(6,114) HOUTER,IN2 WRITE(6,111) LATA,LONA,TOTALN 111 FORMAT('0',5X,'LAT=',I5,5X,'LON=',I6,5X,'N=',F8.2) 112 FORMAT('0','INNER MOST CONTRIBUTION',F8.2,5X,'INNER MOST COUNTS=', &I8) 113 FORMAT('0','INNER ZONE " ',F8.2,5X,'INNER ZONE COUNTS=', &I8) 114 FORMAT('0','OUTER ZONE " ',F8.2,5X,'OUTER ZONE COUNTS=', &I8) C WRITE(6,115) DO 91 I=1,NK 91 WRITE(6,116) I,ABNM(I),ZINMST(I),ZINNER(I),ZOUTER(I) 115 FORMAT('1',5X,'COEFFICIENTS',5X,'INNER MOST',5X,'INNER ZONE',5X,'O &UTER ZONE') 116 FORMAT('0',I5,2X,F8.2,9X,F6.2,8X,F8.2,7X,F8.2) C C**************************(PART III)************************** C CALL REGOID(BIGF,SMLF,BIGL,SMLL,DF,DL,ABNM,NK,TOTALN,SINLA,COSLA, &XLATA,XLONA,NMAX,PSI0) RETURN END SUBROUTINE ISOKNL(A,B,INMAX,MPOINT,FNM,IROW,ICUL) C C ISOKNL COMPUTES THE F (PSI0,PSI) FUNCTIONS (SEE THE MAIN PROGRAM) C NM C OF N,M=0...INMAX FOR A GIVEN SPHERICAL DISTANCE PSI. C INPUTE DATA: A,B- LOWER AND UPPER LIMIT OF THE INTEGRATION C IN RADIANS. C INMAX- THE MAXIMUM LIMIT OF N. C MPOINT- NUMBER OF BASE POINTS USED BY ROUTINE GAUSS. C OUTPUT DATA: FNM(IROW,ICUL)- THE ARRAY CONTAINING THE FUNCTION C VALUES. C M. NAJAFI, 10/1980 IMPLICIT REAL*8(A-H,O-Z) REAL*8 FNM(20,20) EXTERNAL F COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CZ/ ZPP(20,16),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M SCALE=180.D0/PSI0 PI=DARCOS(-1.D0) DR=PI/180.D0 DO 150 I=1,20 DO 150 J=1,20 150 FNM(I,J)=0.D0 C.....PSI AND PSI0 HAVE TO BE GIVEN IN DEGREES. NBP1=INMAX+1 DO 100 IROW=1,NBP1 N=IROW-1 DO 200 ICUL=1,IROW M=ICUL-1 FNM(IROW,ICUL)=GAUSS(A,B,MPOINT,F)*SCALE 200 CONTINUE 100 CONTINUE RETURN END SUBROUTINE ZGMAPP(MP,MMAX,A,B) C C ZGMAPP COMPUTES THE F3 FUNCTION (SEE MAIN PROGRAM) OF M=0...MMAX C FOR A GIVEN PSI AND FOR THE SAME BASE POINTS USED BY THE C ROUTINE GAUSS. C INPUT DATA: MP- NUMBER OF THE BASE POINTS WHICH INTRODUCES C THE BASE POINTS TO THE ROUTINE. C MMAX- THE MAXIMUM ORDER (M). C A,B- LOWER AND UPPER LIMIT OF THE INTEGRATION. C OUTPUT DATA: ZPP(IZP,JZP)- AN ARRAY CONTAINING THE VALUES C OF F3 FUNCTION. C M. NAJAFI 10/1980 IMPLICIT REAL*8 (A-H,O-Z) REAL*8 ZI(24),XI(41) INTEGER NPOINT(7),KEY(8) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CZ/ ZPP(20,16),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M C C.....PRESET NPOINT,KEY,ZI,AND WEIGHT ARRAYS..... C DATA NPOINT/2,3,4,5,6,10,15/ C DATA KEY/1,2,4,6,9,12,17,25/ C DATA ZI/0.577350269D0,0.0D0 ,0.774596669D0,0.339981044D0, & 0.861136312D0,0.0D0 ,0.538469310D0,0.906179846D0, & 0.238619186D0,0.661209387D0,0.932469514D0,0.148874339D0, & 0.433395394D0,0.679409568D0,0.865063367D0,0.973906529D0, &0.0D0,0.201194093997D0,0.394151347078D0,0.570972172609D0,0.72441 &7731360D0,0.848206583410D0,0.937273392401D0,0.987992518020D0/ PI=DARCOS(-1.D0) DR=PI/180.D0 C C.....FIND SUBSCRIPT OF FIRST ZI AND WEIGHT VALUE...... C DO 1 I=1,7 IF(MP.EQ.NPOINT(I)) GO TO 2 1 CONTINUE C C.....INVALID M USED. C WRITE(6,11) 11 FORMAT(' ','******INVALID M HAS BEEN USED.******') GAUSS=0.D0 RETURN C C.....SET UP INITIAL PARAMETERS...... C 2 JFIRST=KEY(I) JLAST=KEY(I+1)-1 C=(B-A)/2.D0 D=(B+A)/2.D0 INFT=INFINT JZP=MMAX+1 K=0 DO 5 J=JFIRST,JLAST K=K+1 IF(ZI(J).NE.0.D0) GO TO 6 XI(K)=D GO TO 5 6 CONTINUE XI(K)=D+ZI(J)*C K=K+1 XI(K)=D-ZI(J)*C 5 CONTINUE IZP=K DO 7 II=1,IZP PSISTR=XI(II) DO 8 JJ=1,JZP M=JJ-1 FOUR=4.D0 IF(M.EQ.0) FOUR=2.D0 I1=M IF(M.LT.2) I1=2 Z=0.D0 P12=0.D0 P22=0.D0 DO 10 I=I1,NCAP P11=P12 P21=P22 P12=P(I,M,PSI) P22=P(I,M,PSISTR/DR) CALL TRAPS(0,0,1000,0,0) Z1=P12*P22/(I-1)/FOUR 10 Z=Z+Z1 KKK=NCAP+1 COSDD=DCOS(PSI*DR) COSTR=DCOS(PSISTR) DO 30 I=KKK,INFT CO1=(2*I-1)*1.D0/(I-M)*H(I,M)/H(I-1,M) CO2=(I+M-1)*1.D0/(2*I-1)*H(I-1,M)/H(I-2,M) HOLD=P12 P12=CO1*(COSDD*P12-CO2*P11) P11=HOLD HOLD=P22 P22=CO1*(COSTR*P22-CO2*P21) P21=HOLD Z1=P12*P22/(I-1)/FOUR 30 Z=Z+Z1 ZPP(II,JJ)=Z 8 CONTINUE 7 CONTINUE RETURN END DOUBLE PRECISION FUNCTION GAUSS(A,B,MP,F) C C GAUSS-LEGENDRE QUADRATURE EVALUATES AN INTEGRAL BETWEEN THE C INTERVAL (A,B),SEE CARNAHAN 1976. C INPUT DATA: A,B- LOWER AND UPPER LIMIT OF THE INTEGRATION. C MP- DEGREE OF LEGENDRE POLYNOMIAL WHOSE ZEROES C ARE USED AS THE BASE POINTS IN THE ROUTINE. C F- VALUE OF THE INTEGRAND AT A POINT. C OUTPUT DATA: GAUSS- VALUE OF THE INTEGRAL. C M. NAJAFI, 10/1980 IMPLICIT REAL*8 (A-H,O-Z) REAL*8 ZI(24),WEIGHT(24) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CZ/ ZPP(20,16),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M INTEGER NPOINT(7),KEY(8) C C.....PRESET NPOINT,KEY,ZI,AND WEIGHT ARRAYS..... C DATA NPOINT/2,3,4,5,6,10,15/ C DATA KEY/1,2,4,6,9,12,17,25/ C DATA ZI/0.577350269D0,0.0D0 ,0.774596669D0,0.339981044D0, & 0.861136312D0,0.0D0 ,0.538469310D0,0.906179846D0, & 0.238619186D0,0.661209387D0,0.932469514D0,0.148874339D0, & 0.433395394D0,0.679409568D0,0.865063367D0,0.973906529D0, &0.0D0,0.201194093997D0,0.394151347078D0,0.570972172609D0,0.72441 &7731360D0,0.848206583410D0,0.937273392401D0,0.987992518020D0/ C DATA WEIGHT/1.D0 ,0.888888889D0,0.555555556D0,0.652145155D0, & 0.347854845D0,0.568888889D0,0.478628671D0,0.236926885D0, & 0.467913935D0,0.360761573D0,0.171324492D0,0.295524225D0, &0.269266719D0,0.219086363D0,0.149451349D0,0.066671344D0,0.2025782 &41926D0,0.198431485327D0,0.186161000116D0,0.166269205817D0,0.1395 &70677926D0,0.107159220467D0,0.070366047488D0,0.030753241996D0/ C C.....FIND SUBSCRIPT OF FIRST ZI AND WEIGHT VALUE...... C DO 1 I=1,7 IF(MP.EQ.NPOINT(I)) GO TO 2 1 CONTINUE C C.....INVALID M USED. C WRITE(6,10) 10 FORMAT(' ','******INVALID M HAS BEEN USED.******') GAUSS=0.D0 RETURN C C.....SET UP INITIAL PARAMETERS...... C 2 JFIRST=KEY(I) JLAST=KEY(I+1)-1 C=(B-A)/2.D0 D=(B+A)/2.D0 C C.....ACCUMULATE THE SUM IN THE M-POINT FORMULA..... SUM=0.D0 INDEX=0 DO 5 J=JFIRST,JLAST IF(ZI(J).NE.0.D0) GO TO 8 SUM=SUM+WEIGHT(J)*F(D) GO TO 5 8 CONTINUE SUM=SUM+WEIGHT(J)*F(ZI(J)*C+D) SUM=SUM+WEIGHT(J)*F(-ZI(J)*C+D) 5 CONTINUE C C.....MAKE INTERVAL CORRECTION AND RETURN...... GAUSS=C*SUM RETURN END DOUBLE PRECISION FUNCTION F(PSISTR) C C THE FUNCTION F COMPUTES THE THREE FUNCTIONS F1,F2 AND F3,FOR C A GIVEN SPHERICAL DISTANCE "PSISTR".THEN, C F=F1*F2*F3. C M. NAJAFI, 10/1980 IMPLICIT REAL*8(A-H,O-Z) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CZ/ ZPP(20,16),INDEX,IZP,JZP,INFINT COMMON/PP/ PSI0,PSI,N,M INDEX=INDEX+1 INFT=INFINT F=0.D0 IF(PSISTR.EQ.0.D0.OR.PSISTR.EQ.PI) RETURN PI=DARCOS(-1.D0) DR=PI/180.D0 SCALE=180.D0/PSI0 PHI=PSISTR*SCALE/DR F1=DSIN(PHI*DR) F2=P(N,M,PHI) F3=ZPP(INDEX,M+1) CALL TRAPS(0,0,1000,0,0) F=F1*F2*F3 RETURN END DOUBLE PRECISION FUNCTION P(NN,MM,TETA) C C THE FUNCTION P TAKES OUT A VALUE OF ASSOCIATED LEGENDRE C FUNCTION CORRESPONDING TO THE GIVEN ANGLE "TETA",DEGREE "NN" C AND ORDER "MM" FROM THE TABLE (ARRAY) OF THE LEGENDRE FUNCTIONS C PNM(SEE THE MAIN PROGRAM) BY A LINEAR INTERPOLATION. C M. NAJAFI, 10/1980 IMPLICIT REAL*8(A-H,O-Z) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM ANG=TETA K=(NN+1)*(NN+2)/2.D0-NN+MM IF(TETA.GT.90.D0) GO TO 99 IR=ANG/PSINT+1 IF(IR.EQ.IPNM) GO TO 88 DANG=ANG-(IR-1)*PSINT P=PNM(IR,K)+(PNM(IR+1,K)-PNM(IR,K))*DANG/PSINT GO TO 77 88 P=PNM(IR,K) 77 RETURN 99 ANG=180.D0-TETA IR=ANG/PSINT+1 DANG=ANG-(IR-1)*PSINT P=PNM(IR,K)+(PNM(IR+1,K)-PNM(IR,K))*DANG/PSINT P=P*(-1)**(NN+MM) RETURN END DOUBLE PRECISION FUNCTION Q(NN,MM,TETA) C C THE FUNCTION Q TAKES OUT A VALUE CORRESPONDING TO A GIVEN C ANGLE "TETA",DEGREE "NN" AND ORDER "MM" FROM THE TABLE OF C THE F FUNCTIONS (THE ARRAY QNM,SEE THE MAIN PROGRAM). C NM C M. NAJAFI, 10/1980 IMPLICIT REAL*8 (A-H,O-Z) COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CQ/ QNM(121,45),QSINT1,QSINT2,PSI00,IQ1,IQ2,IQNM,JQNM K=(NN+1)*(NN+2)*1.D0/2-NN+MM QSINT=QSINT1 IR=TETA/QSINT+1 DTETA=TETA-(IR-1)*QSINT IF(TETA.LE.PSI00) GO TO 100 QSINT=QSINT2 IR=(TETA-PSI00)/QSINT+IQ1 DTETA=TETA-((IR-IQ1)*QSINT+PSI00) 100 Q=QNM(IR,K) IF(IR.EQ.IQNM) RETURN Q=QNM(IR,K)+(QNM(IR+1,K)-QNM(IR,K))*DTETA/QSINT RETURN END DOUBLE PRECISION FUNCTION H(ID,IO) C C THE FUNCTION H COMPUTES A MULTIPLIER,TO NORMALIZE AN C ASSOCIATED LEGENDRE FUNCTION OF DEGREE "ID" AND ORDER "IO". C THAT IS, _ C P =H*P C M. NAJAFI, 10/1980 IMPLICIT REAL*8(A-H,O-Z) IF(IO.NE.0) GO TO 5 H=DSQRT(2.D0*ID+1.D0) RETURN 5 NPM=ID+IO XX=2*(2.D0*ID+1.D0) ITERMS=2*IO DO 10 I=1,ITERMS XX=XX/NPM 10 NPM=NPM-1 H=DSQRT(XX) RETURN END SUBROUTINE REGOID(BIGF,SMLF,BIGL,SMLL,DF,DL,A,KA,CN,SINCF,COSCF, &CF,CL,NMAX,PSI0) C C SUBROUTINE REGOID COMPUTES,AND PRITS OUT THE GEOID. C INPUT DATA: A(KA)- THE COMPUTED COEFFICIENTS IN PART II. C BIGF,SMLF,BIGL,SMLL- SIDES OF THE AREA OF INTEREST C GIVEN IN TWO PARALLELS AND TWO MERIDIANS (DEGREES). C DF,DL- INCREMENTS IN LATITUDE AND LONGITUDE(DEGREES). C SINCF,COSCF- SIN AND COS OF THE CENTRAL LATITUDE. C CF,CL- LATITUDE AND LONGITUDE OF THE CENTRE (RADIANS). C CN- GEOIDAL HEIGHT OF THE CENTRE. C NMAX- MAXIMUM LIMIT UP TO WHICH THE GEOID IS COMPUTED. C PSI0- RADIUS OF THE CAP (DEGREES). C OUTPUT DATA: 1- A LIST OF THE COEFFICIENTS IS PRINTED. C 2- THE GEOIDAL HEIGHTS ARE PRINTED. C M. NAJAFI, 9/1980 IMPLICIT REAL*8 (A-H,O-Z) REAL*8 PN(20,20),SHR(100),GH(20,20),A(KA),LO(20) INTEGER C(100) IPN=20 NP1=NMAX+1 NNP1=NP1 IF(NNP1.LT.5) NNP1=5 PI=DARCOS(-1.D0) TWOPI=2*PI DR=PI/180.D0 PIPSI=180.D0/PSI0 C C COMPUTING THE GEOID FOR THE GIVEN AREA C K=1 C(K)=100 IF(NMAX.EQ.0) GO TO 100 DO 10 N=1,NMAX K=K+1 C(K)=100+N*10 DO 11 M=1,N K=K+1 C(K)=100+N*10+M K=K+1 C(K)=200+N*10+M 11 CONTINUE 10 CONTINUE 100 CONTINUE IF(K.EQ.KA) GO TO 60 WRITE(6,65) RETURN 60 CONTINUE PHI=BIGF IP=0 40 CONTINUE IP=IP+1 XLM=SMLL JP=0 41 CONTINUE JP=JP+1 DXLON=XLM*DR-CL+TWOPI PSI=DARCOS(SINCF*DSIN(PHI*DR)+COSCF*DCOS(PHI*DR)* + DCOS(DXLON)) IF(PSI.GT.1.D-5) GO TO 45 GH(IP,JP)=CN GO TO 46 45 CONTINUE COSAZ=(DSIN(PHI*DR)-DCOS(PSI)*SINCF)/DSIN(PSI)/COSCF IF(COSAZ.LT.0.D0.AND.DABS(COSAZ+1.D0).LT.1.D-10) COSAZ=-1.D0 IF(COSAZ.GT.0.D0.AND.DABS(COSAZ-1.D0).LT.1.D-10) COSAZ=1.D0 AZ=DARCOS(COSAZ) IF(DXLON.GT.PI.AND.DXLON.LT.TWOPI) AZ=TWOPI-AZ FI=(PSI*PIPSI)/DR FII=90.0-FI CALL LEGDRE(1,FII,PN,IPN,NNP1) K=1 SHR(K)=PN(1,1) IF(NMAX.EQ.0) GO TO 200 DO 42 N=1,NMAX K=K+1 SHR(K)=PN(N+1,1) DO 43 M=1,N K=K+1 SHR(K)=DCOS(M*AZ)*PN(N+1,M+1) K=K+1 43 SHR(K)=DSIN(M*AZ)*PN(N+1,M+1) 42 CONTINUE 200 CONTINUE C DN=0.D0 DO 44 I=1,K 44 DN=DN+A(I)*SHR(I) GH(IP,JP)=DN 46 CONTINUE XLM=XLM+DL IF(XLM.LE.BIGL) GO TO 41 PHI=PHI-DF IF(PHI.GE.SMLF) GO TO 40 C C PRINT OUT THE RESULTS C CFD=CF/DR CLD=CL/DR WRITE(6,35) CFD WRITE(6,36) CLD WRITE(6,33) CN WRITE(6,31) DO 17 I=1,K IC=C(I)/10.0 IF(C(I).GE.200) GO TO 61 I1=IC-10 I2=C(I)-IC*10 WRITE(6,32) I,I1,I2,A(I) GO TO 17 61 CONTINUE I1=IC-20 I2=C(I)-IC*10 WRITE(6,66) I,I1,I2,A(I) 17 CONTINUE WRITE(6,53) DO 90 J=1,JP 90 LO(J)=SMLL+(J-1)*DL WRITE(6,50) (LO(J),J=1,JP) DO 51 I=1,IP PHI=BIGF-(I-1)*DF WRITE(6,52) PHI,(GH(I,J),J=1,JP) 51 CONTINUE 31 FORMAT('0',10X,'COEFFICIENTS',15X,'METERS') 32 FORMAT('0',7X,I3,5X,'A',I1,I1,9X,F18.4) 33 FORMAT('0',15X,'N0',10X,F16.2) 35 FORMAT('1',10X,'LATITUDE OF THE CENTRE ',F7.2,1X,'DEG.S') 36 FORMAT('0',10X,'LONGITUDE OF THE CENTRE',F7.2,1X,'DEG.S') 50 FORMAT('0 LONGITUDE',9F8.1/' ---------'/' LATITUDE') 52 FORMAT('0',F9.1,4X,9F8.2) 53 FORMAT('1',' THE REGIONAL GEOID') 65 FORMAT('0','****(K.NE.KA)****') 66 FORMAT('0',7X,I3,5X,'B',I1,I1,9X,F18.4) RETURN END DOUBLE PRECISION FUNCTION SF(SDIS) C C FUNCTION SF COMPUTES THE STOKES FUNCTION FOR A GIVEN SPHERICAL C DISTANCE SDIS IN RADIANS. C M. NAJAFI, 9/1980 IMPLICIT REAL*8 (A-H,O-Z) SSDIS=DSIN(SDIS) CSDIS=DCOS(SDIS) SSDI2=DSIN(SDIS/2.) CSDI2=DCOS(SDIS/2.) SF=1/SSDI2+1-6*SSDI2-CSDIS*(5+3*DLOG(SSDI2+SSDI2**2)) RETURN END SUBROUTINE STKSM(SINLA,COSLA,XLATA,XLONA,XF,XL,MLAT,MLON,SUMA) C C THIS SUBROUTINE COMPUTES: 1- A MEAN VALUE OF THE STOKES FUNCTION C FOR A BLOCK OF SIZE (MLAT X MLON) BY DEVIDING IT INTO SMALL BLOCKS C OF SIZE (IDGXY X IDGXY). 2- SUM OF THE GRAVITY ANOMALIES INSIDE C THE INNERMOST ZONE "DGIN" AND THEIR NUMBER "NINP". C INPUT DATA : C SINLA,COSLA_ ARE SIN , COS OF LATITUDE OF THE COMPUTATION POINT. C C XLATA,XLONA_ ARE LATITUDE,LONGITUDE OF THE COMPUTATION POINT C IN RADIANS. C XF,XL- THE LATITUDE AND LONGITUDE OF THE BLOCK (MLAT X MLON) C IN DEGREES. C C OUTPUT DATA : C SUMA_ THE MEAN VALUE OF THE STOKES FUNCTION FOR THE BLOCK. C DGIN,NINP- SUM OF THE GRAVITY ANOMALIES AND THIER NUMBER. C NOTE: THE IDEA,TO WRITE THE ROUTINE,HAS BEEN ISSUED FROM C LACHAPELLE 1976. M. NAJAFI, 9/1980 IMPLICIT REAL*8 (A-H,O-Z) COMMON/SKS/ STKA(4500),STKB(2000),STKC(900),PSIA,PSIB,PSIC, &HLAT20,HLON20,DG,DGIN,DGINA, &HLTA,HLNA,HLTB,HLNB,HLTC,HLNC,INTSTA,INTSTB,INTSTC, &IDGXYA,IDGXYB,IDGXYC,IDGXYD,NINP,NINPA PI=DARCOS(-1.D0) DR=PI/180.D0 SCRD=206264.80625 DCLT=DABS(XLATA/DR-XF) DCLN=DABS(XLONA/DR-XL) IF(DCLT.GT.HLTA.OR.DCLN.GT.HLNA) GO TO 90 IDGXY=IDGXYA GO TO 93 90 IF(DCLT.GT.HLTB.OR.DCLN.GT.HLNB) GO TO 91 IDGXY=IDGXYB GO TO 93 91 IF(DCLT.GT.HLTC.OR.DCLN.GT.HLNC) GO TO 92 IDGXY=IDGXYC GO TO 93 92 IDGXY=IDGXYD 93 CONTINUE DCLT=DCLT*DR DCLN=DCLN*DR DDGXY=IDGXY/SCRD ISLAT=MLAT*60 ISLON=MLON*60 NLATS=ISLAT*1.D0/IDGXY NLONS=ISLON*1.D0/IDGXY XFR=XF*3600/SCRD XLR=XL*3600/SCRD XLAT=XFR-(ISLAT/2.0)/SCRD+DDGXY/2.0 SLON=XLR-(ISLON/2.0)/SCRD+DDGXY/2.0 SUMA=0.D0 NT=0 IF(DCLT.LT.HLAT20.AND.DCLN.LT.HLON20) GO TO 10 DO 25 I=1,NLATS XLON=SLON DO 20 J=1,NLONS PSI=(DARCOS(SINLA*DSIN(XLAT)+COSLA*DCOS(XLAT)*DCOS(XLONA- + XLON)))*SCRD IF(PSI.GT.PSIA) GO TO 50 K=PSI/INTSTA+0.5D0 SUMA=SUMA+STKA(K) GO TO 52 50 IF(PSI.GT.PSIB) GO TO 51 K=PSI/INTSTB+0.5D0 SUMA=SUMA+STKB(K) GO TO 52 51 K=PSI/INTSTC+0.5D0 SUMA=SUMA+STKC(K) 52 CONTINUE NT=NT+1 20 XLON=XLON+DDGXY 25 XLAT=XLAT+DDGXY SUMA=SUMA/NT RETURN 10 CONTINUE C C COMPUTING THE SUM OF THE GRAVITY ANOMALIES WITHIN THE C INNERMOST ZONE. C DGIN=DGIN+DG NINP=NINP+1 RETURN END SUBROUTINE CSKNM(SINLA,COSLA,XLATA,XLONA,XF,XL,MLAT,MLON,NMAX, &AVKLS,NK) C C S C (CSKNM) COMPUTES MEAN VALUES OF THE KERNELS K ,K (N=0...NMAX) C NM NM C IN A BLOCK OF SIZE (MLAT X MLON), AND SUM OF THE MEAN ANOMALIES C INTO THE INNERMOST ZONE AND THEIR NUMBER. C INPUT DATA : C SINLA,COSLA_ ARE SIN,COS OF THE CENTRE. C XLATA,XLONA_ LATITUDE,LONGITUDE OF THE CENTRE (RADIANS). C XF,XL_ LATITUDE,LONGITUDE OF THE CENTRE OF THE BLOCK GIVEN C IN DEGREES. C MLAT,MLON_ SIDES OF THE BLOCK GIVEN IN MINUTES. C OUTPUT DATA: C AVKLS(NK)_ AN ARRAY CONTAINING THE MEAN VALUES. C M. NAJAFI 10/1980 IMPLICIT REAL*8 (A-H,O-Z) REAL*8 AVKLS(100) COMMON/CF/ COF COMMON/CP/ PNM(541,66),PSINT,NCAP,IPNM,JPNM COMMON/CQ/ QNM(121,45),QSINT1,QSINT2,PSI00,IQ1,IQ2,IQNM,JQNM COMMON/SKS/ STKA(4500),STKB(2000),STKC(900),PSIA,PSIB,PSIC, &HLAT20,HLON20,DG,DGIN,DGINA, &HLTA,HLNA,HLTB,HLNB,HLTC,HLNC,INTSTA,INTSTB,INTSTC, &IDGXYA,IDGXYB,IDGXYC,IDGXYD,NINP,NINPA PI=DARCOS(-1.D0) DR=PI/180.D0 SCRD=206264.80625 DCLT=DABS(XLATA/DR-XF) DCLN=DABS(XLONA/DR-XL) IF(DCLT.GT.HLTA.OR.DCLN.GT.HLNA) GO TO 90 IDGXY=IDGXYA GO TO 93 90 IF(DCLT.GT.HLTB.OR.DCLN.GT.HLNB) GO TO 91 IDGXY=IDGXYB GO TO 93 91 IF(DCLT.GT.HLTC.OR.DCLN.GT.HLNC) GO TO 92 IDGXY=IDGXYC GO TO 93 92 IDGXY=IDGXYD 93 CONTINUE DCLT=DCLT*DR DCLN=DCLN*DR PI2=PI/2 TWOPI=2*PI DDGXY=IDGXY/SCRD ISLAT=MLAT*60 ISLON=MLON*60 NLATS=ISLAT*1.D0/IDGXY NLONS=ISLON*1.D0/IDGXY XFR=XF*3600/SCRD XLR=XL*3600/SCRD XLAT=XFR-(ISLAT/2.0)/SCRD+DDGXY/2.0 SLON=XLR-(ISLON/2.)/SCRD+DDGXY/2.0 DO 20 I=1,100 20 AVKLS(I)=0.D0 NT=0 IF(DCLT.LT.HLAT20.AND.DCLN.LT.HLON20) GO TO 10 DO 24 I=1,NLATS XLON=SLON DO 25 J=1,NLONS DXLON=XLON+2*PI-XLONA PSI=DARCOS(SINLA*DSIN(XLAT)+COSLA*DCOS(XLAT)*DCOS(DXLON)) COSAZ=(DSIN(XLAT)-DCOS(PSI)*SINLA)/DSIN(PSI)/COSLA IF(COSAZ.LT.0.D0.AND.DABS(COSAZ+1.D0).LT.1.D-10) COSAZ=-1.D0 IF(COSAZ.GT.0.D0.AND.DABS(COSAZ-1.D0).LT.1.D-10) COSAZ=1.D0 AZ=DARCOS(COSAZ) IF(DXLON.GT.PI.AND.DXLON.LT.TWOPI) AZ=TWOPI-AZ PSI=PSI/DR K=1 AVKLS(K)=AVKLS(K)+COF*Q(0,0,PSI*COF) IF(NMAX.EQ.0) GO TO 100 DO 26 N=1,NMAX K=K+1 AVKLS(K)=AVKLS(K)+COF*Q(N,0,PSI*COF) DO 27 M=1,N K=K+1 TKR=COF*Q(N,M,PSI*COF) AVKLS(K)=AVKLS(K)+DCOS(M*AZ)*TKR K=K+1 27 AVKLS(K)=AVKLS(K)+DSIN(M*AZ)*TKR 26 CONTINUE 100 CONTINUE C NT=NT+1 25 XLON=XLON+DDGXY 24 XLAT=XLAT+DDGXY NK=K DO 28 I=1,NK 28 AVKLS(I)=AVKLS(I)/NT RETURN 10 CONTINUE C C COMPUTING THE SUM OF THE MEAN GRAVITY ANOMALIES WITHIN THE C INNERMOST ZONE. C DGINA=DGINA+DG NINPA=NINPA+1 RETURN END SUBROUTINE LEGDRE(INORM,PHI,PN,NROW,NP1) C C C SUBROUTINE LEGDRE COMPUTES THE ASSOCIATED LEGENDRE FUNCTIONS C UP TO AND INCLUDING DEGREE AND ORDER N TO BE SPECIFIED. THE C DIMENSION OF PN IS PN(N+1,N+1). THE ASSOCIATED LEGENDRE POLY- C NOMIAL OF DEGREE A AND ORDER B IS STORED IN PN(A+1,B+1) IF A.NE.B C (ZONAL AND TESSERAL) AND IN PN(A+1,A+1) IF A = B (SECTORIAL) C C INORM = FLAG FOR NORMALIZATION: 1 = YES; 0 = NO C PHI = LATITUDE IN DEGREES C NROW = ROW DIMENSION OF PN IN THE CALLING PROGRAM C NP1= N + 1 = THE DIMENSION OF PN IN THE SUBROUTINE. SEE N BELOW C N = DEGREE OF LEGENDRE POLYNOMIAL C C C IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION PN(NROW,NP1) DATA DR/0.017453293D0/ N=NP1-1 N1=NP1 PHI=PHI*DR X =DSIN(PHI) DO 9 I=1,NROW DO 9 J=1,NP1 9 PN(I,J) = 0.0D0 PN(1,1) = 1.0D0 PN(2,1) = X PN(2,2) = 1.0D0 DO 15 I = 2,N I1 = I + 1 DO 10 J = 1,I1 IF(J .NE. 1) GO TO 11 PN(I1,J) = (1./I)*(X*(2.*I-1)*PN(I1-1,J) - (I-1)*PN(I1-2,J)) GO TO 10 11 IF(J.NE.2 .AND. I.NE.2) GO TO 12 PN(I1,J) = PN(I1-2,J) + (2.*I-1)*PN(I1-1,J-1) GO TO 10 12 IF(J.NE.(I1-1).AND.J.NE.I1) GO TO 13 PN(I1,J) = (2.*I-1)*PN(I1-1,J-1) GO TO 10 13 PN(I1,J) = PN(I1-2,J) + (2.*I-1)*PN(I1-1,J-1) 10 CONTINUE 15 CONTINUE DO 16 I = 1,N I1 = I + 1 DO 16 J = 1,I J1 = J + 1 COSPHI = DCOS(PHI) IF(COSPHI.EQ.0.D0) GO TO 20 C CALL TRAPS(0,0,1000,0,0) GO TO 21 20 PN(I1,J1) = 0.D0 GO TO 16 21 PN(I1,J1) = (COSPHI**J)*PN(I1,J1) 16 CONTINUE IF(INORM.EQ.0) GO TO 14 N1 = N + 1 DO 17 I = 1,N1 17 PN(I,1) = DSQRT(2.D0*(I-1) + 1)*PN(I,1) DO 18 I = 1,N I1 = I + 1 DO 18 J = 1,I J1 = J + 1 NMM = I - J NPM = I + J FTERM = FACT(NMM) / FACT(NPM) NPOWER = NMM - NPM IF(MOD(NPOWER,2).EQ.0) GO TO 502 NPOWER = NPOWER + 1 FTERM = FTERM / 10D0 502 FTERM = DSQRT(FTERM) NPOWER = NPOWER / 2 FTERM = FTERM*10.**NPOWER 18 PN(I1,J1) = DSQRT(2.D0*(2.D0*I+1.D0))*FTERM*PN(I1,J1) 14 RETURN END FUNCTION FACT(N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION FA(59),FB(59),FMAN(118),NEXP(118) EQUIVALENCE (FA(1),FMAN(1)),(FB(1),FMAN(60)) DATA FA/0.1,0.2,0.6,0.24,0.12,0.72,0.504,0.4032,0.36288,0.36288, &0.399168,0.4790016,0.62270208,0.871782912,0.1307674368, &0.20922789888,0.355687428096,0.6402373705728,0.1216451004088320, &0.2432902008176640,0.5109094217170944,0.1124000727777608, &0.2585201673888497,0.6204484017332395,0.1551121004333098, &0.4032914611266056,0.1088886945041835,0.3048883446117138, &0.8841761993739699,0.2652528598121910,0.8222838654177920, &0.2631308369336934,0.8683317618811883,0.2952327990396040, &0.1033314796638614,0.3719933267899011,0.1376375309122634, &0.5230226174666009,0.2039788208119743,0.8159152832478973, &0.3345252661316379,0.1405006117752879,0.6041526306337380, &0.2658271574788447,0.1196222208654801,0.5502622159812085, &0.2586232415111680,0.1241391559253606,0.6082818640342671, &0.3041409320171335,0.1551118753287381,0.8065817517094379, &0.4274883284060022,0.2308436973392412,0.1269640335365826, &0.7109985878048627,0.4052691950487717,0.2350561331282876, &0.1386831185456896/ DATA FB/0.8320987112741378,0.5075802138772241,0.3146997326038789, &0.1982608315404437,0.1268869321858840,0.8247650592082458, &0.5443449390774422,0.3647111091818862,0.2480035542436826, &0.1711224524281410,0.1197857166996987,0.8504785885678606, &0.6123445837688597,0.4470115461512676,0.3307885441519380, &0.2480914081139535,0.1885494701666046,0.1451830920282856, &0.1332428117820627,0.8946182130782955,0.7156945704626364, &0.5797126020747355,0.4753643337012831,0.3945523969720649, &0.3314240134565345,0.2817104114380543,0.2422709538367267, &0.2107757298379522,0.1854826422573979,0.1650795516090842, &0.1485715964481757,0.1352001527678399,0.1243841405464127, &0.1156772507081638,0.1087366156656740,0.1032997848823903, &0.9916779348709466,0.9619275968248181,0.9426890448883216, &0.9332621544394384,0.9332621544394384,0.9425947759838328, &0.9614466715035094,0.9902900716486147,0.1029901674514559, &0.1081396758240287,0.1146280563734704,0.1226520203196134, &0.1324641819451824,0.1443859583202488,0.1588245541522737, &0.1762952551090238,0.1974506857221067,0.2231192748659805, &0.2543559733472178,0.2925093693493005,0.3393108684451886, &0.3969937160808705,0.4684525849754272/ DATA NEXP/1,1,1,2,3,3,4,5,6,7,8,9,10,11,13,14,15,16,17,19,20,22, &23,24,26,27,29,30,31,33,34,36,37,39,41,42,44,45,47,48,50,52,53, &55,57,58,60,62,63,65,67,68,70,72,74,75,77,79,81,82,84,86,88,90,91, &93,95,97,99,101,102,104,106,108,110,112,114,116,117,119,121,123, &125,127,129,131,133,135,137,139,141,143,145,147,148,150,152,154, &156,158,160,162,164,167,169,171,173,175,177,179,181,183,185,187, &189,191,193,195/ C IF(N.LT.0) GO TO 51 IF(N.EQ.0) N = 1 IF(N.GT.118) GO TO 52 FACT = FMAN(N) N = NEXP(N) RETURN 51 WRITE(6,301) N 301 FORMAT(/////,'FACTORIAL FUNCTION HAS BEEN GIVEN A NEGATIVE ARGUMEN #T: N=',I6) STOP 52 WRITE(6,302) N 302 FORMAT(/////,'ARGUMENT OF THE FACTORIAL FUNCTION IS GREATER THAN #118: N=',I10) STOP END