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