igrf10-02.f


      PROGRAM IGRF10
C
C     This is a program for synthesising geomagnetic field values from the 
C     International Geomagnetic Reference Field series of models as agreed
c     in Decmember 2004 by IAGA Working Group V-MOD. 
C     It is the 10th generation IGRF, ie the 9th revision. 
C     The main-field models for 1900.0, 1905.0,..1940.0 and 2005.0 are 
C     non-definitive, those for 1945.0, 1950.0,...2000.0 are definitive and
C     the secular-variation model for 2005.0 to 2010.0 is non-definitive.
C
C     Main-field models are to degree and order 10 (ie 120 coefficients)
C     for 1900.0-1995.0 and to 13 (ie 195 coefficients) for 2000.0 onwards. 
C     The predictive secular-variation model is to degree and order 8 (ie 80
C     coefficients).
C
C     Options include values at different locations at different
C     times (spot), values at same location at one year intervals
C     (time series), grid of values at one time (grid); geodetic or
C     geocentric coordinates, latitude & longitude entered as decimal
C     degrees or degrees & minutes (not in grid), choice of main field 
C     or secular variation or both (grid only).
C
c     Adapted from 8th generation version to include new maximum degree for
c     main-field models for 2000.0 and onwards and use WGS84 spheroid instead
c     of International Astronomical Union 1966 spheroid as recommended by IAGA
c     in July 2003. Reference radius remains as 6371.2 km - it is NOT the mean
c     radius (= 6371.0 km) but 6371.2 km is what is used in determining the
c     coefficients. Adaptation by Susan Macmillan, August 2003 (for 
c     9th generation) and December 2004.
c     1995.0 coefficients as published in igrf9coeffs.xls and igrf10coeffs.xls
c     used - (Kimmo Korhonen spotted 1 nT difference in 11 coefficients)
c     Susan Macmillan July 2005
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*11 TYPE
      dimension wwlat(30000),wwlon(30000)
      dimension yr0(30000),gr0(30000),gto0(30000),ga0(30000)
      dimension iyr1(30000),jdy1(30000),min1(30000)
      dimension speed(30000),head(30000),kdep(30000),kflg(30000)
      dimension tmag(30000),nn(30000)

      dtmn=1900.0
      dtmx=2015.0
C
      WRITE(6,*)
      WRITE(6,*)'******************************************************'
      WRITE(6,*)'*              IGRF SYNTHESIS PROGRAM                *'
      WRITE(6,*)'*                                                    *'
      WRITE(6,*)'* A program for the computation of geomagnetic       *'
      WRITE(6,*)'* field elements from the International Geomagnetic  *'
      WRITE(6,*)'* Reference Field (9th generation) as revised in     *'
      WRITE(6,*)'* December 2004 by the IAGA Working Group V-MOD.     *'
      WRITE(6,*)'*                                                    *'
      WRITE(6,*)'* It is valid for dates from 1900.0 to 2010.0,       *'
      WRITE(6,*)'* values up to 2015.0 will be computed but with      *'
      WRITE(6,*)'* reduced accuracy. Values for dates before 1945.0   *'
      WRITE(6,*)'* and after 2000.0 are non-definitive, otherwise the *'
      WRITE(6,*)'* values are definitive.                             *'
      WRITE(6,*)'*                                                    *'
      WRITE(6,*)'* Susan Macmillan          British Geological Survey *'
      WRITE(6,*)'*                     Chair IAGA Working Group V-MOD *'
      WRITE(6,*)'******************************************************'
      WRITE(6,*)
c
      open(20,file='H:\Myhome\gh90a.bs1',
     +   status='old')
      open(30,file='H:\Myhome\gh90a.bs2')
c     if the year is URU, ik=1;others ik=0
      ik=0
c      
      FACT = 180.0/3.141592654
      NCOUNT = 0
C
      itype=1
C
      iopt=1
      idm=2
