PROGRAM IGRFDEC C C 4/91: ~emery/geomag/igrfdec.f C C To compile: f77 -o igrfdec.ob igrfdec.f C To run: igrfdec.ob > igrfdec.out C C IGRFDEC finds the IGRF magnetic field model for a particular geodetic C geographic location, and hence declination, etc. C This version is set to read a master station list and provide the same C list as output along with the computation of the declination angle, C which varies slightly as a function of time. C C RTOD Unit to convert radian to degrees, 45./ATAN(1.). C DTOR Unit to convert degrees to radians ATAN(1.)/45. C RE Radius of the earth, 6371.2 (Km). C REQ Equatorial radius, 6378.165. C COLAT Geographic colatitude of the north pole of the C earth-centered dipole (Deg). C WLON Geographic longitude of the north pole of the C earth-centered dipole (Deg). COMMON /CONST/ RTOD,DTOR,RE,REQ,COLAT,WLON C ICODE Flag set in STOER subprogram. C XI Work array used in subprograms FELDG, SHELG and STOER C H Work array used in subprogram FELDG , SHELG and STOER C DUM1 Work array used in SHELG and STOER. C NAME Not being used any more , used in old version of COFRM C NMAX Order of IGRF model set in COFRM subprogram C TIME Equal to input argument DATE in COFRM routine. C G Array of spherical harmonic coefficients for given C date derived in COFRM subprogram. C DUM2 Array of numbers set in COFRM routine and used in STOER COMMON/MAG/IDUM1,DUM1(307),NAME(4),NMAX,TIME,G(144), + RMIN,RMAX,STEP,STEQ COMMON / LOOPS / IPRNT, JUNIT DIMENSION WK(4) C Put in central dipole pole (Mishin) DATA PLAT/78.5/, PLON/290.5/ DATA RE,REQ /6371.2, 6378.165/ DATA IPRNT,JUNIT/1,6/ character char19*19, char83*83, char36*36 character filein*80, fileout*80 real*4 date PI = 4. * ATAN(1.) DTOR = PI / 180. RTOD = 180. / PI ALT = 130. PI = 3.1415926535898 DTOR = PI / 180. C FIND SOLAR DECLINATION DLE C1 = 23.5 * DTOR DLE=ATAN(TAN(C1)*SIN(2.*PI*(DAY-80.)/365.)) write(6,*) 'Enter master station list file name : ' read(5,'(A80)') filein write(6,*) 'Enter output file name : ' read(5,'(A80)') fileout open (11,file=filein,status='old' ) read (11,"(i3,a19,f7.2)") nsta,char19,DATE write(6,*) 'Enter date (ex : 1991.34)' read(5,*) DATE C IDPOLE = 1,0 if do,not use centered dipole approximation IDPOLE = 0 ICENDP = 1 CALL COFRM (DATE,IDPOLE) open (12,file=fileout,status='new') write (12,"(i3,a19,f7.2,' IDPOLE=',i1,' IGRF (order ',i2, | ') model DECLs (and mag coords if IDPOLE=1)')") nsta,char19, | DATE,IDPOLE,NMAX C Read 2 dummy lines and put into output file read (11,"(a83)") char83 write (12,"(a83)") char83 read (11,"(a83)") char83 write (12,"(a83)") char83 C Read nsta stations, calculate declination for DATE, and output do 100 n=1,nsta read (11,"(a36,2f7.2)") char36,GLAT,GLON CALL FELDG (1,GLAT,GLON,ALT,BN,BE,BD,B) BH = SQRT(BN**2 + BE**2) DINC = ATAN2(BD,BH) * RTOD DECL = ATAN2(BE,BN) * RTOD + 0.0049 IF (ICENDP .EQ. 1) THEN CALL GTM (GLAT*DTOR,GLON*DTOR,PLAT*DTOR,PLON*DTOR,RMLAT,RMLON, | RDIP,RDEC,WK,1) GMLAT = RMLAT / DTOR GMLON = RMLON / DTOR C CALL TSFORM (GLAT, GLON, GMLAT, GMLON, DMAG, 1) write (12,"(a36,2f7.2,1x,2f7.2,2x,f7.2)") char36,GLAT,GLON, | GMLAT,GMLON, DECL GO TO 100 ENDIF write (12,"(a36,2f7.2,17x,f7.2)") char36,GLAT,GLON,DECL 100 CONTINUE C STOP END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Set of routines to calculate the magnetic field and transforms C to various coordinate systems such as the APEX system. Routines C originally from SRI. Were adapted by Harsh Passi at NCAR in 1988 C and then corrected by Cicely Ridley at NCAR in 1989. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE COFRM (DATE,IDPOLE) C***DATE WRITTEN 830415 (YYMMDD) C***LAST REVISION DATE 921005 (YYMMDD) BY B. EMERY, NCAR C***ORIGINAL AUTHOR Wickwar, Vincent B., SRI. INT. C Assign DGRF/IGRF spherical harmonic coefficients for DATE, C yyyy.fraction, into array G. The coefficients are interpolated C from the DGRF/IGRF values for any DATE up to the IGRF initial C year. Coefficients for a later DATE are extrapolated using the IGRF C initial value and the secular change coefficients. A warning C message is issued if DATE is after the last recommended extrapolation. C An DATE input earlier than the first DGRF (EPOCH(1)), results in C an error diagnostic and program STOP. C HISTORY (blame): C The original routine was obtained from SRI. Modifications were to C update the coeff's with DGRF 1980 & IGRF 1985. The new coeff's C were obtained from Eos Vol. 7, No. 24. The common block MAG was C replaced by MAGCOF, thus removing variables not used in subroutine C FELDG. (Jun 1986 - Roy Barnes) C 92 May 5 (ADR) - reactivated MAG instead of MAGCOF. C Modification Apr 1992 (Barnes): Added DGRF 1985 and IGRF 1990 C as described in EOS Vol 73 Number 16 Apr 21 1992. Other changes C were made so future updates should: C (1) Increment NDGY; C (2) Append to EPOCH the next IGRF year; C (3) Append the next DGRF coefficients to G1DIM and H1DIM; and C (4) Replace the IGRF initial values (G0, GT) and rates of C change indices (H0, HT). C Although set up to accomodate second order time derivitives, the C IGRF (GTT, HTT) have been zero. The spherical harmonic coefficients C degree and order is defined by NMAX (currently 10). C C Modification 10/92 (Emery): Added IFRST, IDPOLE, LOOPS C C****************************************************************** C C INPUT C DATE Time (yyyy.fraction) C Added by Cicely Ridley, 1989. C IDPOLE Flag (0,1) to indicate whether to use the full expansion C of B (0), or just the displaced dipole field (1). C OUTPUT C Passed out by the common block MAG see description below. C C***LONG DESCRIPTION C C COMMON Block Used C /LOOPS / IPRNT, JUNIT C /MAG / ICODE,XI(3),H(144),DUM1(2,40),NAME(4),NMAX,TIME,G(144),DUM2(4) C ICODE Flag set in STOER subprogram. C XI Work array used in subprograms FELDG, SHELG and STOER C H Work array used in subprogram FELDG , SHELG and STOER C DUM1 Work array used in SHELG and STOER. C NAME Not being used any more , used in old version of COFRM C NMAX Order of IGRF model set in COFRM subprogram C TIME Equal to input argument DATE in COFRM routine. C G Array of spherical harmonic coefficients for given C date derived in COFRM subprogram. C DUM2 Array of numbers set in COFRM routine and used in STOER C C The following comments accompanied the original code: C****************************************************************** C *****COFRM MUST BE CALLED INITALLY BEFORE FELDG OR SHELG***** C SETS UP COEFFICIENTS (G ARRAY) FOR DATE C C CALLED BY : INITIT C C COEFFICIENTS ARE USED IN STOER ONLY. FIELDG LOADS COEFFS. C******************************************************************** DOUBLE PRECISION F,F0 COMMON / LOOPS / IPRNT, JUNIT C Declare common for communication to FELDG COMMON/MAG/IDUM1,DUM1(307),NAME(4),NMAX,TIME,G(144), + RMIN,RMAX,STEP,STEQ C COMMON/MAGCOF/ NMAX,TIME,G(144) COMMON /CONST/ RTOD,DTOR,RE,REQ,COLAT,WLON C Local declarations PARAMETER (NDGY=5 , NYT = NDGY+1 , NGH = 144*NDGY) C NDGY = Number of DGRF years of sets of coefficients C NYT = Add one for the IGRF set (and point to it). C NGH = Dimension of the equivalenced arrays DIMENSION GYR(12,12,NYT) , HYR(12,12,NYT), EPOCH(NYT) , + G1DIM(NGH) , H1DIM(NGH) , + G0(12,12) , GT(12,12) , GTT(12,12) , + H0(12,12) , HT(12,12) , HTT(12,12) + ,GG(12,12) EQUIVALENCE ( G(1),GG(1,1)) EQUIVALENCE (GYR(1,1,1),G1DIM(1)) , (HYR(1,1,1),H1DIM(1)) , + (GYR(1,1,NYT),G0(1,1)) , (HYR(1,1,NYT),H0(1,1)) SAVE DATEL,GYR,HYR,G0,H0,GT,HT,GTT,HTT,EPOCH,IFRST DATA DATEL /-999./ DATA EPOCH /1965. , 1970. , 1975. , 1980. , 1985. , 1990./ C D_/Dtime2 coefficients are 0 DATA GTT/144*0./,HTT/144*0./ C DGRF g(n,m) for 1965: C The "column" corresponds to "n" and C the "line" corresponds to "m" as indicated in column 6; C e.g., for 1965 g(0,3) = 1297. or g(6,6) = -111. DATA (G1DIM(I),I=1,144) /0., O -30334.,-1662., 1297., 957.,-219., 45., 75., 13., 8.,-2., 2*0., 1 -2119., 2997.,-2038., 804., 358., 61.,-57., 5., 10.,-3., 3*0., 2 1594., 1292., 479., 254., 8., 4., -4., 2., 2., 4*0., 3 856.,-390., -31.,-228.,13.,-14.,-13.,-5., 5*0., 4 252.,-157., 4.,-26., 0., 10.,-2., 6*0., 5 -62., 1., -6., 8., -1., 4., 7*0., 6 -111., 13., -1., -1., 4., 8*0., 7 1., 11., 5., 0., 9*0., 8 4., 1., 2.,10*0., 9 -2., 2.,11*0., A 0.,13*0./ C DGRF g(n,m) for 1970: DATA (G1DIM(I),I=145,288)/ 0., O -30220.,-1781., 1287., 952.,-216., 43., 72., 14., 8.,-3., 2*0., 1 -2068., 3000.,-2091., 800., 359., 64.,-57., 6., 10.,-3., 3*0., 2 1611., 1278., 461., 262., 15., 1., -2., 2., 2., 4*0., 3 838.,-395., -42.,-212.,14.,-13.,-12.,-5., 5*0., 4 234.,-160., 2.,-22., -3., 10.,-1., 6*0., 5 -56., 3., -2., 5., -1., 6., 7*0., 6 -112., 13., 0., 0., 4., 8*0., 7 -2., 11., 3., 1., 9*0., 8 3., 1., 0.,10*0., 9 -1., 3.,11*0., A -1.,13*0./ C DGRF g(n,m) for 1975: DATA (G1DIM(I),I=289,432)/ 0., O -30100.,-1902., 1276., 946.,-218., 45., 71., 14., 7.,-3., 2*0., 1 -2013., 3010.,-2144., 791., 356., 66.,-56., 6., 10.,-3., 3*0., 2 1632., 1260., 438., 264., 28., 1., -1., 2., 2., 4*0., 3 830.,-405., -59.,-198.,16.,-12.,-12.,-5., 5*0., 4 216.,-159., 1.,-14., -8., 10.,-2., 6*0., 5 -49., 6., 0., 4., -1., 5., 7*0., 6 -111., 12., 0., -1., 4., 8*0., 7 -5., 10., 4., 1., 9*0., 8 1., 1., 0.,10*0., 9 -2., 3.,11*0., A -1.,13*0./ C DGRF g(n,m) for 1980: DATA (G1DIM(I),I=433,576)/ 0., O -29992.,-1997., 1281., 938.,-218., 48., 72., 18., 5.,-4., 2*0., 1 -1956., 3027.,-2180., 782., 357., 66.,-59., 6., 10.,-4., 3*0., 2 1663., 1251., 398., 261., 42., 2., 0., 1., 2., 4*0., 3 833.,-419., -74.,-192.,21.,-11.,-12.,-5., 5*0., 4 199.,-162., 4.,-12., -7., 9.,-2., 6*0., 5 -48., 14., 1., 4., -3., 5., 7*0., 6 -108., 11., 3., -1., 3., 8*0., 7 -2., 6., 7., 1., 9*0., 8 -1., 2., 2.,10*0., 9 -5., 3.,11*0., A 0.,13*0./ C DGRF g(n,m) for 1985: DATA (G1DIM(I),I=577,720)/ 0., O -29873.,-2072., 1296., 936.,-214., 53., 74., 21., 5.,-4., 2*0., 1 -1905., 3044.,-2208., 780., 355., 65.,-62., 6., 10.,-4., 3*0., 2 1687., 1247., 361., 253., 51., 3., 0., 1., 3., 4*0., 3 829.,-424., -93.,-185.,24.,-11.,-12.,-5., 5*0., 4 170.,-164., 4., -6., -9., 9.,-2., 6*0., 5 -46., 16., 4., 4., -3., 5., 7*0., 6 -102., 10., 4., -1., 3., 8*0., 7 0., 4., 7., 1., 9*0., 8 -4., 1., 2.,10*0., 9 -5., 3.,11*0., A 0.,13*0./ C DGRF h(n,m) for 1965: DATA (H1DIM(I),I=1,144)/13*0., 1 5776,-2016., -404., 148., 19.,-11.,-61., 7.,-22., 2., 3*0., 2 114., 240.,-269., 128.,100.,-27.,-12., 15., 1., 4*0., 3 -165., 13.,-126., 68., -2., 9., 7., 2., 5*0., 4 -269.,-97.,-32., 6.,-16., -4., 6., 6*0., 5 81., -8., 26., 4., -5.,-4., 7*0., 6 -7.,-23., 24., 10., 0., 8*0., 7 -12., -3., 10.,-2., 9*0., 8 -17., -4., 3.,10*0., 9 1., 0.,11*0., A -6.,13*0./ C DGRF h(n,m) for 1970: DATA (H1DIM(I),I=145,288)/ 13*0., 1 5737.,-2047., -366., 167., 26.,-12.,-70., 7.,-21., 1., 3*0., 2 25., 251.,-266., 139.,100.,-27.,-15., 16., 1., 4*0., 3 -196., 26.,-139., 72., -4., 6., 6., 3., 5*0., 4 -279., -91.,-37., 8.,-17., -4., 4., 6*0., 5 83., -6., 23., 6., -5.,-4., 7*0., 6 1.,-23., 21., 10., 0., 8*0., 7 -11., -6., 11.,-1., 9*0., 8 -16., -2., 3.,10*0., 9 1., 1.,11*0., A -4.,13*0./ C DGRF h(n,m) for 1975: DATA (H1DIM(I),I=289,432)/ 13*0., 1 5675.,-2067., -333., 191., 31.,-13.,-77., 6.,-21., 1., 3*0., 2 -68., 262.,-265., 148., 99.,-26.,-16., 16., 1., 4*0., 3 -223., 39.,-152., 75., -5., 4., 7., 3., 5*0., 4 -288., -83.,-41., 10.,-19., -4., 4., 6*0., 5 88., -4., 22., 6., -5.,-4., 7*0., 6 11.,-23., 18., 10.,-1., 8*0., 7 -12.,-10., 11.,-1., 9*0., 8 -17., -3., 3.,10*0., 9 1., 1.,11*0., A -5.,13*0./ C DGRF h(n,m) for 1980: DATA (H1DIM(I),I=433,576)/ 13*0., 1 5604.,-2129., -336., 212., 46.,-15.,-82., 7.,-21., 1., 3*0., 2 -200., 271.,-257., 150., 93.,-27.,-18., 16., 0., 4*0., 3 -252., 53.,-151., 71., -5., 4., 9., 3., 5*0., 4 -297., -78.,-43., 16.,-22., -5., 6., 6*0., 5 92., -2., 18., 9., -6.,-4., 7*0., 6 17.,-23., 16., 9., 0., 8*0., 7 -10.,-13., 10.,-1., 9*0., 8 -15., -6., 4.,10*0., 9 2., 0.,11*0., A -6.,13*0./ C DGRF h(n,m) for 1985: DATA (H1DIM(I),I=577,720)/ 13*0., 1 5500.,-2197., -310., 232., 47.,-16.,-83., 8.,-21., 1., 3*0., 2 -306., 284.,-249., 150., 88.,-27.,-19., 15., 0., 4*0., 3 -297., 69.,-154., 69., -2., 5., 9., 3., 5*0., 4 -297., -75.,-48., 20.,-23., -6., 6., 6*0., 5 95., -1., 17., 11., -6.,-4., 7*0., 6 21.,-23., 14., 9., 0., 8*0., 7 -7.,-15., 9.,-1., 9*0., 8 -11., -7., 4.,10*0., 9 2., 0.,11*0., A -6.,13*0./ C Initial coefficients g0 (IGRF for 1990): DATA G0/0., O -29775.,-2136., 1315., 939.,-211., 61., 77., 22., 4.,-4., 2*0., 1 -1851., 3058.,-2240., 782., 353., 64.,-64., 5., 10.,-4., 3*0., 2 1693., 1246., 324., 244., 60., 4., -1., 1., 2., 4*0., 3 807.,-423.,-111.,-178.,28.,-11.,-12.,-5., 5*0., 4 142.,-166., 2., 1.,-12., 9.,-2., 6*0., 5 -37., 17., 6., 4., -4., 4., 7*0., 6 -96., 10., 4., -1., 3., 8*0., 7 0., 3., 7., 1., 9*0., 8 -6., 2., 2.,10*0., 9 -6., 3.,11*0., A 0.,13*0./ C D_/Dtime coefficients gt (IGRF for 1990-1995): DATA GT/0., O 18.0, -12.9, 3.3, 0.5, 0.6, 1.3, 0.6, 0.2, 0.0,0.0, 2*0., 1 10.6, 2.4, -6.7, 0.6, -0.1,-0.2,-0.5,-0.7, 0.0,0.0, 3*0., 2 0.0, 0.1, -7.0, -1.6, 1.8,-0.3,-0.2, 0.0,0.0, 4*0., 3 -5.9, 0.5, -3.1, 1.3, 0.6, 0.1, 0.0,0.0, 5*0., 4 -5.5, -0.1,-0.2, 1.6,-1.1, 0.0,0.0, 6*0., 5 2.3, 0.1, 0.2, 0.0, 0.0,0.0, 7*0., 6 1.2, 0.2,-0.1, 0.0,0.0, 8*0., 7 0.3,-0.5, 0.0,0.0, 9*0., 8 -0.6, 0.0,0.0,10*0., 9 0.0,0.0,11*0., A 0.0,13*0./ C Initial coefficients h0 (IGRF for 1990-1995): DATA H0/13*0., 1 5411.,-2278., -287., 248., 47.,-16.,-81., 10.,-21., 1., 3*0., 2 -380., 293.,-240., 153., 83.,-27.,-20., 15., 0., 4*0., 3 -348., 87.,-154., 68., 1., 7., 10., 3., 5*0., 4 -299., -69.,-52., 20.,-22., -6., 6., 6*0., 5 98., 2., 16., 12., -6.,-4., 7*0., 6 27.,-23., 11., 9., 0., 8*0., 7 -5.,-16., 9.,-1., 9*0., 8 -11., -7., 4.,10*0., 9 2., 0.,11*0., A -6.,13*0./ C D_/Dtime coefficients ht (IGRF for 1990-1995): DATA HT/13*0., 1 -16.1, -15.8, 4.4, 2.6, -0.1, 0.2, 0.6, 0.5, 0.0,0.0, 3*0., 2 -13.8, 1.6, 1.8, 0.5,-1.3, 0.2,-0.2, 0.0,0.0, 4*0., 3 -10.6, 3.1, 0.4, 0.0, 0.8, 0.3, 0.0,0.0, 5*0., 4 -1.4, 1.7,-0.9,-0.5, 0.3, 0.0,0.0, 6*0., 5 0.4, 0.5,-0.2, 0.4, 0.0,0.0, 7*0., 6 1.2, 0.0,-0.5, 0.0,0.0, 8*0., 7 0.0,-0.3, 0.0,0.0, 9*0., 8 0.6, 0.0,0.0,10*0., 9 0.0,0.0,11*0., A 0.0,13*0./ DATA IFRST/0/ IFRST = IFRST + 1 IF (IFRST .EQ. 1) |WRITE (6,"(1X,'COFRM includes 1990 coefs: Last update 10/92')") C Set constant variables for subprogram SHELG C RMIN=0.05 C RMAX=1.01 C STEP=0.20 C STEQ=0.03 C Do not need to load new coefficients IF (DATE .EQ. DATEL) GO TO 990 DATEL = DATE C Trap out of range date: IF (DATE .LT. EPOCH(1)) GO TO 9100 CC IF (DATE .GT. EPOCH(NYT)+5.) WRITE(6,9200) DATE IF (DATE .GT. EPOCH(NYT)+6.) WRITE(6,9200) DATE DO 100 I=1,NYT IF (DATE .LT. EPOCH(I)) GO TO 110 IY = I 100 CONTINUE 110 CONTINUE NMAX = 10 WRITE( JUNIT, "(' Order of IGRF model NMAX= 'I3,' for DATE=', | F7.2,' IDPOLE =',I2)") NMAX,DATE,IDPOLE TIME = DATE IS = 0 T = TIME-EPOCH(IY) G(1) = 0.0 I = 2 F0 = 1.0D-5 IF (IS) 200,190,200 190 F0 = -F0 200 DO 600 N=1,NMAX FN = N F0 = F0*FN*FN/(4.0*FN-2.0) IF (IS) 300,290,300 290 F0 = F0*(2.0*FN-1.0)/FN 300 F = F0*0.5 IF (IS) 400,390,400 390 F = F*SQRT(2.0) 400 NN = N+1 MM = 1 IF (IY .EQ. NYT) THEN C Extrapolate coefficients G(I) = ((GTT(NN,MM)*T + GT(NN,MM))*T + G0(NN,MM)) * F0 ELSE C Interpolate coefficients G(I) = (GYR(NN,MM,IY) + + T/5.0 * (GYR(NN,MM,IY+1)-GYR(NN,MM,IY))) * F0 ENDIF I = I+1 DO 600 M=1,N TMP1 = N+M TMP2 = N-M+1 F = F*TMP1/TMP2 IF (IS) 500,490,500 490 F = F*SQRT(TMP2/TMP1) 500 NN = N+1 MM = M+1 I1 = I+1 IF (IY .EQ. NYT) THEN C Extrapolate coefficients G(I) = ((GTT(NN,MM)*T + GT(NN,MM))*T + G0(NN,MM)) * F G(I1) = ((HTT(NN,MM)*T + HT(NN,MM))*T + H0(NN,MM)) * F ELSE C Interpolate coefficients G(I) = (GYR(NN,MM,IY) + + T/5.0 * (GYR(NN,MM,IY+1)-GYR(NN,MM,IY))) * F G(I1) = (HYR(NN,MM,IY) + + T/5.0 * (HYR(NN,MM,IY+1)-HYR(NN,MM,IY))) * F ENDIF 600 I = I+2 IF(IDPOLE.EQ.1)THEN C **** SET ALL BUT DIPOLE COEFFICIENTS TO ZERO G(1) = 0. DO 10 I = 5,144 G(I) = 0. 10 CONTINUE ENDIF 990 CONTINUE C COMPUTE GEOGRAPHIC COLATITUDE AND LONGITUDE OF THE NORTH POLE OF C EARTH CENTERED DIPOLE C********************** C COLAT = ACOS(GG(2,1)/SQRT(GG(2,1)**2+GG(3,1)**2+GG(4,1)**2)) COLAT = ACOS( G(2 )/SQRT( G(2 )**2+ G(3 )**2+ G(4 )**2)) 1 *RTOD C WLON = ATAN2(GG(4,1),GG(3,1))*RTOD WLON = ATAN2( G(4 ), G(3 ))*RTOD PRINT *, ' MAGNETIC POLE ', COLAT,WLON IF (GG(2,1)+GG(4,1) .NE. GG(4,1)) THEN COLATRV = ATAN2(SQRT(GG(3,1)**2+GG(4,1)**2),GG(2,1)**2)*RTOD C ELSE COLATRV = 90. END IF IF (GG(3,1)+GG(4,1) .NE. GG(4,1)) THEN WLONRV = ATAN2(GG(4,1),GG(3,1))*RTOD C ELSE WLONRV = 90. END IF C********************** PRINT *, ' MAGNETIC POLE ', COLATRV,WLONRV RETURN C Error trap diagnostics: 9100 WRITE(6,'('' '',/, +'' COFRM: DATE'',F9.3,'' preceeds DGRF coefficients'', + '' presently coded.'')') DATE STOP 'mor cod' 9200 FORMAT(' ',/, +' COFRM: DATE',F9.3,' is after the maximum', + ' recommended for extrapolation.') END SUBROUTINE OLDCOFRM (DATE,IDPOLE) SAVE C***BEGIN PROLOGUE COFRM C***DATE WRITTEN 830415 (YYMMDD) C***REVISION DATE 880201 (YYMMDD) C***AUTHOR Wickwar, Vincent B., SRI. INT. C***PURPOSE This assigns DGRF/IGRF spherical harmonic coefficients C for date TIME (yyyy.fraction) into array G. The coeff's C are interpolated from the 5 year DGRF/IGRF values, or C extrapolated if TIME is more recent than 1985.0 . The first C year of DGRF coeff's is 1965; if TIME is earlier than this C an error message is generated. C***DESCRIPTION C C The origonal routine was obtained from SRI. Modifications C were to update the coeff's with DGRF 1980 & IGRF 1985. C The new coeff's were obtained from Eos Vol. 7, No. 24, C 17 Jun 1986. (April 1987 - Roy Barnes) C COFRM must be called before FELDG OR SHELG subprograms. C It sets up coefficients ( G array ) for date. C C INPUT C DATE Time (yyyy.fraction) C Added by Cicely Ridley, 1989. C IDPOLE Flag (0,1) to indicate whether to use the full expansion C of B (0), or just the displaced dipole field (1). C OUTPUT C Passed out by the common block MAG see description below. C C***LONG DESCRIPTION C C COMMON Block Used C /LOOPS / IPRNT, JUNIT C /MAG / ICODE,XI(3),H(144),DUM1(2,40),NAME(4),NMAX,TIME,G(144),DUM2(4) C ICODE Flag set in STOER subprogram. C XI Work array used in subprograms FELDG, SHELG and STOER C H Work array used in subprogram FELDG , SHELG and STOER C DUM1 Work array used in SHELG and STOER. C NAME Not being used any more , used in old version of COFRM C NMAX Order of IGRF model set in COFRM subprogram C TIME Equal to input argument DATE in COFRM routine. C G Array of spherical harmonic coefficients for given C date derived in COFRM subprogram. C DUM2 Array of numbers set in COFRM routine and used in STOER C C***REFERENCES reference 1 C continuation of reference 1 C C***ROUTINES CALLED None C***COMMON BLOCKS MAG C***END PROLOGUE COFRM DOUBLE PRECISION F,F0 COMMON / LOOPS / IPRNT, JUNIT COMMON/MAG/IDUM1,DUM1(307),NAME(4),NMAX,TIME,G(144), + RMIN,RMAX,STEP,STEQ COMMON /CONST/ RTOD,DTOR,RE,REQ,COLAT,WLON DIMENSION G0(12,12),GT(12,12),GTT(12,12),H0(12,12),HT(12,12), + HTT(12,12) DIMENSION GYR(12,12,5),HYR(12,12,5), + G1DIM(576), H1DIM(576),GG(12,12) EQUIVALENCE ( G(1),GG(1,1)) EQUIVALENCE ( GYR(1,1,1), G1DIM(1) ),( HYR(1,1,1), H1DIM(1) ), + (G0(1,1),GYR(1,1,5)) , (H0(1,1),HYR(1,1,5)) SAVE DATEL,GYR,HYR,GO,HO,GT,HT,GTT,HTT DATA DATEL/-999./ C D_/Dtime2 coefficients are 0 (for extrapolation T>1985.0) DATA GTT/144*0./,HTT/144*0./ C DGRF g(n,m) for 1965: DATA (G1DIM(I),I=1,144) /0., 1 -30334.,-1662., 1297., 957.,-219., 45., 75., 13., 8.,-2., 2*0., 1 -2119., 2997.,-2038., 804., 358., 61.,-57., 5., 10.,-3., 3*0., 2 1594., 1292., 479., 254., 8., 4., -4., 2., 2., 4*0., 3 856.,-390., -31.,-228.,13.,-14.,-13.,-5., 5*0., 4 252.,-157., 4.,-26., -0., 10.,-2., 6*0., 5 -62., 1., -6., 8., -1., 4., 7*0., 6 -111., 13., -1., -1., 4., 8*0., 7 1., 11., 5., 0., 9*0., 8 4., 1., 2.,10*0., 9 -2., 2.,11*0., A 0.,13*0./ C DGRF g(n,m) for 1970: DATA (G1DIM(I),I=145,288)/ 0., 1 -30220.,-1781., 1287., 952.,-216., 43., 72., 14., 8.,-3., 2*0., 1 -2068., 3000.,-2091., 800., 359., 64.,-57., 6., 10.,-3., 3*0., 2 1611., 1278., 461., 262., 15., 1., -2., 2., 2., 4*0., 3 838.,-395., -42.,-212.,14.,-13.,-12.,-5., 5*0., 4 234.,-160., 2.,-22., -3., 10.,-1., 6*0., 5 -56., 3., -2., 5., -1., 6., 7*0., 6 -112., 13., 0., 0., 4., 8*0., 7 -2., 11., 3., 1., 9*0., 8 3., 1., 0.,10*0., 9 -1., 3.,11*0., A -1.,13*0./ C DGRF g(n,m) for 1975: DATA (G1DIM(I),I=289,432)/ 0., 1 -30100.,-1902., 1276., 946.,-218., 45., 71., 14., 7.,-3., 2*0., 1 -2013., 3010.,-2144., 791., 356., 66.,-56., 6., 10.,-3., 3*0., 2 1632., 1260., 438., 264., 28., 1., -1., 2., 2., 4*0., 3 830.,-405., -59.,-192.,16.,-12.,-12.,-5., 5*0., 4 216.,-159., 1.,-14., -8., 10.,-2., 6*0., 5 -49., 6., 0., 4., -1., 5., 7*0., 6 -111., 12., 0., -1., 4., 8*0., 7 -5., 10., 4., 1., 9*0., 8 1., 1., 0.,10*0., 9 -2., 3.,11*0., A -1.,13*0./ C DGRF g(n,m) for 1980: DATA (G1DIM(I),I=433,576)/ 0., 1 -29992.,-1997., 1281., 938.,-218., 48., 72., 18., 5.,-4., 2*0., 1 -1956., 3027.,-2180., 782., 357., 66.,-59., 6., 10.,-4., 3*0., 2 1663., 1251., 398., 261., 42., 2., 0., 1., 2., 4*0., 3 833.,-419., -74.,-198.,21.,-11.,-12.,-5., 5*0., 4 199.,-162., 4.,-12., -7., 9.,-2., 6*0., 5 -48., 14., 1., 4., -3., 5., 7*0., 6 -108., 11., 3., -1., 3., 8*0., 7 -2., 6., 7., 1., 9*0., 8 -1., 2., 2.,10*0., 9 -5., 3.,11*0., A 0.,13*0./ C DGRF h(n,m) for 1965: DATA (H1DIM(I),I=1,144)/13*0., 1 5776,-2016., -404., 148., 19.,-11.,-61., 7.,-22., 2., 3*0., 2 114., 240.,-269., 128.,100.,-27.,-12., 15., 1., 4*0., 3 -165., 13.,-126., 68., -2., 9., 7., 2., 5*0., 4 -269.,-97.,-32., 6.,-16., -4., 6., 6*0., 5 81., -8., 26., 4., -5.,-4., 7*0., 6 -7.,-23., 24., 10., 0., 8*0., 7 -12., -3., 10.,-2., 9*0., 8 -17., -4., 3.,10*0., 9 1., 0.,11*0., A -6.,13*0./ C DGRF h(n,m) for 1970: DATA (H1DIM(I),I=145,288)/ 13*0., 1 5737.,-2047., -366., 167., 26.,-12.,-70., 7.,-21., 1., 3*0., 2 25., 251.,-266., 139.,100.,-27.,-15., 16., 1., 4*0., 3 -196., 26.,-139., 72., -4., 6., 6., 3., 5*0., 4 -279., -91.,-37., 8.,-17., -4., 4., 6*0., 5 83., -6., 23., 6., -5.,-4., 7*0., 6 1.,-23., 21., 10., 0., 8*0., 7 -11., -6., 11.,-1., 9*0., 8 -16., -2., 3.,10*0., 9 1., 1.,11*0., A -4.,13*0./ C DGRF h(n,m) for 1975: DATA (H1DIM(I),I=289,432)/ 13*0., 1 5675.,-2067., -333., 191., 31.,-13.,-77., 6.,-21., 1., 3*0., 2 -68., 262.,-265., 148., 99.,-26.,-16., 16., 1., 4*0., 3 -223., 39.,-152., 75., -5., 4., 7., 3., 5*0., 4 -288., -83.,-41., 10.,-19., -4., 4., 6*0., 5 88., -4., 22., 6., -5.,-4., 7*0., 6 11.,-23., 18., 10.,-1., 8*0., 7 -12.,-10., 11.,-1., 9*0., 8 -17., -3., 3.,10*0., 9 1., 1.,11*0., A -5.,13*0./ C DGRF h(n,m) for 1980: DATA (H1DIM(I),I=433,576)/ 13*0., 1 5604.,-2129., -336., 212., 46.,-15.,-82., 7.,-21., 1., 3*0., 2 -200., 271.,-257., 150., 93.,-27.,-18., 16., 0., 4*0., 3 -252., 53.,-151., 71., -5., 4., 9., 3., 5*0., 4 -297., -78.,-43., 16.,-22., -5., 6., 6*0., 5 92., -2., 18., 9., -6.,-4., 7*0., 6 17.,-23., 16., 9., 0., 8*0., 7 -10.,-13., 10.,-1., 9*0., 8 -15., -6., 4.,10*0., 9 2., 0.,11*0., A -6.,13*0./ C Initial coefficients g0 (for extrapolation T>1985.0): DATA G0/0., 1 -29877.,-2073., 1300., 937.,-215., 52., 75., 21., 5.,-4., 2*0., 1 -1903., 3045.,-2208., 780., 356., 65.,-61., 6., 10.,-4., 3*0., 2 1691., 1244., 363., 253., 50., 2., 0., 1., 2., 4*0., 3 835.,-426., -94.,-186.,24.,-11.,-12.,-5., 5*0., 4 169.,-161., 4., -6., -9., 9.,-2., 6*0., 5 -48., 17., 4., 2., -3., 5., 7*0., 6 -102., 9., 4., -1., 3., 8*0., 7 0., 4., 7., 1., 9*0., 8 -6., 2., 2.,10*0., 9 -5., 3.,11*0., A 0.,13*0./ C D_/Dtime coefficients gt (for extrapolation T>1985.0): DATA GT/0., 1 23.2, -13.7, 5.1, 0.1, 1.3, 1.4, 0.2, 0.7, 0.0,0.0, 2*0., 1 10.0, 3.4, -4.6, -0.6, 0.1,-0.3,-0.6, 0.0, 0.0,0.0, 3*0., 2 7.0, -0.6, -7.8, -1.5, 1.7,-0.5, 0.3, 0.0,0.0, 4*0., 3 0.1, -1.4, -3.2, 0.6, 0.8, 0.4, 0.0,0.0, 5*0., 4 -6.8, 0.1, 0.0, 1.0,-0.3, 0.0,0.0, 6*0., 5 -0.1, 0.9, 0.4,-0.3, 0.0,0.0, 7*0., 6 1.2,-0.5, 0.1, 0.0,0.0, 8*0., 7 -0.1,-0.5, 0.0,0.0, 9*0., 8 -0.8, 0.0,0.0,10*0., 9 0.0,0.0,11*0., A 0.0,13*0./ C Initial coefficients h0 (for extrapolation T>1985.0): DATA H0/13*0., 1 5497.,-2191., -312., 233., 47.,-16.,-82., 7.,-21., 1., 3*0., 2 -309., 284.,-250., 148., 90.,-26.,-21., 16., 0., 4*0., 3 -296., 68.,-155., 69., -1., 5., 9., 3., 5*0., 4 -298., -75.,-50., 23.,-25., -5., 6., 6*0., 5 95., -4., 17., 11., -6.,-4., 7*0., 6 20.,-21., 12., 9., 0., 8*0., 7 -6.,-16., 10.,-1., 9*0., 8 -10., -6., 4.,10*0., 9 2., 0.,11*0., A -6.,13*0./ C D_/Dtime coefficients ht (for extrapolation T>1985.0): DATA HT/13*0., 1 -24.5, -11.5, 5.3, 3.8, 0.1,-0.4, 0.2, 0.1, 0.0,0.0, 3*0., 2 -20.2, 2.3, 2.2, -0.2,-1.1, 1.0,-1.0, 0.0,0.0, 4*0., 3 -10.8, 2.5, -0.1,-0.8, 1.1, 0.1, 0.0,0.0, 5*0., 4 0.9, 0.6,-2.3, 1.9,-0.8, 0.0,0.0, 6*0., 5 0.0,-0.5, 0.3, 0.2, 0.0,0.0, 7*0., 6 -0.1, 0.2,-0.8, 0.0,0.0, 8*0., 7 0.9,-0.1, 0.0,0.0, 9*0., 8 1.3, 0.0,0.0,10*0., 9 0.0,0.0,11*0., A 0.0,13*0./ C C SET CONSTANT VARIABLES FOR SUBPROGRAM SHELG C RMIN=0.05 RMAX=1.01 STEP=0.20 STEQ=0.03 C C Trap out of range date: IF(DATE .LT. 1965.) GO TO 9100 IF(DATE .GT. 1990.)WRITE(6,9200) C Do not need to load new coefficients IF (DATE .EQ. DATEL) RETURN DATEL=DATE NMAX=10 WRITE( JUNIT, "(' Order of IGRF model NMAX= 'I3,' for DATE=', | F7.2)") NMAX,DATE TIME=DATE IF (DATE .LT. 1970.0)THEN INT=1 EPOCH=1965.0 ELSEIF(DATE .LT. 1975.0)THEN INT=2 EPOCH=1970.0 ELSEIF(DATE .LT. 1980.0)THEN INT=3 EPOCH=1975.0 ELSEIF(DATE .LT. 1985.0)THEN INT=4 EPOCH=1980.0 ELSE INT=5 EPOCH=1985.0 ENDIF IS=0 T=TIME-EPOCH G(1)=0.0 I=2 F0=1.0D-5 IF (IS) 5,4,5 4 F0=-F0 5 DO 9 N=1,NMAX FN=N F0=F0*FN*FN/(4.0*FN-2.0) IF (IS) 502,501,502 501 F0=F0*(2.0*FN-1.0)/FN 502 F=F0*0.5 IF (IS) 6,503,6 503 F=F*SQRT(2.0) 6 NN=N+1 MM=1 IF (INT.EQ.5)THEN C Extrapolate coefficients G(I)=((GTT(NN,MM)*T+GT(NN,MM))*T+G0(NN,MM))*F0 ELSE C Interpolate coefficients G(I)=(GYR(NN,MM,INT)+T/5.0*(GYR(NN,MM,INT+1)-GYR(NN,MM,INT))) + *F0 ENDIF I=I+1 DO 9 M=1,N TMP1=N+M TMP2=N-M+1 F=F*TMP1/TMP2 IF (IS) 602,601,602 601 F=F*SQRT(TMP2/TMP1) 602 NN=N+1 MM=M+1 I1=I+1 IF (INT .EQ. 5)THEN C Extrapolate coefficients G(I) =((GTT(NN,MM)*T+GT(NN,MM))*T+G0(NN,MM))*F G(I1)=((HTT(NN,MM)*T+HT(NN,MM))*T+H0(NN,MM))*F ELSE C Interpolate coefficients G(I) =(GYR(NN,MM,INT)+T/5.0*(GYR(NN,MM,INT+1)-GYR(NN,MM,INT)))*F G(I1)=(HYR(NN,MM,INT)+T/5.0*(HYR(NN,MM,INT+1)-HYR(NN,MM,INT)))*F ENDIF 9 I=I+2 IF(IDPOLE.EQ.1)THEN C **** SET ALL BUT DIPOLE COEFFICIENTS TO ZERO G(1) = 0. DO 10 I = 5,144 G(I) = 0. 10 CONTINUE ENDIF C COMPUTE GEOGRAPHIC COLATITUDE AND LONGITUDE OF THE NORTH POLE OF C EARTH CENTERED DIPOLE C********************** COLAT = ACOS(GG(2,1)/SQRT(GG(2,1)**2+GG(3,1)**2+GG(4,1)**2)) 1 *RTOD WLON = ATAN2(GG(4,1),GG(3,1))*RTOD PRINT *, ' MAGNETIC POLE ', COLAT,WLON IF (GG(2,1)+GG(4,1) .NE. GG(4,1)) THEN COLATRV = ATAN2(SQRT(GG(3,1)**2+GG(4,1)**2),GG(2,1)**2)*RTOD C ELSE COLATRV = 90. END IF IF (GG(3,1)+GG(4,1) .NE. GG(4,1)) THEN WLONRV = ATAN2(GG(4,1),GG(3,1))*RTOD C ELSE WLONRV = 90. END IF C********************** PRINT *, ' MAGNETIC POLE ', COLATRV,WLONRV RETURN C Error trap diagnostics: 9100 WRITE(6,"('0Requested year ',F8.3,' preceeds IGRF/DGRF ', + 'coefficients presently coded. See',/, + ' subroutine COFRM for modifications.')")DATE STOP 'mor cod' 9200 FORMAT('0Warning: IGRF model coefficients are extrapolated from', + ' 1985.') END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE FELDG(IENTY,GLAT,GLON,ALT,BNRTH,BEAST,BDOWN,BABS) SAVE C***BEGIN PROLOGUE FELDG C***DATE WRITTEN 830415 (YYMMDD) C***REVISION DATE 880201 (YYMMDD) C***AUTHOR Wickwar, Vincent B., SRI INT. C***PURPOSE It computes the DGRF/IGRF coeffs. at the point GLAT, C GLON,ALT. COFRM must be called to establish coeffs C G (array) for correct date prior to calling FELDG. C***DESCRIPTION C C INPUT C IENTY If IENTY=1, (GLAT,GLON,ALT) Geodetic coordinates C If IENTY=2, (GLAT,GLON,ALT) Earth centered cartesian C coordinates C If IENTY=3, Entry point used for L computation. C INPUT IF IENTY=1 C GLAT LATITUDE OF POINT (DEG) C GLON LONGITUDE (EAST=+) OF POINT (DEG) C ALT HT OF POINT (KM) C OUTPUT IF IENTY=1 C BNRTH NORTH COMPONENT OF FIELD VECTOR (Gauss) C BEAST EAST COMPONENT OF FIELD VECTOR (Gauss) C BDOWN DOWNWARD COMPONENT OF FIELD VECTOR (Gauss) C BABS MAGNITUDE OF FIELD VECTOR (Gauss) C INPUT IF IENTY=2 C GLAT X coordinate (in units of earth radii 6371.2 km ) C GLON Y coordinate (in units of earth radii 6371.2 km ) C ALT Z coordinate (in units of earth radii 6371.2 km ) C OUTPUT IF IENTY=2 C BNRTH X COMPONENT OF FIELD VECTOR (Gauss) C BEAST Y COMPONENT OF FIELD VECTOR (Gauss) C BDOWN Z COMPONENT OF FIELD VECTOR (Gauss) C BABS MAGNITUDE OF FIELD VECTOR (Gauss) C C INPUT IF IENTY=3 C GLAT .Dummy arg C GLON Dummy arg C ALT Dummy arg C OUTPUT IF IENTY=3 C Passed through the common block MAG. C See the description below. C***LONG DESCRIPTION C C COMMON Block Used C /MAG / ICODE,XI(3),H(144),DUM1(2,40),NAME(4),NMAX,TIME,G(144),DUM2(4) C ICODE Flag set in STOER subprogram. C XI Work array used in subprograms FELDG, SHELG and STOER C Also used for communication. C H Work array used in subprogram FELDG , SHELG and STOER C Also used for communication. C DUM1 Work array used in SHELG and STOER. C NAME Not being used any more , used in old version of COFRM C NMAX Order of IGRF model set in COFRM subprogram C TIME Equal to input argument DATE in COFRM routine. C G Array of spherical harmonic coefficients for given C date derived in COFRM subprogram. C DUM2 Array of numbers set in COFRM routine and used in STOER C C /CONST/ RTOD,DTOR,RE,REQ,COLAT,WLON C RTOD Unit to convert radian to degrees, 45./ATAN(1.). C DTOR Unit to convert degrees to radians ATAN(1.)/45. C RE Radius of the earth, 6371.2 (Km). C REQ Equatorial radius, 6378.165. C COLAT Geographic colatitude of the north pole of the C earth-centered dipole (Deg). C WLON Geographic longitude of the north pole of the C earth-centered dipole (Deg). C***REFERENCES reference 1 C continuation of reference 1 C C***ROUTINES CALLED None C***COMMON BLOCKS MAG C***END PROLOGUE FELDG COMMON/MAG/ICODE,XI(3),H(144),DUM1(8,20),NAME(4), + NMAX,TIME,G(144),DUM2(4) COMMON / CONST / RTOD,DTOR,RE,REQ,COLAT,WLON DIMENSION V(3),B(3) GO TO(100,200,300), IENTY 100 CONTINUE IS=1 RLAT = GLAT*DTOR CT = SIN(RLAT) ST = COS(RLAT) D = SQRT(40680925.E0-272336.E0*CT**2) RLON = GLON*DTOR CP = COS(RLON) SP = SIN(RLON) ZZZ = (ALT+40408589.E0/D)*CT/RE RHO = (ALT+40680925.E0/D)*ST/RE XXX = RHO*CP YYY = RHO*SP GO TO 10 C ENTRY FELDC 200 CONTINUE C*****ENTRY POINT FELDC TO BE USED WITH CARTESIAN CO-ORDINATES IS = 2 XXX = GLAT V(1) = GLAT YYY = GLON V(2) = GLON ZZZ = ALT V(3) = ALT 10 RQ = 1./(XXX**2+YYY**2+ZZZ**2) XI(1) = XXX*RQ XI(2) = YYY*RQ XI(3) = ZZZ*RQ GO TO 20 C ENTRY FELDI 300 CONTINUE C*****ENTRY POINT FELDI USED FOR L COMPUTATION IS = 3 20 IHMAX = NMAX*NMAX+1 LAST = IHMAX+NMAX+NMAX IMAX = NMAX+NMAX-1 DO 8 I=IHMAX,LAST 8 H(I) = G(I) MK = 3 IF ( IMAX .EQ. 1) MK=1 DO 6 K=1,MK,2 I = IMAX IH = IHMAX 1 IL = IH-I F = 2./FLOAT(I-K+2) X = XI(1)*F Y = XI(2)*F Z = XI(3)*(F+F) I = I-2 IF(I-1)5,4,2 2 DO 3 M=3,I,2 IHM = IH+M ILM = IL+M H(ILM+1) = G(ILM+1)+ Z*H(IHM+1) + X*(H(IHM+3)-H(IHM-1)) + -Y*(H(IHM+2)+H(IHM-2)) 3 H(ILM) = G(ILM) + Z*H(IHM) + X*(H(IHM+2)-H(IHM-2)) + +Y*(H(IHM+3)+H(IHM-1)) 4 H(IL+2) = G(IL+2) + Z*H(IH+2) + X*H(IH+4) - Y*(H(IH+3)+H(IH)) H(IL+1) = G(IL+1) + Z*H(IH+1) + Y*H(IH+4) + X*(H(IH+3)-H(IH)) 5 H(IL) = G(IL) + Z*H(IH) + 2.*(X*H(IH+1)+Y*H(IH+2)) IH = IL IF (I-K) 6,1,1 6 CONTINUE IF(IS .EQ. 3)RETURN C Have completed the L shell calculation S = .5*H(1)+2.*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2)) T = (RQ+RQ)*SQRT(RQ) BXXX = T*(H(3)-S*XXX) BYYY = T*(H(4)-S*YYY) BZZZ = T*(H(2)-S*ZZZ) BABS = SQRT(BXXX**2+BYYY**2+BZZZ**2) C Convert from gauss to nT C BABS = BABS * 1.E5 IF(IS .EQ. 1)THEN C Convert back to geodetic BEAST = BYYY*CP-BXXX*SP BRHO = BYYY*SP+BXXX*CP BNRTH = BZZZ*ST-BRHO*CT BDOWN = -BZZZ*CT-BRHO*ST C Leave in earth centered cartesian ELSEIF(IS .EQ. 2)THEN BNRTH = BXXX B(1) = BXXX BEAST = BYYY B(2) = BYYY BDOWN = BZZZ B(3) = BZZZ ENDIF RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE TSFORM(CLAT, CLON, GMLAT, GMLON, DMAG, I) SAVE C***BEGIN PROLOGUE TSFORM C***DATE WRITTEN 731029 (YYMMDD) C***REVISION DATE 880201 (YYMMDD) C***AUTHOR CLARK, W, N.O.A.A. C***PURPOSE Converts geocentric to geomagnetic coordinates, or C vice-versa. It is used to find apex longitude. C i.e. geomagnetic longitude. C***DESCRIPTION C C C INPUT C I If I=1, convert from geocentric to geomagnetic point C If I=2, convert from geomagnetic to geocentric point C INPUT IF I .EQ. 1 C CLAT Geocentric latitude (Deg) C CLON Geocentric longitude (Deg) C OUTPUT IF I .EQ. 1 C GMLAT Geomagnetic latitude (Deg) C GMLON Geomagnetic longitude (Deg) C SIND Constant which may be used in conversion of C X and Y field components from geocentric C to geomagnetic directions. C =(cos(colat)-sin(gmlat)sin(clat))/cos(gmlat)cos(clat) C where colat is the geocentric colat of the north pole C of the earth-centered dipole. C COSD Constant which may be used in conversion of C X and Y field components from geocentric C to geomagnetic directions. C =(sin(colat)*sin(clon-wlon))/cos(gmlat) C where colat and wlon are the geocentric colat and lon. C of the north pole of the earth-centered dipole. C INPUT IF I .EQ. 1 C GMLAT Geomagnetic latitude (Deg) C GMLON Geomagnetic longitude (Deg) C OUTPUT IF I .EQ. 1 C CLAT Geocentric latitude (Deg) C CLON Geocentric longitude (Deg) C SIND Constant which may be used in conversion of C X and Y field components from geocentric C to geomagnetic directions. C =(cos(colat)-sin(gmlat)sin(clat))/cos(gmlat)cos(clat) C where colat is the geocentric colat of the north pole C of the earth-centered dipole. C COSD Constant which may be used in conversion of C X and Y field components from geocentric C to geomagnetic directions. C =(sin(colat)*sin(clon-wlon))/cos(gmlat) C where colat and wlon are the geocentric colat and lon. C of the north pole of the earth-centered dipole. C***LONG DESCRIPTION C C C COMMON Block Used C /MAG / ICODE,XI(3),H(144),DUM1(2,40),NAME(4),NMAX,TIME,G(144),DUM2(4) C ICODE Flag set in STOER subprogram. C XI Work array used in subprograms FELDG, SHELG and STOER C H Work array used in subprogram FELDG , SHELG and STOER C DUM1 Work array used in SHELG and STOER. C NAME Not being used any more , used in old version of COFRM C NMAX Order of IGRF model set in COFRM subprogram C TIME Equal to input argument DATE in COFRM routine. C G Array of spherical harmonic coefficients for given C date derived in COFRM subprogram. C DUM2 Array of numbers set in COFRM routine and used in SHELG C C C /CONST/ RTOD,DTOR,RE,REQ,COLAT,WLON C RTOD Unit to convert radian to degrees, 45./ATAN(1.). C DTOR Unit to convert degrees to radians ATAN(1.)/45. C RE Radius of the earth, 6371.2 (Km). C REQ Equatorial radius, 6378.165. C COLAT Geographic colatitude of the north pole of the C earth-centered dipole (Deg). C WLON Geographic longitude of the north pole of the C earth-centered dipole (Deg). C C***REFERENCES reference 1 C C***ROUTINES CALLED None C***COMMON BLOCKS CONST,MAG C***END PROLOGUE TSFORM COMMON / CONST / RTOD,DTOR,RE,REQ,COLAT,WLON COMMON/MAG/IDUM1,DUM1(307),NAME(4),NMAX,TIME,G(12,12),DUM2(4) C COMMON/MAGCOF/ NMAX,TIME,G(144) C STO = SIN(COLAT*DTOR) CTO = COS(COLAT*DTOR) GO TO (10, 20) I C C ** GEOCENTRIC TO GEOMAGNETIC C ** 10 CTG = SIN(CLAT*DTOR) STG = COS(CLAT*DTOR) CTD = CTO*CTG + STO*STG*COS((CLON-WLON)*DTOR) STD = SQRT(1-(CTD*CTD)) GMLAT = RTOD*ATAN2(CTD,STD) C**************** C CLD = (CTO*CTD-CTG)/(STO*STD) C SLGLO = SIN((CLON-WLON)*DTOR) C SLD = STG*SLGLO/STD C GMLON = 180. - RTOD*ATAN2(SLD, -1*CLD) C IF (GMLON .GT. 180.) GMLON = GMLON - 360. CLD = (CTO*CTD-CTG) SLGLO = SIN((CLON-WLON)*DTOR) SLD = STG*SLGLO*STO GMLON = RTOD*ATAN2(SLD,CLD) C**************** GO TO 30 C C ** GEOMAGETIC TO GEOCENTRIC C ** 20 CTD = SIN(GMLAT*DTOR) STD = COS(GMLAT*DTOR) CTG = CTO*CTD - STO*STD*COS(GMLON*DTOR) STG = SQRT(1-(CTG*CTG)) CLAT = RTOD*ATAN2(CTG,STG) CLGLO = (CTD-CTO*CTG) SLGLO = STD*SIN(GMLON*DTOR)*STO SWLON = SIN(SWLON*DTOR) CWLON = COS(SWLON*DTOR) CLON = RTOD*ATAN2(SLGLO*CWLON+CLGLO*SWLON,CLGLO*CWLON-SLGLO* 1 SWLON) C IF (CLON .LT. 0.) CLON = CLON + 360. C 30 COSD = (CTO-CTD*CTG) SIND = STO*SLGLO*STG DMAG = RTOD*ATAN2(-SIND,COSD) RETURN C END C SUBROUTINE GTM(RLATG,RLONG,RLATMP,RLONMP,RLATM,RLONM,DIP,DEC,W,N) C **** TRANSFORMS GEOGRAPHIC COORDINATES TO GEOMAGNETIC, C **** CALCULATES DIP AND DECLINATION. C **** C **** RLATG(N), RLONG(N) -- N-DIMENSIONAL ARRAYS OF GEOGRAPHIC C **** LATITUDES AND CORRESPONDING LONGITUDES. (RADIANS) C **** C **** RLATMP, RLONMP -- GEOGRAPHIC COORDINATES OF GEOMAGNETIC C **** NORTH POLE. (RADIANS) C **** C **** RLATM(N), RLONM(N) -- OUTPUT ARRAYS OF GEOMAGNETIC C **** COORDINATES. (RADIANS) C **** C **** DIP(N), DEC(N) -- OUTPUT ARRAYS OF DIP AND DECLINATION. C **** C **** W(N,4) -- WORKING ARRAY. C **** C **** N -- NUMBER OF (LAT,LON) POINTS TO BE TRANSFORMED. C **** DIMENSION RLATG(N),RLONG(N),RLATM(N),RLONM(N),DIP(N),DEC(N),W(N,4) DATA E/1.E-10/ SINLAM=SIN(RLATMP) COSLAM=COS(RLATMP) DO 1 I=1,N W(I,1)=SIN(RLATG(I)) W(I,2)=COS(RLATG(I)) W(I,3)=SIN(RLONG(I)+E-RLONMP) W(I,4)=COS(RLONG(I)+E-RLONMP) W(I,4)=SINLAM*W(I,1)+COSLAM*W(I,2)*W(I,4) RLONM(I)=ATAN2(W(I,3)*W(I,2),(W(I,4)*SINLAM-W(I,1))/COSLAM) RLATM(I)=ASIN(W(I,4)) C DEC(I)=ATAN2(COSLAM,(SINLAM-W(I,1)*W(I,4))/(W(I,2)*W(I,3))) DEC(I)=ATAN2(COSLAM*W(I,2)*W(I,3),SINLAM-W(I,1)*W(I,4)) DIP(I)=ATAN(2.*TAN(RLATM(I))) 1 CONTINUE RETURN ENTRY MTG(RLATG,RLONG,RLATMP,RLONMP,RLATM,RLONM,W,N) C **** INVERSE TRANSFORMATION C **** C **** PARAMETERS AS ABOVE EXCEPT INPUT -- RLATM,RLONM C **** OUTPUT -- RLATG,RLONG C **** SINLAM=SIN(RLATMP) COSLAM=COS(RLATMP) SINLOM=SIN(RLONMP) COSLOM=COS(RLONMP) DO 2 I=1,N W(I,1)=SIN(RLATM(I)) W(I,2)=COS(RLATM(I)) W(I,3)=SIN(RLONM(I)) W(I,4)=COS(RLONM(I)) RLATG(I)=SINLAM*W(I,1)-COSLAM*W(I,2)*W(I,4) W(I,3)=W(I,3)*W(I,2) W(I,4)=-(RLATG(I)*SINLAM-W(I,1))/COSLAM RLATG(I)=ASIN(RLATG(I)) RLONG(I)=ATAN2(SINLOM*W(I,4)+COSLOM*W(I,3),COSLOM*W(I,4)-SINLOM A*W(I,3)) 2 CONTINUE RETURN END