subroutine wmm2005(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), c 2004-03-25 (external cosd, sind), 2006-03-27 (wmm2005); adapted c from original routine shval3, from the igrf program suite geomag 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=2453371.5d0) !epoch 2005.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 / &-29556.8d0, -1671.7d0, 5079.8d0, -2340.6d0, 3046.9d0, -2594.7d0, & 1657.0d0, -516.7d0, 1335.4d0, -2305.1d0, -199.9d0, 1246.7d0, & 269.3d0, 674.0d0, -524.2d0, 919.8d0, 798.1d0, 281.5d0, & 211.3d0, -226.0d0, -379.4d0, 145.8d0, 100.0d0, -304.7d0, & -227.4d0, 354.6d0, 42.4d0, 208.7d0, 179.8d0, -136.5d0, & -123.0d0, -168.3d0, -19.5d0, -14.1d0, 103.6d0, 73.2d0, & 69.7d0, -20.3d0, 76.7d0, 54.7d0, -151.2d0, 63.6d0, & -14.9d0, -63.4d0, 14.6d0, -0.1d0, -86.3d0, 50.4d0, & 80.1d0, -74.5d0, -61.5d0, -1.4d0, -22.4d0, 38.5d0, & 7.2d0, 12.4d0, 25.4d0, 9.5d0, 11.0d0, 5.7d0, & -26.4d0, 1.8d0, -5.1d0, 24.9d0, 7.7d0, 11.2d0, & -11.6d0, -21.0d0, -6.9d0, 9.6d0, -18.2d0, -19.8d0, & 10.0d0, 16.1d0, 9.2d0, 7.7d0, -11.6d0, -12.9d0, & -5.2d0, -0.2d0, 5.6d0, 9.9d0, -20.1d0, 3.5d0, & 12.9d0, -7.0d0, 12.6d0, 5.1d0, -6.7d0, -10.8d0, & -8.1d0, -1.3d0, 8.0d0, 8.8d0, 2.9d0, -6.7d0, & -7.9d0, -9.1d0, 6.0d0, -2.3d0, -6.3d0, 2.4d0, & 1.6d0, 0.2d0, -2.6d0, 4.4d0, 0.0d0, 4.8d0, & 3.1d0, -6.5d0, 0.4d0, -1.1d0, 2.1d0, -3.4d0, & 3.9d0, -0.8d0, -0.1d0, -2.3d0, -2.3d0, -7.9d0, & 2.8d0, -1.6d0, 0.3d0, -1.7d0, 1.2d0, 1.7d0, & -0.8d0, -0.1d0, -2.5d0, 0.1d0, 0.9d0, -0.7d0, & -0.6d0, 0.7d0, -2.7d0, 1.8d0, -0.9d0, 0.0d0, & -1.3d0, 1.1d0, -2.0d0, 4.1d0, -1.2d0, -2.4d0, & -0.4d0, -0.4d0, 0.2d0, 0.3d0, 0.8d0, 2.4d0, & -0.3d0, -2.6d0, 1.1d0, 0.6d0, -0.5d0, 0.3d0, & 0.4d0, 0.0d0, -0.3d0, 0.0d0, -0.3d0, 0.3d0, & -0.1d0, -0.9d0, -0.3d0, -0.4d0, -0.1d0, 0.8d0/ ! data dgh / & 8.0d0, 10.6d0, -20.9d0, -15.1d0, -7.8d0, -23.2d0, & -0.8d0, -14.6d0, 0.4d0, -2.6d0, 5.0d0, -1.2d0, & -7.0d0, -6.5d0, -0.6d0, -2.5d0, 2.8d0, 2.2d0, & -7.0d0, 1.6d0, 6.2d0, 5.8d0, -3.8d0, 0.1d0, & -2.8d0, 0.7d0, 0.0d0, -3.2d0, 1.7d0, -1.1d0, & 2.1d0, 0.1d0, 4.8d0, -0.8d0, -1.1d0, -0.7d0, & 0.4d0, -0.6d0, -0.3d0, -1.9d0, 2.3d0, -0.4d0, & -2.1d0, -0.5d0, -0.6d0, -0.3d0, 1.4d0, 0.7d0, & 0.2d0, -0.1d0, 0.6d0, -0.3d0, 0.4d0, 1.1d0, & 0.2d0, 0.6d0, 0.3d0, 0.5d0, -0.8d0, -0.4d0, & -0.2d0, 0.6d0, 0.1d0, 0.1d0, 0.3d0, -0.2d0, & -0.4d0, 0.1d0, 0.3d0, 0.3d0, -0.3d0, 0.4d0, & 0.2d0, 0.1d0, 0.4d0, -0.2d0, -0.7d0, 0.4d0, & 0.4d0, 0.4d0, 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