c      
      i=0
    1 continue
      read(20,*,end=99) no,iyr,jdy,min,alt,aln,
     +   spd0,hdg0,idp,mgt,mga,
     +   gr,gto,ga,iflg
c      write(6,*) no
           i=i+1
           nn(i)=no
           iyr1(i)=iyr
           jdy1(i)=jdy
           min1(i)=min
           wwlat(i)=dble(alt)
           wwlon(i)=dble(aln)
           speed(i)=spd0
           head(i)=hdg0
           kdep(i)=idp
           tmag(i)=float(mgt)
           gr0(i)=gr
           gto0(i)=gto
           ga0(i)=ga
           kflg(i)=iflg
           iyr=iyr+1900
           yr0(i)=float(iyr)+float(jdy)/float(365+ik)
         go to 1
   99 close(20)
      jnum=i
      write(6,*) 'read end-30',jnum
      k=0
C
   50 ncount=1
      ALT=0
C
      do 235 is=1,jnum
       xlt=wwlat(is)
       xln=wwlon(is)
       date=yr0(is)
       IF (XLT.LT.-90.0.OR.XLT.GT.90.0) then
          GO TO 202
       end if
       IF (XLN.LT.-360.0.OR.XLN.GT.360.0) then
         GO TO 203
       ENDIF
C
      CLT = 90.0 - XLT
      IF (CLT.LT.0.0.OR.CLT.GT.180.0) GO TO 204
      IF (XLN.LE.-360.0.OR.XLN.GE.360.0) GO TO 205
C
      CALL IGRF10SYN (0,DATE,ITYPE,ALT,CLT,XLN,X,Y,Z,F)
      D = FACT*ATAN2(Y,X)
      H = SQRT(X*X + Y*Y)
      S = FACT*ATAN2(Z,H)
      CALL DDECDM (D,IDEC,IDECM)
      CALL DDECDM (S,INC,INCM)
C
      CALL IGRF10SYN (1,DATE,ITYPE,ALT,CLT,XLN,DX,DY,DZ,F1)
      DD = (60.0*FACT*(X*DY - Y*DX))/(H*H)
      DH = (X*DX + Y*DY)/H
      DS = (60.0*FACT*(H*DZ - Z*DH))/(F*F)
      DF = (H*DH + Z*DZ)/F
C
      amag=tmag(is)-f
c
      write(6,*) is,iyr1(is),jdy1(is),min1(is),tmag(is),amag
c
      write(30,207)  nn(is),iyr1(is),jdy1(is),min1(is),xlt,xln,
     +   speed(is),head(is),kdep(is),tmag(is),amag,
     +   gr0(is),gto0(is),ga0(is),kflg(is)
  207 format(i7,i3,i4,i5,f10.5,f11.5,f7.2,f6.1,i6,f9.1,f8.1,f7.1,
     +   f10.1,f8.1,i2)
c
C
  235 continue
C
   60 CONTINUE
C
C
  159 continue
      close(30)
      STOP
C
  209 WRITE(6,972) DATE
  972 FORMAT (' ***** Error *****'/' DATE =',F9.3,
     1        ' - out of range')
      STOP
C
  210 WRITE(6,973) ALT,ITYPE
  973 FORMAT (' ***** Error *****'/' A value of ALT =',F10.3,
     1        ' is not allowed when ITYPE =',I2)
      STOP
C
  202 WRITE(6,966) XLT
  966 FORMAT (' ***** Error *****'/' XLT =',F9.3,
     1        ' - out of range')
      STOP
C
  203 WRITE(6,967) XLN
  967 FORMAT (' ***** Error *****'/' XLN =',F10.3,
     1        ' - out of range')
      STOP
C
  204 WRITE(6,968) LTD,LTM
  968 FORMAT (' ***** Error *****'/' Latitude out of range',
     1        ' - LTD =',I6,5X,'LTM =',I4)
      STOP
