//GPOT JOB NOTIFY=6461 /*SERVICE -4 /*JOBPARM L=19,S=29,R=2048,PRINT=ALL //STEP1 EXEC FORTVCLG,REGION=2048K //FORT.SYSIN DD * C--------------------------------------------------------------------- C ****************** PROGRAM GPOT *********************** C TEST PROGRAM USING TSCHERNING'S HARMONIC MODEL SUBROUTINE C PROGRAM TO COMPUTE GEOID UNDULATIONS , DEFLECTIONS OF THE C VERTICAL AND GRAVITY DISTURBANCES FROM POTENTIAL COEFFICIENTS C********************************************************************* C C REFERENCE : MANUSCRIPTA GEODAETICA VOL. 8 (1983) 249-272 C C********************************************************************* C C GRID GEOID COMPUTATIONS C C INPUT: LAT1,LAT2 - MAXIMUM AND MINIMUM LATITUDE OF THE AREA C OF INTEREST (IN INTEGER DEGREES) C LON1,LON2 - MAXIMUM AND MINIMUM LONGITUDE OF THE AREA C OF INTEREST (IN INTEGER DEGREES) C INT - GRID INTERVAL (IN INTEGER MINUTES OF ARC) C NMAX - DEGREE AND ORDER OF SPHERICAL HARMONIC C EXPANSION DESIRED C C OUTPUT: UND - GEOIDAL HEIGHT C XI, ETA - N-S AND E-W DEFLECTIONS OF THE VERTICAL C DG - GRAVITY DISTURBANCE C C--------------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) C--------------------------------------------------------------------- C INPUT LOGICAL UNIT FOR EXTERNAL INPUT IS : 5 C INPUT LOGICAL UNIT FOR POTENTIAL COEFFICIENTS IS : 12 C--------------------------------------------------------------------- ICR = 5 C--------------------------------------------------------------------- C READ THE MAXIMUM AND MINIMUM LATITUDE AND LONGITUDE C OF THE AREA OF GEOID COMPUTATION IN INTEGER DEGREES C READ THE GRID INTERVAL IN INTEGER MINUTES OF ARC C--------------------------------------------------------------------- READ(ICR,100) LAT1,LAT2,LON1,LON2,INT 100 FORMAT(5I4) C LATR =((LAT1-LAT2)*(60/INT))+1 LONR =((LON1-LON2)*(60/INT))+1 WRITE(6,100) LAT1,LAT2,LON1,LON2,LATR,LONR DINCR=INT/60.0D0 C LATST = 1 DO 20 LAT=LATST,LATR PHI = DFLOAT(LAT2)+((LAT-1)*DINCR) C DO 20 LON=1,LONR DLAM = DFLOAT(LON2)+((LON-1)*DINCR) C CALL POT1(PHI,DLAM,0.D0,UN1,XI1,ETA1,DG1) C--------------------------------------------------------------------- C IF YOU WISH YOU CAN CONVERT THE GRAVITY DISTURBANCE INTO C GRAVITY ANOMALY USING THE FOLLOWING FORMULA: C DG1=DG1-0.3086D0*UN1 C--------------------------------------------------------------------- WRITE (6,111) PHI,DLAM,UN1,DG1 111 FORMAT(4F10.4) C 20 CONTINUE STOP END C********************************************************************* BLOCK DATA IMPLICIT REAL*8(A-H,O-Z) LOGICAL INIT LOGICAL FIRST INTEGER OLDORD REAL*4 C,C0 COMMON/GPOTCM/OLDT,OLDR,IZ,FIRST,OLDORD,I1,I2,I3,I4, * I5,I6,I7,I8,I9,NMAXSV COMMON/POT1CM/SU(1810),DJN(20),GM,FLAT,AE,OMEGA,INIT,IORDER,NMAX, * NEGN COMMON/PIDTR/PI,DTR COMMON/TRANCM/TOL,MAXIT COMMON/CM/C20IN,G1(3),G2(3,3),CM3,CM2,CM1,C0,C(32760) DATA IZ/0/ DATA FIRST/.FALSE./ DATA OLDT/0.D0/,OLDR/0.D0/ DATA OLDORD,I1,I2,I3,I4,I5,I6,I7,I8,I9/10*0/ DATA NMAXSV/0/ DATA DJN/20*0.D0/ DATA INIT/.TRUE./ DATA IORDER/1/ DATA NMAX/180/ DATA PI/3.141592653589793D0/ DATA DTR/0.1745329251994330D-1/ DATA TOL/1.D-14/,MAXIT/10/ END C********************************************************************* SUBROUTINE POT1(PHI,DLON,HT,UN,XI,ETA,DIST) C IMPLICIT REAL*8 (A-H,O-Z) REAL*4 C,C0 LOGICAL INIT DIMENSION P(6) COMMON/PIDTR/PI,DTR COMMON/POT1CM/SU(1810),DJN(20),GM,FLAT,AE,OMEGA,INIT,IORDER,NMAX, & NEGN COMMON/CM/C20IN,G1(3),G2(3,3),CM3,CM2,CM1,C0,C(32760) C************************************************************** C * C INPUT * C * C PHI LATITUDE (GEODETIC) IN DEGREES * C DLON LONGITUDE (WEST) IN DEGREES * C HT HEIGHT IN METERS * C * C OUTPUT * C * C UN HEIGHT ANOMALY IN METERS * C XI N-S DEFLECTION IN SEC OF ARC * C ETA E-W DEFLECTION IN SECONDS OF ARC (WEST POSITIVE) * C DIST GRAVITY DISTURBANCE IN MGALS * C * C************************************************************** IF(.NOT.INIT) GO TO 500 INIT=.FALSE. NEGN=-NMAX C C SET CONSTANTS (GRS80) C DJN(2)=0.00108263D0 OMEGA=7.292115D-5 AE=6378137.D0 GM=3.986005D+14 C C OBTAIN FLATTENING OF REFERENCE ELLIPSOID TO BE USED LATER WHEN C TRANSFORMING FROM GEODETIC LAT,LON,HT TO ECG X,Y,Z C CALL REFVAL(AE,GM,DJN(2),OMEGA,FLATI) FLAT=1.D0/FLATI C C GENERATE NORMAL POTENTIAL VALUES C CAPESQ=AE**2*(2.D0-FLAT)*FLAT ESQ=CAPESQ/AE**2 C C COMPUTE NORMAL EVEN ZONALS FROM 4 TO 20 C DO 200 N2=4,20,2 N=N2/2 DJN(N2)=(-1.D0)**(N+1)*3.D0/(N2+1.D0)*ESQ**N/(N2+3.D0)* & (1.D0-(1.D0-5.D0*DJN(2)/ESQ)*N) 200 CONTINUE CM3=GM CM2=AE CM1=0.D0 C0=0.E0 DO 10 I=1,32760 C(I)=0.E0 10 CONTINUE CALL LOADCS(NMAX) C0=0.D0 C C COMPUTE DELC20 IN DOUBLE PRECISION BEFORE STORING IN SINGLE C PRECISION RATHER THAN STORING IN SINGLE PRECISION AND THEN C DIFFERENCING C DELC20=C20IN+DJN(2)/DSQRT(5.D0) C CALL STORC(0,0,0.D0,0.D0) CALL STORC(1,0,0.D0,0.D0) CALL STORC(1,1,0.D0,0.D0) CALL STORC(2,0,DELC20,0.D0) C C NOW SUBTRACT OFF GRS80 EVEN ZONAL HARMONICS ( J2 ALREADY DONE ) C TO GET THE ANOMALOUS POTENTIAL FOR HEIGHT ANOMALY AND DERIVATIVE C COMPUTATIONS C DO 400 N=4,20,2 C NORMALIZE ZONALS OF NORMAL POTENTIAL BEFORE SUBTRACTING DNJ=DJN(N)/DSQRT(N+N+1.D0) CALL MODC(N,0,DNJ,0.D0) 400 CONTINUE C CALL SETCM(NMAX) C C CONVERT FROM GEODETIC TO EARTH-CENTERED-FIXED X,Y,Z COORDINATES C 500 CALL TRANF(1,PHI,DLON,HT,X,Y,Z,AE,FLAT) C C SET UP INPUT ARRAY TO GPOTDR C P(1)=DSQRT(X**2+Y**2) P(2)=DSQRT(X**2+Y**2+Z**2) P(3)=Z/P(2) P(4)=P(1)/P(2) P(5)=Y/P(1) P(6)=X/P(1) C TP=GPOTDR(P,NEGN,IORDER,SU) C CALL NORMAL(GM,DJN,AE,OMEGA,X,Y,Z,UP,GAMMA,GRAD) C THE HEIGHT ANOMALY IN METRES IS C UN=TP/GAMMA C C DEFLECTIONS OF THE VERTICAL IN SEC OF ARC C N-S XI =-G1(1)*206264.8D0 / GAMMA C E-W ETA =-G1(2)*206264.8D0 / GAMMA C GRAVITY DISTURBANCE IN MILLIGALS DIST =-G1(3)*1.D+05 RETURN END C********************************************************************* FUNCTION GPOTDR(PO,NMAX,ORDER,SU) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C GI REG.NO. 80039 AUTHOR -C.C.TSCHERNING, NOV 1980 IN ALGOL C -C.C.GOAD, MAR 1981 TRANSLATED TO FORTRAN C C REFERENCES: C (1) TSCHERNING, C.C.:ON THE CHAIN-RULE METHOD FOR COMPUTING C POTENTIAL DERIVATIVES. MANUSCRIPTA GEODAETICA, VOL.1, C PP. 125-141, 1976 C C (2) GERSTL,M.:VERGLEICH VON ALGORITMEN ZUR SUMMATION VON C KUGELFLAECHENFUNKTIONEN. VEROEFFENTL. DER BAYER. KOMM. C F.D. INT. ERDMESSUNG DER BAYER. AKADEMIE DER WISSEN., C HEFT NR. 38, PP. 81-88, 1978. C THE PROCEDURE COMPUTES THE VALUE AND UP TO THE SECOND ORDER C DERIVATIVES OF THE POTENTIAL OF THE EARTH (W) OR OF ITS C CORRESPONDING AMONALOUS POTENTIAL(T). (THE COMPUTATION OF C THE SECOND ORDER DERIVATIVES HAS NOT YET BEEN IMPLEMENTED). C C THE POTENTIAL IS REPRESENTED BY A SERIES IN SOLID SPHERICAL C HARMONICS, WITH UN-NORMALIZED OR QUASI-NORMALIZED COEFFICIENTS. C THE CHAIN-RULE IS USED COMBINED WITH THE CLENSHAW ALGORITHM. C THE ARRAY C MUST HOLD THE COEFFICIENTS, C(0,0)=1.D0 FOR W AND C 0.0 FOR T, C(1)=C(1,0),C(2)=C(1,1),C(3)=S(1,1) ETC. UP TO C C((N+1)**2-1) = S(N,N). C C C PARAMETERS: C C (A) INPUT VALUES: C C NMAX C THE ABSOLUTE VALUES OF NMAX IS EQUAL TO THE MAXIMAL DEGREE AND C ORDER OF THE SERIES. NEGATIVE NMAX INDICATES THAT THE COEFFICIENTS C ARE QUASI-NORMALIZED. C C IORDER C THE MAXIMAL ORDER OF THE DERIVATIVES (< 2 P.T.). C C PO C ARRAY HOLDING POSITION INFORMATION. PO(6) C PO(1)=P, THE DISTANCE FROM THE Z (ROTATION) AXIS, C PO(2)=R, THE DISTANCE FROM THE ORIGIN, C PO(3),PO(4) COS AND SIN OF THE GEOCENTRIC POLAR ANGLE(COLATITUDE), C PO(5),PO(6) SIN AND COS OF THE LONGITUDE. C C C C C MUST BE DECLARED WITH BOUNDS (-3: (N+1)**2-1) WHEN THE C COEFFICIENTS ARE UN-NORMALIZED AND WITH BOUNDS (-3: (N+3)**2-2) C WHEN THE COEFFICIENTS ARE QUASI-NORMALIZED. C(1) TO C((N+1)**2-1) C CONTAIN THE COEFFICIENTS AND WE MUST HAVE C C(-3)=GM C C(-2)=A THE SEMI-MAJOR AXIS OF THE REF ELLIPSOID C C(-1)=THE ANGULAR VELOCITY (=0, WHEN DEALING WITH T). C C C SQUARE ROOT ARRAY C C C((N+1)**2+K) = SQRT(K), 0.LE.K.LE.2(ABS(N)+1)-1 WHEN N <0 C C C MOD--APRIL,1981 C WITH THE USE OF THE TABLE OF SQUARE ROOTS STORED IN ARRAY ROOT, C THE DIMENSION OF THE C ARRAY ONLY NEEDS TO BE (N+1)**2-1 FOR BOTH C UN-NORMALIZED AND NORMALIZED COEFFICIENTS. ARRAY ROOT MUST CONTAIN C SQUARE ROOT OF 0 TO N+2 C C C (B) RETURN VALUES C C G C THE RESULT IS STORED IN G AS FOLLOWS: C C G1(1)=DW/DX, G1(2)=DW/DY, G1(3)=DW/DZ C G(1,1)=DDW/DDX, G(1,2)=G(2,1)=DDW/DXDY, C G(1,3)=G(3,1)=DDW/DXDZ, G(2,2)=DDW/DDY, C G(2,3)=G(3,2)=DDW/DYDZ AND G(3,3)=DDW/DDZ C WHERE W MAY BE INTERCHANGED WITH T AND C VARIABLES X, Y, Z ARE THE CARTESIAN COORDINATES C IN A LOCAL (FIXED) FRAME WITH ORIGIN IN THE POINT C OF EVALUATION, X POSITIVE NORTH, Y POSITIVE EAST, C AND Z POSITIVE IN THE DIRECTION OF THE RADIUS C VECTOR, (CF. REF.(1),EQ (4) AND (5)). C THE VALUES OF W OR T WILL BE RETURNED IN POTCC. C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) INTEGER CAPN,ORDER,CAPN21,OLDORD LOGICAL QUASI,DERIV1,DERIV2,POLE LOGICAL FIRST,NEW,OLD,NPOLE REAL*4 C,C0 REAL*8 M21,M21T,M21U,M21U0 DIMENSION SML(181),CML(181),SMLP1(182),CMLP1(182),PO(6) DIMENSION SU(1810) COMMON/SQROOT/DZERO,ROOT(362) COMMON/GPOTCM/OLDT,OLDR,IZ,FIRST,OLDORD,I1,I2,I3,I4, & I5,I6,I7,I8,I9,NMAXSV COMMON/CM/C20IN,G1(3),G2(3,3),CM3,CM2,CM1,C0,C(32760) EQUIVALENCE(SML(1),SMLP1(2)),(CML(1),CMLP1(2)) IF(NMAXSV.NE.NMAX)FIRST=.FALSE. NMAXSV=NMAX IF(FIRST) GO TO 100 FIRST=.TRUE. OLDT=2.D0 J=IABS(NMAX) I=J+1 I1=I+1 I2=I1+I I3=I2+I I4=I3+I I5=I4+I I6=I5+I I7=I6+I I8=I7+I I9=I8+I 100 CAPN=NMAX P=PO(1) R=PO(2) T=PO(3) U=PO(4) SL=PO(5) CL=PO(6) T2=T+T POLE=DABS(U).LE.1.D-9 NEW=DABS(OLDR-R).GT.1.D-3.OR.DABS(OLDT-T).GT.1.D-9.OR. & OLDORD.NE.ORDER.OR.POLE OLD=.NOT.NEW NPOLE=.NOT.POLE IF(OLD) GO TO 200 OLDR=R OLDT=T OLDORD=ORDER C 200 QUASI=.FALSE. IF(CAPN.LT.0)QUASI=.TRUE. IF(QUASI)CAPN=-CAPN S=CM2/R S2=S**2 CMLP1(1)=1.D0 C CML(0)=1.D0 SMLP1(1)=0.D0 C SML(0)=0.D0 DERIV1=.FALSE. IF(ORDER.GT.0)DERIV1=.TRUE. DERIV2=.FALSE. IF( ORDER.GT.1)DERIV2=.TRUE. C C SML(M) AND CML(M) ARE THE SINE AND COSINE OF M*LONGITUDE C C WRITE(6,DBUG) C SML(1)=SL CML(1)=CL C M1=1 DO 300 M=2,CAPN SML(M)=SML(M1)*CL+CML(M1)*SL CML(M)=CML(M1)*CL-SML(M1)*SL 300 M1=M C CAPN21=CAPN+CAPN+1 VM=0.D0 VXM=0.D0 VYM=0.D0 VZM=0.D0 SQNM1=1.D0 SQNPM1=1.D0 IF(.NOT.DERIV2) GO TO 400 VXXM=0.D0 VYYM=0.D0 VZZM=0.D0 VXYM=0.D0 VXZM=0.D0 VYZM=0.D0 400 KM=(CAPN+1)**2 MAX2=CAPN21 C C WE NOW USE THE CLENSHAW ALGORITHM, CF. REF.(2), EQ(9), C MODIFIED IN AN OBVIOUS WAY FOLLOWING REF.(1). C ITWO=2 DO 1700 IM=IZ,CAPN M=CAPN-IM MPLUS1=M+1 IF(M.EQ.0)ITWO=1 KM=KM-ITWO K=KM N21=CAPN21 VS=0.D0 VC=0.D0 VS1=0.D0 VC1=0.D0 VXS1=0.D0 VXC1=0.D0 VZS=0.D0 VZC=0.D0 VZS1=0.D0 VZC1=0.D0 VXC=0.D0 VXS=0.D0 C IF(.NOT.DERIV2) GO TO 500 VXXC=0.D0 VXXS=0.D0 VXXC1=0.D0 VXXS1=0.D0 VZZC=0.D0 VZZS=0.D0 VZZC1=0.D0 VZZS1=0.D0 VXZC=0.D0 VXZS=0.D0 VXZC1=0.D0 VXZS1=0.D0 500 CM=CMLP1(MPLUS1) SM=SMLP1(MPLUS1) C NM1=CAPN-M+2 N1=CAPN+1 NPM1=CAPN+M+2 C IF(DERIV2)M2=M*M IF(OLD) GO TO 1300 C N=CAPN+1 DO 1000 IN=M,CAPN N=N-1 NM2=NM1 NM1=NM1-1 NPM1=NPM1-1 IF(.NOT.QUASI) GO TO 600 SQNM2=SQNM1 SQNM1=ROOT(NM1) SQNPM2=SQNPM1 SQNPM1=ROOT(NPM1) SQ1=SQNM1*SQNPM1 A1=S*N21/SQ1 B2=-S2*SQ1/(SQNM2*SQNPM2) GO TO 700 600 A1=S*N21/NM1 B2=-(S2*NPM1)/NM2 700 A1T=A1*T A1U=A1*U N21=N21-2 CK=C(K) CK1=C(K+1) K=K-N21 V2=VC1 VC1=VC VC=VC1*A1T+V2*B2+CK V2=VS1 VS1=VS VS=VS1*A1T+V2*B2+CK1 C IF(.NOT.DERIV1) GO TO 1000 CKZ=CK*N1 CK1Z=CK1*N1 V2=VXC1 VXC1=VXC VXC=VXC1*A1T+VC1*A1U+V2*B2 V2=VXS1 VXS1=VXS VXS=VXS1*A1T+VS1*A1U+V2*B2 V2=VZC1 VZC1=VZC VZC=VZC1*A1T+V2*B2-CKZ V2=VZS1 VZS1=VZS VZS=VZS1*A1T+V2*B2-CK1Z N1=N IF(.NOT.DERIV2) GO TO 1000 N2=N+2 V2=VZZC1 VZZC1=VZZC VZZC=VZZC1*A1T+V2*B2+N2*CKZ V2=VZZS1 VZZS1=VZZS VZZS=VZZS1*A1T+V2*B2+N2*CK1Z IF(NPOLE) GO TO 800 V2=VXXC1 VXXC1=VXXC VXXC=(VXXC1-VC1)*A1T+(A1U+A1U)*VXC1+V2*B2 V2=VXXS1 VXXS1=VXXS VXXS=(VXXS1-VS1)*A1T+(A1U+A1U)*VXS1+V2*B2 800 V2=VXZC1 VXZC1=VXZC VXZC=(VXZC1)*A1T+(A1U)*VZC1+V2*B2 V2=VXZS1 VXZS1=VXZS VXZS=(VXZS1)*A1T+(A1U)*VZS1+V2*B2 1000 CONTINUE SU(M+1)=VC SU(M+I1)=VS IF(.NOT.DERIV1) GO TO 1500 SU(M+I2)=VXC SU(M+I3)=VXS SU(M+I4)=VZC SU(M+I5)=VZS IF(.NOT.DERIV2) GO TO 1500 SU(M+I6)=VZZC SU(M+I7)=VZZS SU(M+I8)=VXZC SU(M+I9)=VXZS GO TO 1500 1300 VC=SU(M+1) VS=SU(M+I1) IF(.NOT.QUASI) GO TO 1400 SQNPM1=ROOT(MAX2) SQNPM2=ROOT(MAX2+1) 1400 NPM1=MAX2 MAX2=MAX2-2 IF(.NOT.DERIV1) GO TO 1500 VXC=SU(M+I2) VXS=SU(M+I3) VZC=SU(M+I4) VZS=SU(M+I5) IF(.NOT.DERIV2) GO TO 1500 VZZC=SU(M+I6) VZZS=SU(M+I7) VXZC=SU(M+I8) VXZS=SU(M+I9) C 1500 U0=U IF(M.EQ.0) U0=1.D0 AUX=NPM1 IF(QUASI) AUX=SQNPM1/SQNPM2 M21=S*AUX M21U=M21*U IF(.NOT.DERIV1) GO TO 1700 M21T=M21*T M21U0=M21*U0 IF(.NOT.DERIV2) GO TO 1600 VZZM=VZZC*CM+VZZS*SM+M21U*VZZM IF(M.GT.0) VXYM=M*(VXS*CM-VXC*SM)+M21U*VXYM-M21T*VYM VXZM=VXZC*CM+VXZS*SM-M21T*VZM+M21U*VXZM VYZM=(VZS*CM-VZC*SM)*M+M21U0*VYZM IF(POLE) VXXM=VXXC*CM+VXXS*SM+M21*(U*(VXXM-VM)-T2*VXM) IF(NPOLE) VYYM=-(VC*CM+VS*SM)*M2+M21U0*VYYM 1600 VXM=VXC*CM+VXS*SM-M21T*VM+M21U*VXM VYM=M*(VS*CM-VC*SM)+M21U0*VYM VZM=(VZC*CM+VZS*SM)+M21U*VZM 1700 VM=VC*CM+VS*SM+M21U*VM C C NOW ADD THE CONTRIBUTIONS FROM THE ROTATIONAL POTENTIAL C OM2=CM1**2 S=CM3/R GPOTDR=S*VM+OM2*P**2*.5D0 IF(.NOT.DERIV1) RETURN S=S/R G1(1)=S*VXM-T*P*OM2 G1(2)=S*VYM G1(3)=VZM*S+U**2*OM2*R IF(.NOT.DERIV2) RETURN S=S/R IF(NPOLE) GO TO 1900 VXXM=VXXM+VZM VYYM=-(VXXM+VZZM) GO TO 2000 1900 VYYM=VZM+(VYYM-T*VXM)/U VXXM=-(VZZM+VYYM) 2000 G2(1,1)=VXXM*S+OM2*T**2 G2(1,2)=S*VXYM G2(2,1)=G2(1,2) G2(1,3)=S*(VXZM-VXM)-U*T*OM2 G2(3,1)=G2(1,3) G2(2,2)=VYYM*S+OM2 G2(2,3)=S*(VYZM-VYM) G2(3,2)=G2(2,3) G2(3,3)=S*VZZM+U**2*OM2 RETURN END C********************************************************************* SUBROUTINE SETCM(CAPN) C IMPLICIT REAL*8(A-H,O-Z) INTEGER CAPN REAL*4 C,C0 COMMON/SQROOT/DZERO,ROOT(362) COMMON/CM/C20IN,G1(3),G2(3,3),CM3,CM2,CM1,C0,C(32760) C C THIS ROUTINE SETS THE SQUARE ROOT TABLE IN COMMON C CM AND CREATES QUASI-NORMALIZED COEFFICIENTS C FROM NORMALIZED COEFFICIENTS C C C APRIL 1981 CCG C THIS VERSION STORES SQUARE ROOT TABLE IN COMMON SQROOT C DZERO=0.D0 DO 22 I=1,362 22 ROOT(I)=DSQRT(DFLOAT(I)) G1(1)=0.D0 G1(2)=0.D0 G1(3)=0.D0 SMALLC=1.D0 IF(C0.NE.0.D0)SMALLC=1.D0/C0 SQ2=DSQRT(2.D0) DO 200 N=1,CAPN N2=N+N S21=DSQRT(N2+1.D0) K=N**2 C C D IS THE QUASI-NORMALIZATION FACTOR FOR ZONAL TERMS C D=SMALLC*S21 C(K)=C(K)*D C C GG IS THE QUASI-NORMALIZATION FACTOR FOR NON-ZONAL TERMS C GG=D*SQ2 DO 100 J=1,N KJ2=J+J+K C(KJ2-1)=C(KJ2-1)*GG C(KJ2)=C(KJ2)*GG 100 CONTINUE 200 CONTINUE RETURN END C********************************************************************* SUBROUTINE STORC(N,M,CNM,SNM) C C STORE INDIVIDUAL C AND S TERMS C IMPLICIT REAL*8 (A-H,O-Z) REAL*4 C,C0 COMMON/CM/C20IN,G1(3),G2(3,3),CM3,CM2,CM1,C0,C(32760) C C SUM OF THE PREVIOUS NUMBER OF TERMS C J=(N-1)*(N+1) IF(M.EQ.0) GO TO 10 K=M+M C(J+K)=CNM C(J+K+1)=SNM RETURN 10 C(J+1)=CNM RETURN END C********************************************************************* SUBROUTINE LOADCS(NMAX) C IMPLICIT REAL *8(A-H,O-Z) REAL *4 C,C0 COMMON/CM/C20IN,G1(3),G2(3,3),CM3,CM2,CM1,C0,C(32760) 100 READ(12,END=200)N,M,CNM,SNM IF(N.EQ.2.AND.M.EQ.0) C20IN=CNM IF(N.GT.180) GO TO 200 IF(N.GT.NMAX) GO TO 100 IF(N.LT.5) WRITE(6,55) N,M,CNM,SNM 55 FORMAT(1X,4G20.12) CALL STORC(N,M,CNM,SNM) GO TO 100 200 CONTINUE REWIND 12 RETURN END C********************************************************************* SUBROUTINE MODC(N,M,CNM,SNM) C IMPLICIT REAL *8(A-H,O-Z) REAL *4 C,C0 COMMON/CM/C20IN,G1(3),G2(3,3),CM3,CM2,CM1,C0,C(32760) J=(N-1)*(N+1) IF(M.EQ.0) GO TO 10 K = M+M C(J+K) = C(J+K) + CNM C(J+K+1) = C(J+K+1) + SNM RETURN 10 C(J+1) = C(J+1) + CNM RETURN END C********************************************************************* SUBROUTINE NORMAL(GM,DJN,AE,OMEGA,X,Y,Z,U,GAMMA,GRAD) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DJN(20) PSQ = X**2+Y**2 RSQ = PSQ + Z**2 R = DSQRT(RSQ) GMOR = GM/R AOR = AE/R SINE=Z/R PNL1 = SINE PN=1.5D0*SINE**2-0.5D0 F = 1.D0 FP = 0.D0 FPP = 0.D0 AORN = 1.D0 AOR2 = AOR**2 DO 10 N=2,20,2 AORN = AORN*AOR2 TERM=-DJN(N)*AORN*PN F = F + TERM FP = FP -N*TERM/R FPP = FPP+N*(N+1)*TERM/RSQ SAVE =(N+N+1.D0)/(N+1.D0)*SINE*PN-N/(N+1.D0)*PNL1 PNL1 = PN PN = SAVE SAVE = (N+N+3.D0)/(N+2.D0)*SINE*PN-(N+1.D0)/(N+2.D0)*PNL1 PNL1 = PN 10 PN=SAVE U = GMOR*F GAMMA =GMOR*FP-U/R OMEGA2 = OMEGA**2 GRAD=(U+U)/RSQ-2.D0*GMOR*FP/R+GMOR*FPP+OMEGA2*PSQ/RSQ GRAD = DABS(GRAD) U = U+0.5D0*OMEGA2*PSQ GAMMA = DABS(GAMMA+PSQ/R*OMEGA2) RETURN END C********************************************************************* SUBROUTINE REFVAL(A,GM,J2,OMEGA,FLATI) C IMPLICIT REAL*8(A-H,O-Z) REAL*8 J2 PWR=(OMEGA*A)**2*A/GM ESQSV=3.D0*J2+PWR ESQ = ESQSV ITIME = 0 5 ITIME = ITIME + 1 EP2 = ESQ/(1.D0-ESQ) FACT=1.D0/(1.D0-ESQ)**1.5D0 DS=1.D0 TWOQP = 0.D0 DO 10 N = 1, 20 TWOQP=TWOQP+DS*4.D0*N/((N+N+1.D0)*(N+N+3.D0))*FACT DS = -DS FACT = FACT*EP2 10 CONTINUE ESQ=ESQSV+PWR*(4.D0/15.D0/TWOQP-1.D0) FLATI=1.D0/(1.D0-DSQRT(1.D0-ESQ)) TEST = DABS(FLATI-298.257D0) IF(TEST.GT.1.D0) GO TO 100 IF(ITIME.LT.10) GO TO 5 WRITE(6,20) A,GM,J2,OMEGA 20 FORMAT(' REFERENCE VALUES OF NORMAL ELLIPSOID A GM J2 OMEGA'/, * 3X,4G20.12) WRITE(6,30) FLATI 30 FORMAT(' COMPUTED VALUE OF FLATTENING INVERSE ',G20.12) RETURN 100 WRITE(6,110) 110 FORMAT(' SOMETHING WRONG IN REFVAL FLATTENING NOT CONVERGING') STOP END C********************************************************************* SUBROUTINE TRANF(ISWICH,GLAT,ELON,HT,X,Y,Z,AE,FLAT) C IMPLICIT REAL* 8(A-H,O-Z) COMMON/PIDTR/PI,DTR COMMON/TRANCM/TOL,MAXIT FLATFN=(2.D0-FLAT)*FLAT FUNSQ=(1.D0-FLAT)**2 IF(ISWICH.GT.1) GO TO 1000 SPHI = DSIN(GLAT*DTR) G1=AE/DSQRT(1.D0-FLATFN*SPHI**2) G2=G1*FUNSQ+HT G1=G1+HT X=G1*DCOS(GLAT*DTR) Y=X*DSIN(ELON*DTR) X=X*DCOS(ELON*DTR) Z=G2*SPHI RETURN 1000 RSQ=X**2+Y**2 R=DSQRT(RSQ) E=DATAN2(Y,X) IF(E.LT.0.D0) E=E+PI+PI ELON=E/DTR RHO=DSQRT(Z**2+RSQ) SPHI=Z/RHO GLATR=DARSIN(SPHI) HT=RHO-AE*(1.D0-FLAT*SPHI**2) ITER=0 1100 SPHI=DSIN(GLATR) CPHI=DCOS(GLATR) G1=AE/DSQRT(1.D0-FLATFN*SPHI**2) G2=G1*FUNSQ+HT G1 = G1+HT DR=R-G1*CPHI DZ=Z-G2*SPHI DHT=DR*CPHI+DZ*SPHI HT=HT+DHT DLATR = (DZ*CPHI-DR*SPHI)/(AE+HT) GLATR = GLATR+DLATR ITER=ITER+1 IF(ITER.GT.MAXIT) GO TO 1200 IF(DABS(DLATR).GT.TOL) GO TO 1100 IF(DABS(DHT)/(AE+HT).GT.TOL) GO TO 1100 1200 GLAT=GLATR/DTR RETURN END //GO.FT12F001 DD DSN=A.M66774.WENZEL.COEF200.UNFMD,DISP=SHR //GO.SYSIN DD * 67 65 290 288 10 //