subroutine wmm2000(jd, lat, lon, r, xyz) c c ================================================================ c c version 1.01 c c calculates field components from spherical harmonic (sh) models c c c input: c jd - julian date (Julian day number and fraction) c lat - north geocentric latitude, in degrees c lon - east geocentric longitude, in degrees c r - radial distance from earth's center (km) c c output: c xyz - [northward, eastward, downward] components c c based on subroutine 'igrf' by d. r. barraclough and c s. r. c. malin, report no. 71/1, institute of geological c sciences, u.k. c c norman w. peddie, u.s. geological survey, mail stop 964, c federal center, box 25046, denver, colorado 80225 c c b. g. knapp, 1998-11-16 (wmm1995), 2001-04-18 (wmm2000), adapted c from original routine shval3, from the igrf program suite geomag; c 2004-03-25 (external cosd, sind) c c ================================================================ c c the required sizes of the arrays used in this subroutine c depend on the value of nmax. the minimum dimensions c needed are indicated in the table below. c c nmax - maximum degree and order of coefficients c gh - schmidt quasi-normal internal spherical c harmonic coefficients c c n m g h c ---------------------- c / 1 0 gh(1) - c / 1 1 gh(2) gh(3) c npq = nmax*(nmax+3)/2 / 2 0 gh(4) - c records / 2 1 gh(5) gh(6) c / 2 2 gh(7) gh(8) c \ 3 0 gh(9) - c ngh = nmax*(nmax+2) \ . . . . c elements in gh \ . . . . c \ . . . . c \ nmax nmax . . c c n and m are, respectively, the degree and order of the c coefficient. ! C C RCS DATA C C $Header$ C C $Log$ C C implicit none ! ! Input real*8 jd, lat, lon, r ! ! Output real*8 xyz(3) ! ! Constants integer*4 nmax, ngh, npq parameter (nmax=12, ngh=nmax*(nmax+2), npq=(nmax*(nmax+3))/2) real*8 erad parameter (erad=6371.2d0) !mean radius of Earth (km) real*8 jd0 parameter (jd0=2451544.5d0) !epoch 2000.0 ! ! Local real*8 gh0(ngh), dgh(ngh), gh(ngh) real*8 p(npq), q(npq) real*8 sl(nmax), cl(nmax) real*8 rr, ratio, slat, clat, sd, cd, aa, bb, cc, f, fn, fm real*8 x, y, z integer*4 i,j,k,l,m,n save gh0, dgh ! data gh0 / &-29616.0d0, -1722.7d0, 5194.5d0, -2266.7d0, 3070.2d0, -2484.8d0, & 1677.6d0, -467.9d0, 1322.4d0, -2291.5d0, -224.7d0, 1255.9d0, & 293.0d0, 724.8d0, -486.5d0, 932.1d0, 786.3d0, 273.3d0, & 250.6d0, -227.9d0, -401.5d0, 120.9d0, 106.2d0, -302.7d0, & -211.9d0, 351.6d0, 42.0d0, 220.8d0, 173.8d0, -134.5d0, & -135.0d0, -168.8d0, -38.6d0, -13.3d0, 105.2d0, 73.8d0, & 68.2d0, -17.4d0, 74.1d0, 61.2d0, -163.5d0, 63.2d0, & -3.8d0, -62.9d0, 17.1d0, 0.2d0, -85.1d0, 43.0d0, & 77.4d0, -73.9d0, -62.3d0, 2.2d0, -24.5d0, 35.7d0, & 8.9d0, 7.3d0, 23.4d0, 5.2d0, 15.0d0, 8.4d0, & -27.6d0, -1.5d0, -7.8d0, 23.3d0, 7.3d0, 12.4d0, & -8.5d0, -20.8d0, -6.6d0, 8.4d0, -16.9d0, -21.2d0, & 8.6d0, 15.5d0, 4.9d0, 9.1d0, -7.8d0, -15.5d0, & -7.6d0, -5.4d0, 5.7d0, 8.5d0, -20.4d0, 2.0d0, & 13.9d0, -9.8d0, 12.0d0, 7.6d0, -6.2d0, -7.0d0, & -8.6d0, -2.0d0, 9.4d0, 9.2d0, 5.0d0, -2.2d0, & -8.4d0, -6.6d0, 3.2d0, -2.2d0, -5.7d0, 0.9d0, & 1.6d0, -0.7d0, -3.7d0, 3.9d0, -0.6d0, 4.8d0, & 4.1d0, -5.3d0, 2.2d0, -1.0d0, 2.2d0, -2.4d0, & 4.6d0, 1.3d0, 2.3d0, -2.3d0, 0.1d0, -6.4d0, & 3.3d0, -1.1d0, -1.5d0, -2.4d0, 0.7d0, 2.6d0, & -1.1d0, -1.3d0, -2.3d0, -1.7d0, 1.3d0, -0.6d0, & -0.6d0, 0.4d0, -2.8d0, 0.7d0, -1.6d0, -0.3d0, & -0.1d0, 2.3d0, -1.9d0, 4.2d0, 1.4d0, -1.5d0, & -0.2d0, -1.0d0, -0.3d0, 0.7d0, 0.5d0, 2.2d0, & 0.2d0, -2.5d0, 0.9d0, -0.2d0, -1.4d0, 0.0d0, & 0.6d0, -0.2d0, -0.6d0, 0.0d0, -1.0d0, 0.2d0, & -0.3d0, -0.9d0, 0.3d0, -0.2d0, 0.4d0, 1.0d0/ ! data dgh / & 14.7d0, 11.1d0, -20.4d0, -13.6d0, -0.7d0, -21.5d0, & -1.8d0, -9.6d0, 0.3d0, -4.3d0, 6.4d0, 0.9d0, & -1.3d0, -8.4d0, -13.3d0, -1.6d0, 0.9d0, 2.3d0, & -7.6d0, 0.7d0, 2.2d0, 3.7d0, -3.2d0, -0.5d0, & -0.9d0, -0.2d0, 0.0d0, -2.5d0, 2.1d0, -2.7d0, & 2.3d0, -0.9d0, 3.1d0, 1.7d0, 0.0d0, 1.2d0, & 0.2d0, -0.3d0, 1.7d0, -1.7d0, 1.6d0, -0.9d0, & -0.1d0, -1.0d0, -0.3d0, -0.1d0, 0.8d0, 1.9d0, & -0.4d0, -0.8d0, 1.4d0, -0.2d0, 0.2d0, 1.1d0, & 0.7d0, 0.4d0, 0.4d0, 0.0d0, -0.3d0, -0.2d0, & -0.8d0, -0.2d0, -0.1d0, -0.3d0, 0.6d0, -0.5d0, & -0.8d0, 0.1d0, 0.3d0, -0.2d0, -0.2d0, 0.0d0, & 0.5d0, 0.1d0, 0.0d0, -0.1d0, -0.6d0, 0.3d0, & 0.1d0, 0.2d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, & 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0/ c c External functions real*8 cosd, sind external cosd, sind c c ================================================================ c c Apply secular variation (dgh) for requested date f = (jd-jd0)/365.25d0 do l=1,ngh gh(l) = gh0(l)+f*dgh(l) enddo c c Initialize slat = sind(lat) clat = cosd(lat) sl(1) = sind(lon) cl(1) = cosd(lon) x = 0.d0 y = 0.d0 z = 0.d0 sd = 0.d0 cd = 1.d0 n = 0 l = 1 m = 1 c ratio = erad/r aa = sqrt(3.d0) p(1) = 2.d0*slat p(2) = 2.d0*clat p(3) = 4.5d0*slat**2 - 1.5d0 p(4) = 3.d0*aa*clat*slat q(1) = -clat q(2) = slat q(3) = -3.d0*clat*slat q(4) = aa*(slat+clat)*(slat-clat) c c Evaluate the spherical harmonic function at the given arguments do k = 1,npq if (n .lt. m) then m = 0 n = n + 1 rr = ratio**(n+2) fn = n endif fm = m if (k .ge. 5) then if (m .eq. n) then aa = sqrt(1.d0-.5d0/fm) j = k - n - 1 p(k) = (1.d0+1.d0/fm)*aa*clat*p(j) q(k) = aa*(clat*q(j) + slat/fm*p(j)) sl(m) = sl(m-1)*cl(1) + cl(m-1)*sl(1) cl(m) = cl(m-1)*cl(1) - sl(m-1)*sl(1) else aa = sqrt((fn+fm)*(fn-fm)) bb = sqrt((fn-1.d0+fm)*(fn-1.d0-fm))/aa cc = (2.d0*fn-1.d0)/aa i = k - n j = k - 2*n + 1 p(k) = (fn+1.d0)*(cc*slat/fn*p(i) - bb/(fn-1.d0)*p(j)) q(k) = cc*(slat*q(i) - clat/fn*p(i)) - bb*q(j) endif endif aa = rr*gh(l) if (m .eq. 0) then x = x + aa*q(k) z = z - aa*p(k) l = l + 1 else bb = rr*gh(l+1) cc = aa*cl(m) + bb*sl(m) x = x + cc*q(k) z = z - cc*p(k) if (clat .gt. 0.d0) then y = y + (aa*sl(m) - bb*cl(m))*fm*p(k)/((fn+1.d0)*clat) else y = y + (aa*sl(m) - bb*cl(m))*q(k)*slat endif l = l + 2 endif m = m + 1 enddo xyz(1) = x*cd + z*sd xyz(2) = y xyz(3) = z*cd - x*sd return end