C
  205 WRITE(6,969) LND,LNM
  969 FORMAT (' ***** Error *****'/' Longitude out of range',
     1        ' - LND =',I8,5X,'LNM =',I4)
      STOP
C
      END
C
      SUBROUTINE DMDDEC (I,M,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DE = I
      EM = M
      IF (I.LT.0) EM = -EM
      X = DE + EM/60.0
      RETURN
      END
C
      SUBROUTINE DDECDM (X,I,M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      if(x.ge.0) then
         sig=1.1
       else
         sig=-1.1
      end if
c     SIG = SIGN(1.1,X)
      DR = ABS(X)
      I = INT(DR)
      T = I
      M = NINT(60.*(DR - T))
      IF (M.EQ.60) THEN
       M = 0
       I = I + 1
      ENDIF
      ISIG = INT(SIG)
      IF (I.NE.0) THEN
       I = I * ISIG
      ELSE
       IF (M.NE.0) M = M * ISIG
      ENDIF
      RETURN
      END
c
      subroutine igrf10syn (isv,date,itype,alt,colat,elong,x,y,z,f)
c
c     This is a synthesis routine for the 10th generation IGRF as agreed 
c     in December 2004 by IAGA Working Group V-MOD. It is valid 1900.0 to
c     2010.0 inclusive. Values for dates from 1945.0 to 2000.0 inclusive are 
c     definitve, otherwise they are non-definitive.
c   INPUT
c     isv   = 0 if main-field values are required
c     isv   = 1 if secular variation values are required
c     date  = year A.D. Must be greater than or equal to 1900.0 and 
c             less than or equal to 2015.0. Warning message is given 
c             for dates greater than 2010.0. Must be double precision.
c     itype = 1 if geodetic (spheroid)
c     itype = 2 if geocentric (sphere)
c     alt   = height in km above sea level if itype = 1
c           = distance from centre of Earth in km if itype = 2 (>3485 km)
c     colat = colatitude (0-180)
c     elong = east-longitude (0-360)
c     alt, colat and elong must be double precision.
c   OUTPUT
c     x     = north component (nT) if isv = 0, nT/year if isv = 1
c     y     = east component (nT) if isv = 0, nT/year if isv = 1
c     z     = vertical component (nT) if isv = 0, nT/year if isv = 1
c     f     = total intensity (nT) if isv = 0, rubbish if isv = 1
c
c     To get the other geomagnetic elements (D, I, H and secular
c     variations dD, dH, dI and dF) use routines ptoc and ptocsv.
c
c     Adapted from 8th generation version to include new maximum degree for
c     main-field models for 2000.0 and onwards and use WGS84 spheroid instead
c     of International Astronomical Union 1966 spheroid as recommended by IAGA
c     in July 2003. Reference radius remains as 6371.2 km - it is NOT the mean
c     radius (= 6371.0 km) but 6371.2 km is what is used in determining the
c     coefficients. Adaptation by Susan Macmillan, August 2003 (for 
c     9th generation) and December 2004.
c     1995.0 coefficients as published in igrf9coeffs.xls and igrf10coeffs.xls
c     used - (Kimmo Korhonen spotted 1 nT difference in 11 coefficients)
c     Susan Macmillan July 2005
c
      implicit double precision (a-h,o-z)
      dimension gh(3060),g0(120),g1(120),g2(120),g3(120),g4(120),
     1          g5(120),g6(120),g7(120),g8(120),g9(120),ga(120),
     2          gb(120),gc(120),gd(120),ge0(120),gf(120),gg(120),
     3          gi(120),gj(120),gk(195),gl(195),gm(195),gp(195),
     4          p(105),q(105),cl(13),sl(13)
      equivalence (g0,gh(  1)),(g1,gh(121)),(g2,gh(241)),(g3,gh(361)),
     1            (g4,gh(481)),(g5,gh(601)),(g6,gh(721)),(g7,gh(841)),
     2            (g8,gh(961)),(g9,gh(1081)),(ga,gh(1201)),
     3            (gb,gh(1321)),(gc,gh(1441)),(gd,gh(1561)),
     4            (ge0,gh(1681)),(gf,gh(1801)),(gg,gh(1921)),
     5            (gi,gh(2041)),(gj,gh(2161)),(gk,gh(2281)),
     6            (gl,gh(2476)),(gm,gh(2671)),(gp,gh(2866))
c
      open(10,file='H:\Myhome\programs\igrf10coeffs.data',status='old')
      do 120 kkk=1,120
      read(10,*) g0(kkk)
  120 continue
      do 121 kkk=1,120
      read(10,*) g1(kkk)
  121 continue
      do 122 kkk=1,120
      read(10,*) g2(kkk)
  122 continue
      do 123 kkk=1,120
      read(10,*) g3(kkk)
  123 continue
      do 124 kkk=1,120
      read(10,*) g4(kkk)
  124 continue
      do 125 kkk=1,120
      read(10,*) g5(kkk)
  125 continue
      do 126 kkk=1,120
      read(10,*) g6(kkk)
  126 continue
      do 127 kkk=1,120
      read(10,*) g7(kkk)
  127 continue
      do 128 kkk=1,120
      read(10,*) g8(kkk)
  128 continue
      do 129 kkk=1,120
      read(10,*) g9(kkk)
  129 continue
      do 130 kkk=1,120
      read(10,*) ga(kkk)
  130 continue
      do 131 kkk=1,120
      read(10,*) gb(kkk)
  131 continue
      do 132 kkk=1,120
      read(10,*) gc(kkk)
  132 continue
      do 133 kkk=1,120
      read(10,*) gd(kkk)
  133 continue
      do 134 kkk=1,120
      read(10,*) ge0(kkk)
  134 continue
      do 135 kkk=1,120
      read(10,*) gf(kkk)
  135 continue
      do 136 kkk=1,120
      read(10,*) gg(kkk)
  136 continue
      do 137 kkk=1,120
      read(10,*) gi(kkk)
  137 continue
      do 138 kkk=1,120
      read(10,*) gj(kkk)
  138 continue
      do 139 kkk=1,120
      read(10,*) gk(kkk)
  139 continue
      do 140 kkk=1,195
      read(10,*) gl(kkk)
  140 continue
       do 141 kkk=1,195
      read(10,*) gm(kkk)
  141 continue
       do 142 kkk=1,80
      read(10,*) gp(kkk)
  142 continue
      do 143 lll=121, 195
      gk(lll)=0
  143 continue
      do 144 lll=81, 195
      gp(lll)=0.0
  144 continue
c
      close(10)
c     write(6,*) gp(1),g0(1)
c
c
c     set initial values
c
      x     = 0.0
      y     = 0.0
      z     = 0.0
      if (date.lt.1900.0.or.date.gt.2015.0) go to 11
      if (date.gt.2010.0) write (6,960) date
  960 format (/' This version of the IGRF is intended for use up',
     1        ' to 2010.0.'/' values for',f9.3,' will be computed',
     2        ' but may be of reduced accuracy'/)
      if (date.ge.2005.0) go to 1
      t     = 0.2*(date - 1900.0)                                             
      ll    = t
      one   = ll
      t     = t - one
      if (date.lt.1995.0) then
       nmx   = 10
       nc    = nmx*(nmx+2)
       ll    = nc*ll
       kmx   = (nmx+1)*(nmx+2)/2
      else
       nmx   = 13
       nc    = nmx*(nmx+2)
       ll    = 0.2*(date - 1995.0)
       ll    = 120*19 + nc*ll
       kmx   = (nmx+1)*(nmx+2)/2
      endif
      tc    = 1.0 - t
      if (isv.eq.1) then
       tc = -0.2
       t = 0.2
      end if
      go to 2
c
    1 t     = date - 2005.0
      tc    = 1.0
      if (isv.eq.1) then
       t = 1.0
       tc = 0.0
      end if
      ll    = 2670
      nmx   = 13
      nc    = nmx*(nmx+2)
      kmx   = (nmx+1)*(nmx+2)/2
    2 r     = alt
      one   = colat*0.017453292
      ct    = cos(one)
      st    = sin(one)
      one   = elong*0.017453292
      cl(1) = cos(one)
      sl(1) = sin(one)
      cd    = 1.0
      sd    = 0.0
      l     = 1
      m     = 1
      n     = 0
      if (itype.eq.2) go to 3
c
c     conversion from geodetic to geocentric coordinates 
c     (using the WGS84 spheroid)
c
      a2    = 40680631.6
      b2    = 40408296.0
      one   = a2*st*st
      two   = b2*ct*ct
      three = one + two
      rho   = sqrt(three)
      r     = sqrt(alt*(alt + 2.0*rho) + (a2*one + b2*two)/three)
      cd    = (alt + rho)/r
      sd    = (a2 - b2)/rho*ct*st/r
      one   = ct
      ct    = ct*cd -  st*sd
      st    = st*cd + one*sd
c
    3 ratio = 6371.2/r
      rr    = ratio*ratio
c
c     computation of Schmidt quasi-normal coefficients p and x(=q)
c
      p(1)  = 1.0
      p(3)  = st
      q(1)  = 0.0
      q(3)  =  ct
      do 10 k=2,kmx
      if (n.ge.m) go to 4
      m     = 0
      n     = n + 1
      rr    = rr*ratio
      fn    = n
      gn    = n - 1
    4 fm    = m
      if (m.ne.n) go to 5
      if (k.eq.3) go to 6
      one   = sqrt(1.0 - 0.5/fm)
      j     = k - n - 1
      p(k)  = one*st*p(j)
      q(k)  = one*(st*q(j) + ct*p(j))
      cl(m) = cl(m-1)*cl(1) - sl(m-1)*sl(1)
      sl(m) = sl(m-1)*cl(1) + cl(m-1)*sl(1)
      go to 6                                                           
    5 gmm    = m*m
      one   = sqrt(fn*fn - gmm)
      two   = sqrt(gn*gn - gmm)/one
      three = (fn + gn)/one
      i     = k - n
      j     = i - n + 1
      p(k)  = three*ct*p(i) - two*p(j)
      q(k)  = three*(ct*q(i) - st*p(i)) - two*q(j)
c
c     synthesis of x, y and z in geocentric coordinates
c
    6 lm    = ll + l
      one   = (tc*gh(lm) + t*gh(lm+nc))*rr
      if (m.eq.0) go to 9
      two   = (tc*gh(lm+1) + t*gh(lm+nc+1))*rr
      three = one*cl(m) + two*sl(m)
      x     = x + three*q(k)
      z     = z - (fn + 1.0)*three*p(k)
      if (st.eq.0.0) go to 7
      y     = y + (one*sl(m) - two*cl(m))*fm*p(k)/st
      go to 8
    7 y     = y + (one*sl(m) - two*cl(m))*q(k)*ct
    8 l     = l + 2
      go to 10
    9 x     = x + one*q(k)
      z     = z - (fn + 1.0)*one*p(k)
      l     = l + 1
   10 m     = m + 1
c
c     conversion to coordinate system specified by itype
c
      one   = x
      x     = x*cd +   z*sd
      z     = z*cd - one*sd
      f     = sqrt(x*x + y*y + z*z)
c
      return
c
c     error return if date out of bounds
c
   11 f     = 1.0d8
      write (6,961) date
  961 format (/' This subroutine will not work with a date of',
     1        f9.3,'.  Date must be in the range 1900.0.ge.date',
     2        '.le.2015.0.  On return f = 1.0d8., x = y = z = 0.')
      return
      end