73,76c73,79 < C Assign DGRF/IGRF spherical harmonic coefficients, to degree and < C order NMAX, for DATE, yyyy.fraction, into array G. Coefficients < C are interpolated from the DGRF dates through the current IGRF year. < C Coefficients for a later DATE are extrapolated using the IGRF --- > C***DATE WRITTEN 830415 (YYMMDD) > C***LAST REVISION DATE 920430 (YYMMDD) BY R. BARNES, 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 78,87c81,84 < C message is issued if DATE is later than the last recommended < C (5 yrs later than the IGRF). An DATE input earlier than the < C first DGRF (EPOCH(1)), results in a diagnostic and a STOP. < C < C Output in COMMON /MAGCOF/ NMAX,GB(144),GV(144),ICHG < C NMAX = Maximum order of spherical harmonic coefficients used < C GB = Coefficients for magnetic field calculation < C GV = Coefficients for magnetic potential calculation < C ICHG = Flag indicating when GB,GV have been changed in COFRM < C --- > 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. > 89,96c86,88 < C COFRM and FELDG originated 15 Apr 83 by Vincent B. Wickwar < C (formerly at SRI. Int., currently at Utah State). Although set < C up to accomodate second order time derivitives, the IGRF < C (GTT, HTT) have been zero. The spherical harmonic coefficients < C degree and order is defined by NMAX (currently 10). < C < C Jun 86: Updated coefficients adding DGRF 1980 & IGRF 1985, which < C were obtained from Eos Vol. 7, No. 24. Common block MAG was --- > 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 98,102c90,95 < C FELDG. (Roy Barnes) < C < C Apr 1992 (Barnes): Added DGRF 1985 and IGRF 1990 as described < C in EOS Vol 73 Number 16 Apr 21 1992. Other changes were made so < C future updates should: --- > 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: 108,121c101,103 < C Apr 94 (Art Richmond): Computation of GV added, for finding < C magnetic potential. < C < C Aug 95 (Barnes): Added DGRF for 1990 and IGRF for 1995, which were < C obtained by anonymous ftp geomag.gsfc.nasa.gov (cd pub, mget table*) < C as per instructions from Bob Langel (langel@geomag.gsfc.nasa.gov), < C but, problems are to be reported to baldwin@geomag.gsfc.nasa.gov < < C Oct 95 (Barnes): Correct error in IGRF-95 G 7 6 and H 8 7 (see < C email in folder). Also found bug whereby coefficients were not being < C updated in FELDG when IENTY did not change. ICHG was added to flag < C date changes. Also, a vestigial switch (IS) was removed from COFRM: < C It was always 0 and involved 3 branch if statements in the main < C polynomial construction loop (now numbered 200). --- > 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). 123c105,113 < C REAL :: RTOD=5.72957795130823E1,DTOR=1.745329251994330E-2 --- > 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******************************************************************** 125,127d114 < COMMON /MAGCOF/ NMAX,G(144),GV(144),ICHG < COMMON /CONST/ RTOD,DTOR,RE,REQ,COLAT,WLON < DATA NMAX,ICHG /10,-99999/ 129c116,123 < PARAMETER (NDGY=6 , NYT = NDGY+1 , NGH = 144*NDGY) --- > 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) 140,141c134 < SAVE DATEL,GYR,HYR,EPOCH,G1DIM,H1DIM,G0,H0,GT,HT,GTT,HTT < C SAVE DATEL,GYR,HYR,EPOCH, G0,H0,GT,HT,GTT,HTT --- > SAVE DATEL,GYR,HYR,G0,H0,GT,HT,GTT,HTT,EPOCH 143c136 < DATA EPOCH /1965. , 1970. , 1975. , 1980. , 1985. , 1990. , 1995./ --- > DATA EPOCH /1965. , 1970. , 1975. , 1980. , 1985. , 1990./ 150d142 < 215,227d206 < C DGRF g(n,m) for 1990: < DATA (G1DIM(I),I=721,864)/ 0., < O -29775.,-2131., 1314., 939., -214., 61., 77., 23., 4., -3., 2*0., < 1 -1848., 3059.,-2239., 780., 353., 65.,-64., 5., 9., -4., 3*0., < 2 1686., 1248., 325., 245., 59., 2., -1., 1., 2., 4*0., < 3 802.,-423., -109.,-178.,26.,-10.,-12.,-5., 5*0., < 4 141., -165., 3., -1.,-12., 9., -2., 6*0., < 5 -36., 18., 5., 3.,-4., 4., 7*0., < 6 -96., 9., 4.,-2., 3., 8*0., < 7 0., 2., 7., 1., 9*0., < 8 -6., 1., 3.,10*0., < 9 -6., 3.,11*0., < A 0.,13*0./ 277d255 < C DGRF h(n,m) for 1985: 289,303c267 < C DGRF h(n,m) for 1990: < DATA (H1DIM(I),I=721,864)/ 13*0., < 1 5406., -2279., -284., 247., 46.,-16.,-80., 10.,-20., 2., 3*0., < 2 -373., 293.,-240., 154., 82.,-26.,-19., 15., 1., 4*0., < 3 -352., 84.,-153., 69., 0., 6., 11., 3., 5*0., < 4 -299., -69.,-52., 21.,-22., -7., 6., 6*0., < 5 97., 1., 17., 12., -7.,-4., 7*0., < 6 24.,-23., 12., 9., 0., 8*0., < 7 -4.,-16., 8.,-2., 9*0., < 8 -10., -7., 3.,10*0., < 9 2.,-1.,11*0., < A -6.,13*0./ < < < C Initial coefficients g0 (IGRF for 1995): --- > C Initial coefficients g0 (IGRF for 1990): 305,313c269,277 < O -29682., -2197., 1329., 941.,-210., 66., 78., 24., 4.,-3., 2*0., < 1 -1789., 3074.,-2268., 782., 352., 64.,-67., 4., 9.,-4., 3*0., < 2 1685., 1249., 291., 237., 65., 1., -1., 1., 2., 4*0., < 3 769.,-421.,-122.,-172.,29., -9.,-12.,-5., 5*0., < 4 116.,-167., 2., 4.,-14., 9.,-2., 6*0., < 5 -26., 17., 8., 4., -4., 4., 7*0., < 6 -94., 10., 5., -2., 3., 8*0., < 7 -2., 0., 7., 1., 9*0., < 8 -7., 0., 3.,10*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., 316c280 < C D_/Dtime coefficients gt (IGRF for 1995-2000): --- > C D_/Dtime coefficients gt (IGRF for 1990-1995): 318,326c282,290 < O 17.6, -13.2, 1.5, .8, .8, .5, -.2, .3, 0.0,0.0, 2*0., < 1 13.0, 3.7, -6.4, .9, .1, -.4, -.8, -.2, 0.0,0.0, 3*0., < 2 -.8, -.2, -6.9, -1.5, .6, -.6, .1, 0.0,0.0, 4*0., < 3 -8.1, .5, -2.0, 1.9, .6, .4, 0.0,0.0, 5*0., < 4 -4.6, -.1, -.2, 1.2,-1.1, 0.0,0.0, 6*0., < 5 2.3, -.2, .1, .3, 0.0,0.0, 7*0., < 6 .0, .2, .2, 0.0,0.0, 8*0., < 7 -.6, -.9, 0.0,0.0, 9*0., < 8 -.3, 0.0,0.0,10*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., 329c293 < C Initial coefficients h0 (IGRF for 1995-2000): --- > C Initial coefficients h0 (IGRF for 1990-1995): 331,339c295,303 < 1 5318., -2356., -263., 262., 44.,-16.,-77., 12.,-19., 2., 3*0., < 2 -425., 302.,-232.,157., 77.,-25.,-20., 15., 1., 4*0., < 3 -406., 98.,-152., 67., 3., 7., 11., 3., 5*0., < 4 -301.,-64.,-57., 22.,-21., -7., 6., 6*0., < 5 99., 4., 16., 12., -7.,-4., 7*0., < 6 28.,-23., 10., 9., 0., 8*0., < 7 -3.,-17., 7.,-2., 9*0., < 8 -10., -8., 3.,10*0., < 9 1.,-1.,11*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., 341c305 < C D_/Dtime coefficients ht (IGRF for 1995-2000): --- > C D_/Dtime coefficients ht (IGRF for 1990-1995): 343,350c307,314 < 1 -18.3, -15.0, 4.1, 1.8, .2, .3, .8, .4, 0.0,0.0, 3*0., < 2 -8.8, 2.2, 1.2,1.2,-1.6, .2, -.2, 0.0,0.0, 4*0., < 3 -12.1, 2.7, .3, -.2, .6, .2, 0.0,0.0, 5*0., < 4 -1.0,1.8, -.9, -.4, .7, 0.0,0.0, 6*0., < 5 .9, 1.0, .0, .0, 0.0,0.0, 7*0., < 6 2.2, -.3,-1.2, 0.0,0.0, 8*0., < 7 .0, -.7, 0.0,0.0, 9*0., < 8 -.6, 0.0,0.0,10*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., 353,359c317,323 < < RTOD=5.72957795130823E1 < DTOR=1.745329251994330E-2 < < C Do not need to load new coefficients if date has not changed < ICHG = 0 < IF (DATE .EQ. DATEL) GO TO 300 --- > 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 361,362c325 < ICHG = 1 < --- > 365,367c328 < c IF (DATE .GT. EPOCH(NYT)+5.) WRITE(6,9200) DATE < IF (DATE .GT. EPOCH(NYT)+5.) < + WRITE(6,"('COFRM: DATE',F9.3,' is after the maximum')") DATE --- > IF (DATE .GT. EPOCH(NYT)+5.) WRITE(6,9200) DATE 374a336 > NMAX = 10 375a338 > IS = 0 379,384c342,353 < F0 = -1.0D-5 < DO 200 N=1,NMAX < RN = REAL (N) < F0 = F0*RN*(2.*RN-1.) / (4.*RN-2.) < F = 0.5*SQRT(2.0)*F0 < NN = N+1 --- > 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 388c357 < G(I) = ((GTT(NN,MM)*T + GT(NN,MM))*T + G0(NN,MM)) * F0 --- > G(I) = ((GTT(NN,MM)*T + GT(NN,MM))*T + G0(NN,MM)) * F0 391c360 < G(I) = (GYR(NN,MM,IY) + --- > G(I) = (GYR(NN,MM,IY) + 394d362 < GV(I) = G(I)/FLOAT(NN) 396,400c364,370 < DO 200 M=1,N < TMP1 = REAL (N+M) < TMP2 = REAL (N-M+1) < F = F*SQRT(TMP2/TMP1)*TMP1/TMP2 < NN = N+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 405,406c375,376 < 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 --- > 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 409c379 < G(I) = (GYR(NN,MM,IY) + --- > G(I) = (GYR(NN,MM,IY) + 411c381 < G(I1) = (HYR(NN,MM,IY) + --- > G(I1) = (HYR(NN,MM,IY) + 414,418c384,386 < GV(I) = G(I) / REAL(NN) < GV(I1) = G(I1) / REAL(NN) < 200 I = I+2 < < 300 CONTINUE --- > 600 I = I+2 > 990 CONTINUE > C COMPUTE GEOGRAPHIC COLATITUDE AND LONGITUDE OF THE NORTH POLE OF 426d393 < 428c395 < --- > 438d404 < C