bfield2.f
676 lines
| 20.8 KiB
| text/x-fortran
|
FortranFixedLexer
| r1601 | subroutine geobfield(tm,r,theta,phi,br,bt,bp,b) | |||
| c*****evaluate bfield at geocentric coords r,theta,phi | ||||
| c and return components of field in r,theta and phi directions | ||||
| c and magnitude | ||||
| c*****theta and phi are in radians, r is in km, b is in gauss | ||||
| c =============================================================== | ||||
| c | ||||
| c essentially, this is an adaptation of igrfdemo. it interpolates | ||||
| c smoothly between the coefficients from 1945 to 1985. | ||||
| c | ||||
| c subroutines: | ||||
| c getshc, interpshc, extrapshc, shval3 | ||||
| c | ||||
| c igrfdemo is due to: | ||||
| c a. zunde | ||||
| c usgs, ms 964, box 25046 federal center, denver, co 80225 | ||||
| c | ||||
| c =============================================================== | ||||
| character*8 filmod(17) | ||||
| c character*31 fqual | ||||
| dimension gh1(220), gh2(220), gha(224), ext(3), dtemod(17) | ||||
| data ext /3*0./ | ||||
| data filmod / 'dgrf45', 'dgrf50', | ||||
| 1 'dgrf55', 'dgrf60', 'dgrf65', | ||||
| 2 'dgrf70', 'dgrf75', 'dgrf80', | ||||
| 3 'dgrf85', 'dgrf90', 'dgrf95', | ||||
| 4 'dgrf00', 'dgrf05', 'dgrf10', | ||||
| 5 'dgrf15', 'igrf20', 'igrf20s'/ | ||||
| c data fqual/"/usr/local/lib/faraday/bfmodel/"/ | ||||
| data dtemod / 1945., 1950., 1955., 1960., | ||||
| 1 1965., 1970., 1975., 1980., 1985., 1990., 1995., 2000., | ||||
| 2 2005.,2010.,2015.,2020.,2025./ | ||||
| data iu/98/,ndt/17/ | ||||
| c data a2/40680925./, b2/40408588./ | ||||
| data a2/40680631.6/, b2/40408296.0/ ! updated 2010 | ||||
| c | ||||
| data rtd/57.29577951/ | ||||
| data tmp/0./,lp/0/ | ||||
| character(1024) :: fqual_temp | ||||
| character(:), allocatable :: fqual | ||||
| call get_path(fqual_temp) | ||||
| c write(*,*) "L_BEF: ", fqual_temp, "L_BEF_end" | ||||
| fqual = TRIM(fqual_temp) | ||||
| c write(*,*) "L: ", fqual, "L_end" | ||||
| flat=90. - theta*rtd | ||||
| flon=rtd*phi | ||||
| c*****if previous time is not equal to current time | ||||
| if(tm .ne. tmp)then | ||||
| l=0 | ||||
| C for i=1,ndt | ||||
| DO 1010 i=1,ndt | ||||
| C exit for if(tm .lt. dtemod(i)) | ||||
| if(tm .lt. dtemod(i)) GO TO 1011 | ||||
| l=i | ||||
| C end for | ||||
| 1010 CONTINUE | ||||
| 1011 CONTINUE | ||||
| c write(*,-)tm,l | ||||
| if(l .eq. 0)then | ||||
| write(*,fmt='(" geobfield: time is before earliest model.")') | ||||
| stop | ||||
| end if | ||||
| if(l .eq. ndt)then | ||||
| l=l-1 | ||||
| write(*,fmt='(" geobfield: warning - extrapolating beyond ", | ||||
| 1 "last set of coefficients.")') | ||||
| end if | ||||
| if(l .ne. lp)then | ||||
| c*********if previous epoch not the same, read in new coefs | ||||
| c write(*,fmt='(" read coefs"))') | ||||
| c write(*,*) "filmod", fqual//filmod(l) | ||||
| call getshc2 (iu, fqual//filmod(l), nmax1, erad, gh1, ier) | ||||
| c write(*,*) "AA: ", ier, filmod(l) | ||||
| if(ier .ne. 0)then | ||||
| write(*,fmt='(" geobfield: read error=",i2," on ",a)') | ||||
| 1 ier,filmod(l) | ||||
| stop | ||||
| end if | ||||
| call getshc2 (iu, fqual//filmod(l+1), nmax2, erad, gh2, ier) | ||||
| if(ier .ne. 0)then | ||||
| write(*,fmt='(" geobfield: read error=",i2," on ",a)') | ||||
| 1 ier,filmod(l+1) | ||||
| stop | ||||
| end if | ||||
| end if | ||||
| if (l .lt. ndt-1) then | ||||
| call interpshc (tm, dtemod(l), nmax1, gh1, dtemod(l+1), | ||||
| 1 nmax2, gh2, nmax, gha) | ||||
| c write(*,fmt='(" interpolate"))') | ||||
| else | ||||
| call extrapshc (tm, dtemod(l), nmax1, gh1, nmax2, | ||||
| 1 gh2, nmax, gha) | ||||
| c write(*,fmt='(" extrapolate"))') | ||||
| end if | ||||
| end if | ||||
| c tmp=tm | ||||
| c lp=l | ||||
| call shval3(2,flat,flon,r,erad,a2,b2,nmax,gha,0,ext,x,y,z) | ||||
| br=-z/100000. | ||||
| bp=y/100000. | ||||
| bt=-x/100000. | ||||
| b=sqrt(br**2 + bp**2 + bt**2) | ||||
| C write(*,*)flat,flon,r,erad,a2,b2,nmax,gha,ext | ||||
| C write(*,*)tm,r,theta,phi,br,bt,bp,b | ||||
| return | ||||
| end | ||||
| c | ||||
| subroutine convrt(i, gdlat, gdalt, gclat, rkm) | ||||
| c convrt converts between geodetic and geocentric coordinates. the | ||||
| c reference geoid is that adopted by the iau in 1964. a=6378.16, | ||||
| c b=6356.7746, f=1/298.25. the equations for conversion from | ||||
| c geocentric to geodetic are from astron. j., vol 66, 1961, p. 15. | ||||
| c i=1 geodetic to geocentric | ||||
| c i=2 geocentric to geodetic | ||||
| c gdlat geodetic latitude (degrees) | ||||
| c gdalt altitude above geoid (km) | ||||
| c gclat geocentric latitude (degrees) | ||||
| c rkm geocentric radial distance (km) | ||||
| c | ||||
| c data a/6378.16/, ab2/1.0067397/, ep2/.0067397/ | ||||
| c update 2010 | ||||
| data a/6378.137/, ab2/1.0067396/, ep2/.0067396/ | ||||
| data dtr/.0174532925199/ | ||||
| if (i .eq. 1) then | ||||
| c .....geodetic to geocentric..... | ||||
| gdl = dtr*gdlat | ||||
| sinlat = sin(gdl) | ||||
| coslat = cos(gdl) | ||||
| sl2 = sinlat*sinlat | ||||
| cl2 = ab2*coslat | ||||
| cl2 = cl2*cl2 | ||||
| sinbet = sinlat/sqrt(sl2+cl2) | ||||
| sb2 = amin1(sinbet*sinbet, 1.) | ||||
| cosbet = sqrt(1. - sb2) | ||||
| rgeoid = a/sqrt(1. + ep2*sb2) | ||||
| x = rgeoid*cosbet + gdalt*coslat | ||||
| y = rgeoid*sinbet + gdalt*sinlat | ||||
| rkm = sqrt(x*x + y*y) | ||||
| gclat = atan2(y,x)/dtr | ||||
| return | ||||
| else if (i .eq. 2) then | ||||
| c .....geocentric to geodetic..... | ||||
| rer = rkm/a | ||||
| a2 = ((-1.4127348e-8/rer + .94339131e-8)/rer + | ||||
| 1 .33523288e-2)/rer | ||||
| a4 = (((-1.2545063e-10/rer + .11760996e-9)/rer + | ||||
| 1 .11238084e-4)/rer - .2814244e-5)/rer | ||||
| a6 = ((54.939685e-9/rer - 28.301730e-9)/rer + | ||||
| 1 3.5435979e-9)/rer | ||||
| a8 = (((320./rer - 252.)/rer + 64.)/rer - 5.) | ||||
| 1 /rer*0.98008304e-12 | ||||
| gcl = dtr*gclat | ||||
| ccl = cos(gcl) | ||||
| scl = sin(gcl) | ||||
| s2cl = 2.*scl*ccl | ||||
| c2cl = 2.*ccl*ccl - 1. | ||||
| s4cl = 2.*s2cl*c2cl | ||||
| c4cl = 2.*c2cl*c2cl - 1. | ||||
| s8cl = 2.*s4cl*c4cl | ||||
| s6cl = s2cl*c4cl + c2cl*s4cl | ||||
| dltcl = s2cl*a2 + s4cl*a4 + s6cl*a6 + s8cl*a8 | ||||
| gdlat = gclat + dltcl/dtr | ||||
| gdalt = rkm - a/sqrt(1.+ep2*scl*scl) | ||||
| return | ||||
| end if | ||||
| end | ||||
| c | ||||
| subroutine getshc2 (iu, fspec, nmax, erad, gh, ier) | ||||
| c =============================================================== | ||||
| c | ||||
| c version 1.01 | ||||
| c | ||||
| c reads spherical harmonic coefficients from the specified | ||||
| c file into an array. | ||||
| c | ||||
| c input: | ||||
| c iu - logical unit number | ||||
| c fspec - file specification | ||||
| c | ||||
| c output: | ||||
| c nmax - maximum degree and order of model | ||||
| c erad - earth's radius associated with the spherical | ||||
| c harmonic coefficients, in the same units as | ||||
| c elevation | ||||
| c gh - schmidt quasi-normal internal spherical | ||||
| c harmonic coefficients | ||||
| c ier - error number: = 0, no error | ||||
| c = -2, records out of order | ||||
| c = fortran run-time error number | ||||
| c | ||||
| c a. zunde | ||||
| c usgs, ms 964, box 25046 federal center, denver, co 80225 | ||||
| c | ||||
| c =============================================================== | ||||
| character fspec*(*) | ||||
| dimension gh(*) | ||||
| c --------------------------------------------------------------- | ||||
| c open coefficient file. read past first header record. | ||||
| c read degree and order of model and earth's radius. | ||||
| c --------------------------------------------------------------- | ||||
| open (iu, file=fspec, status='old', iostat=ier, err=999 ) | ||||
| read (iu, *, iostat=ier, err=999) | ||||
| read (iu, *, iostat=ier, err=999) nmax, erad | ||||
| c --------------------------------------------------------------- | ||||
| c read the coefficient file, arranged as follows: | ||||
| c | ||||
| c n m g h | ||||
| c ---------------------- | ||||
| c / 1 0 gh(1) - | ||||
| c / 1 1 gh(2) gh(3) | ||||
| c / 2 0 gh(4) - | ||||
| c / 2 1 gh(5) gh(6) | ||||
| c nmax*(nmax+3)/2 / 2 2 gh(7) gh(8) | ||||
| c records \ 3 0 gh(9) - | ||||
| c \ . . . . | ||||
| c \ . . . . | ||||
| c nmax*(nmax+2) \ . . . . | ||||
| c elements in gh \ nmax nmax . . | ||||
| c | ||||
| c n and m are, respectively, the degree and order of the | ||||
| c coefficient. | ||||
| c --------------------------------------------------------------- | ||||
| i = 0 | ||||
| C for nn = 1, nmax | ||||
| DO 1010 nn = 1, nmax | ||||
| C for mm = 0, nn | ||||
| DO 1020 mm = 0, nn | ||||
| read (iu, *, iostat=ier, err=999) n, m, g, h | ||||
| c write(*,*) n,m,g,h,nn,nmax,mm,ier | ||||
| if (nn .ne. n .or. mm .ne. m) then | ||||
| ier = -2 | ||||
| goto 999 | ||||
| endif | ||||
| i = i + 1 | ||||
| gh(i) = g | ||||
| if (m .ne. 0) then | ||||
| i = i + 1 | ||||
| gh(i) = h | ||||
| endif | ||||
| C end for | ||||
| 1020 CONTINUE | ||||
| 1021 CONTINUE | ||||
| C end for | ||||
| 1010 CONTINUE | ||||
| 1011 CONTINUE | ||||
| 999 close (iu) | ||||
| return | ||||
| end | ||||
| c | ||||
| c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | ||||
| c | ||||
| subroutine interpshc (date, dte1, nmax1, gh1, dte2, | ||||
| 1 nmax2, gh2, nmax, gh) | ||||
| c =============================================================== | ||||
| c | ||||
| c version 1.01 | ||||
| c | ||||
| c interpolates linearly, in time, between two spherical | ||||
| c harmonic models. | ||||
| c | ||||
| c input: | ||||
| c date - date of resulting model (in decimal year) | ||||
| c dte1 - date of earlier model | ||||
| c nmax1 - maximum degree and order of earlier model | ||||
| c gh1 - schmidt quasi-normal internal spherical | ||||
| c harmonic coefficients of earlier model | ||||
| c dte2 - date of later model | ||||
| c nmax2 - maximum degree and order of later model | ||||
| c gh2 - schmidt quasi-normal internal spherical | ||||
| c harmonic coefficients of later model | ||||
| c | ||||
| c output: | ||||
| c gh - coefficients of resulting model | ||||
| c nmax - maximum degree and order of resulting model | ||||
| c | ||||
| c a. zunde | ||||
| c usgs, ms 964, box 25046 federal center, denver, co 80225 | ||||
| c | ||||
| c =============================================================== | ||||
| dimension gh1(*), gh2(*), gh(*) | ||||
| c --------------------------------------------------------------- | ||||
| c the coefficients (gh) of the resulting model, at date | ||||
| c date, are computed by linearly interpolating between the | ||||
| c coefficients of the earlier model (gh1), at date dte1, | ||||
| c and those of the later model (gh2), at date dte2. if one | ||||
| c model is smaller than the other, the interpolation is | ||||
| c performed with the missing coefficients assumed to be 0. | ||||
| c --------------------------------------------------------------- | ||||
| factor = (date - dte1) / (dte2 - dte1) | ||||
| if (nmax1 .eq. nmax2) then | ||||
| k = nmax1 * (nmax1 + 2) | ||||
| nmax = nmax1 | ||||
| else if (nmax1 .gt. nmax2) then | ||||
| k = nmax2 * (nmax2 + 2) | ||||
| l = nmax1 * (nmax1 + 2) | ||||
| C for i = k + 1, l | ||||
| DO 1010 i = k + 1, l | ||||
| gh(i) = gh1(i) + factor * (-gh1(i)) | ||||
| C end for | ||||
| 1010 CONTINUE | ||||
| 1011 CONTINUE | ||||
| nmax = nmax1 | ||||
| else | ||||
| k = nmax1 * (nmax1 + 2) | ||||
| l = nmax2 * (nmax2 + 2) | ||||
| C for i = k + 1, l | ||||
| DO 1020 i = k + 1, l | ||||
| gh(i) = factor * gh2(i) | ||||
| C end for | ||||
| 1020 CONTINUE | ||||
| 1021 CONTINUE | ||||
| nmax = nmax2 | ||||
| endif | ||||
| C for i = 1, k | ||||
| DO 1030 i = 1, k | ||||
| gh(i) = gh1(i) + factor * (gh2(i) - gh1(i)) | ||||
| C end for | ||||
| 1030 CONTINUE | ||||
| 1031 CONTINUE | ||||
| return | ||||
| end | ||||
| c | ||||
| c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | ||||
| c | ||||
| subroutine extrapshc (date, dte1, nmax1, gh1, nmax2, | ||||
| 1 gh2, nmax, gh) | ||||
| c =============================================================== | ||||
| c | ||||
| c version 1.01 | ||||
| c | ||||
| c extrapolates linearly a spherical harmonic model with a | ||||
| c rate-of-change model. | ||||
| c | ||||
| c input: | ||||
| c date - date of resulting model (in decimal year) | ||||
| c dte1 - date of base model | ||||
| c nmax1 - maximum degree and order of base model | ||||
| c gh1 - schmidt quasi-normal internal spherical | ||||
| c harmonic coefficients of base model | ||||
| c nmax2 - maximum degree and order of rate-of-change | ||||
| c model | ||||
| c gh2 - schmidt quasi-normal internal spherical | ||||
| c harmonic coefficients of rate-of-change model | ||||
| c | ||||
| c output: | ||||
| c gh - coefficients of resulting model | ||||
| c nmax - maximum degree and order of resulting model | ||||
| c | ||||
| c a. zunde | ||||
| c usgs, ms 964, box 25046 federal center, denver, co 80225 | ||||
| c | ||||
| c =============================================================== | ||||
| dimension gh1(*), gh2(*), gh(*) | ||||
| c --------------------------------------------------------------- | ||||
| c the coefficients (gh) of the resulting model, at date | ||||
| c date, are computed by linearly extrapolating the coef- | ||||
| c ficients of the base model (gh1), at date dte1, using | ||||
| c those of the rate-of-change model (gh2), at date dte2. if | ||||
| c one model is smaller than the other, the extrapolation is | ||||
| c performed with the missing coefficients assumed to be 0. | ||||
| c --------------------------------------------------------------- | ||||
| factor = (date - dte1) | ||||
| if (nmax1 .eq. nmax2) then | ||||
| k = nmax1 * (nmax1 + 2) | ||||
| nmax = nmax1 | ||||
| else if (nmax1 .gt. nmax2) then | ||||
| k = nmax2 * (nmax2 + 2) | ||||
| l = nmax1 * (nmax1 + 2) | ||||
| C for i = k + 1, l | ||||
| DO 1010 i = k + 1, l | ||||
| gh(i) = gh1(i) | ||||
| C end for | ||||
| 1010 CONTINUE | ||||
| 1011 CONTINUE | ||||
| nmax = nmax1 | ||||
| else | ||||
| k = nmax1 * (nmax1 + 2) | ||||
| l = nmax2 * (nmax2 + 2) | ||||
| C for i = k + 1, l | ||||
| DO 1020 i = k + 1, l | ||||
| gh(i) = factor * gh2(i) | ||||
| C end for | ||||
| 1020 CONTINUE | ||||
| 1021 CONTINUE | ||||
| nmax = nmax2 | ||||
| endif | ||||
| C for i = 1, k | ||||
| DO 1030 i = 1, k | ||||
| gh(i) = gh1(i) + factor * gh2(i) | ||||
| C end for | ||||
| 1030 CONTINUE | ||||
| 1031 CONTINUE | ||||
| return | ||||
| end | ||||
| c | ||||
| c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | ||||
| c | ||||
| subroutine shval3 (igdgc, flat, flon, elev, erad, a2, b2, | ||||
| 1 nmax, gh, iext, ext, x, y, z) | ||||
| c ================================================================ | ||||
| c | ||||
| c version 1.01 | ||||
| c | ||||
| c calculates field components from spherical harmonic (sh) | ||||
| c models. | ||||
| c | ||||
| c input: | ||||
| c igdgc - indicates coordinate system used; set equal to | ||||
| c 1 if geodetic, 2 if geocentric | ||||
| c flat - north latitude, in degrees | ||||
| c flon - east longitude, in degrees | ||||
| c elev - elevation above mean sea level (igdgc=1), or | ||||
| c radial distance from earth's center (igdgc=2) | ||||
| c erad - value of earth's radius associated with the sh | ||||
| c coefficients, in same units as elev | ||||
| c a2,b2 - squares of semi-major and semi-minor axes of | ||||
| c the reference spheroid used for transforming | ||||
| c between geodetic and geocentric coordinates or | ||||
| c components | ||||
| c nmax - maximum degree and order of coefficients | ||||
| c gh - schmidt quasi-normal internal spherical | ||||
| c harmonic coefficients | ||||
| c iext - external coefficients flag (= 0 if none) | ||||
| c ext - the three 1st-degree external coefficients | ||||
| c (not used if iext = 0) | ||||
| c | ||||
| c output: | ||||
| c x - northward component | ||||
| c y - eastward component | ||||
| c z - vertically-downward component | ||||
| 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 ================================================================ | ||||
| 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. (note that this | ||||
| c version is dimensioned for nmax of 14 or less). | ||||
| c | ||||
| c minimum dimension | ||||
| c -------------------------- | ||||
| c nmax | ||||
| dimension sl(14), cl(14) | ||||
| c (nmax * (nmax + 3)) / 2 | ||||
| dimension p(119), q(119) | ||||
| c nmax * (nmax + 2) | ||||
| dimension gh(224) | ||||
| c 3 | ||||
| dimension ext(3) | ||||
| c ================================================================ | ||||
| parameter (dtr = .01745329) | ||||
| r = elev | ||||
| slat = sin (flat * dtr) | ||||
| if (90. - flat .lt. .001) then | ||||
| c 300 ft from n. pole | ||||
| aa = 89.999 | ||||
| else if (90. + flat .lt. .001) then | ||||
| c 300 ft from s. pole | ||||
| aa = -89.999 | ||||
| else | ||||
| aa = flat | ||||
| endif | ||||
| clat = cos (aa * dtr) | ||||
| sl(1) = sin (flon * dtr) | ||||
| cl(1) = cos (flon * dtr) | ||||
| x = 0. | ||||
| y = 0. | ||||
| z = 0. | ||||
| sd = 0. | ||||
| cd = 1. | ||||
| n = 0 | ||||
| l = 1 | ||||
| m = 1 | ||||
| npq = (nmax * (nmax + 3)) / 2 | ||||
| if (igdgc .eq. 1) then | ||||
| aa = a2 * clat * clat | ||||
| bb = b2 * slat * slat | ||||
| cc = aa + bb | ||||
| dd = sqrt (cc) | ||||
| r = sqrt (elev * (elev + 2. * dd) | ||||
| 1 + (a2 * aa + b2 * bb) / cc) | ||||
| cd = (elev + dd) / r | ||||
| sd = (a2 - b2) / dd * slat * clat / r | ||||
| aa = slat | ||||
| slat = slat * cd - clat * sd | ||||
| clat = clat * cd + aa * sd | ||||
| endif | ||||
| ratio = erad / r | ||||
| aa = sqrt (3.) | ||||
| p(1) = 2. * slat | ||||
| p(2) = 2. * clat | ||||
| p(3) = 4.5 * slat * slat - 1.5 | ||||
| p(4) = 3. * aa * clat * slat | ||||
| q(1) = -clat | ||||
| q(2) = slat | ||||
| q(3) = -3. * clat * slat | ||||
| q(4) = aa * (slat * slat - clat * clat) | ||||
| C for k = 1, npq | ||||
| DO 1010 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. - .5 / fm) | ||||
| j = k - n - 1 | ||||
| p(k) = (1. + 1. / 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 * fn - fm * fm) | ||||
| bb = sqrt ((fn - 1.)**2 - fm * fm) / aa | ||||
| cc = (2. * fn - 1.) / aa | ||||
| i = k - n | ||||
| j = k - 2 * n + 1 | ||||
| p(k) = (fn + 1.) * (cc * slat / fn * p(i) - bb | ||||
| 1 / (fn - 1.) * p(j)) | ||||
| q(k) = cc * (slat * q(i) - clat / fn * p(i)) | ||||
| 1 - 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.) then | ||||
| y = y + (aa * sl(m) - bb * cl(m)) * fm * p(k) | ||||
| 1 / ((fn + 1.) * clat) | ||||
| else | ||||
| y = y + (aa * sl(m) - bb * cl(m)) * q(k) | ||||
| 1 * slat | ||||
| endif | ||||
| l = l + 2 | ||||
| endif | ||||
| m = m + 1 | ||||
| C end for | ||||
| 1010 CONTINUE | ||||
| 1011 CONTINUE | ||||
| if (iext .ne. 0) then | ||||
| aa = ext(2) * cl(1) + ext(3) * sl(1) | ||||
| x = x - ext(1) * clat + aa * slat | ||||
| y = y + ext(2) * sl(1) - ext(3) * cl(1) | ||||
| z = z + ext(1) * slat + aa * clat | ||||
| endif | ||||
| aa = x | ||||
| x = x * cd + z * sd | ||||
| z = z * cd - aa * sd | ||||
| return | ||||
| end | ||||
| c | ||||
| c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | ||||
| c | ||||
| subroutine dihf (x, y, z, d, i, h, f) | ||||
| c =============================================================== | ||||
| c | ||||
| c version 1.01 | ||||
| c | ||||
| c computes the geomagnetic elements d, i, h, and f from | ||||
| c x, y, and z. | ||||
| c | ||||
| c input: | ||||
| c x - northward component | ||||
| c y - eastward component | ||||
| c z - vertically-downward component | ||||
| c | ||||
| c output: | ||||
| c d - declination | ||||
| c i - inclination | ||||
| c h - horizontal intensity | ||||
| c f - total intensity | ||||
| c | ||||
| c a. zunde | ||||
| c usgs, ms 964, box 25046 federal center, denver, co 80225 | ||||
| c | ||||
| c =============================================================== | ||||
| real i | ||||
| parameter ( sn = 0.0001 ) | ||||
| c --------------------------------------------------------------- | ||||
| c if d and i cannot be determined, set equal to 999.0. | ||||
| c --------------------------------------------------------------- | ||||
| h2 = x*x + y*y | ||||
| h = sqrt (h2) | ||||
| f = sqrt (h2 + z*z) | ||||
| if (f .lt. sn) then | ||||
| d = 999. | ||||
| i = 999. | ||||
| else | ||||
| i = atan2 (z, h) | ||||
| if (h .lt. sn) then | ||||
| d = 999. | ||||
| else | ||||
| hpx = h + x | ||||
| if (hpx .lt. sn) then | ||||
| d = 180. | ||||
| else | ||||
| d = 2. * atan2 (y, hpx) | ||||
| endif | ||||
| endif | ||||
| endif | ||||
| return | ||||
| end | ||||
