subroutine dgrf70_95 (year,rre,thet,phi,br,bt,bp,b,ier)
c*
c***********************************************************************
c*
c*PRO MAGLIB
c*
c*VER 99.03.31 - V 1.0
c*VER 01.06.05 - V 2.0
c*VER 03.01.06 - V 2.1
c*VER 01.01.07 - V 3.0
c*
c*AUT spec. N A S A
c*AUT spec. CNES - JC KOSIK - janvier 1996
c*AUT port. CISI
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul du champ magnetique pour la periode 
c*ROL         comprise entre 1970 et 1995.
c*
c*PAR year (I) : annee fractionnaire (entre 1970. et 1995.)
c*
c*PAR rre  (I) : distance radiale geocentrique (rayons terrestres)
c*PAR thet (I) : colatitude geocentrique (radians), la colatitude
c*PAR          : est calculee du pole nord au pole sud 
c*PAR phi  (I) : longitude geocentrique (radians), calculee 
c*PAR          : positive vers l'est de 0 a 2*PI
c*
c*PAR br   (O) : composante radiale du champ magnetique le long du
c*             : meridien positive vers l'exterieur (gauss)
c*PAR bt   (O) : composante tangentielle du champ magnetique le long
c*PAR          : du meridien positive vers le sud (gauss)
c*PAR bp   (O) : composante azimuthale du champ magnetique, positive
c*PAR          : vers l'est (gauss)
c*
c*PAR b    (O) : module du champ magnetique (gauss)
c*
c*PAR ier  (O) : code de retour
c*
c*NOT ier      : sans objet
c
c*NOT Un test est effectue pour les colatitudes inferieurs a 1.d-15
c*NOT ou superieures a PI - 1.d-15 et le cosinus est remplace par 1.d0
c*NOT et le sinus par 1.d-15 ou bien le cosinus est remplace par -1.d0
c*NOT et le sinus par 1.d-15.
c*
c*NOT common   : util, util2
c*
c*INF utilise  : sans objet
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.05 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - coefficients definitifs du champ chp95
c*HST                        - remplace dgrfnew de la version 2.1
c*
c***********************************************************************
c*
      implicit none
c
c     ---------------------------------
c*FON Declaration identificateur rcs_id
c     ---------------------------------
c
      character rcs_id*100
c
c     --------------------------
c*FON Declaration des parametres
c     --------------------------
c
      double precision year
      double precision rre, thet, phi
      double precision br, bt, bp, b
      integer ier
c
c     ----------------------------------
c*FON Declaration des variables communes
c     ----------------------------------
c
      double precision pi,dpi,rad,deg,pid,xmu,rayt
c
c*COM pi   : constante pi (obtenue a partir de acos(-1.))
c*COM dpi  : constante 2 * pi
c*COM pid  : constante pi / 2
c*COM rad  : facteur de conversion degres  ----> radians
c*COM deg  : facteur de conversion radians ----> degres
c*COM xmu  : constante de gravitation terrestre (km**3/sec**2)
c*COM rayt : rayon equatorial terrestre (km)
c
      common/util/pi,dpi,rad,deg,pid,xmu,rayt
c
      double precision rgmt
c
c*COM rgmt : rayon geocentrique terrestre (km)
c
      common/util2/rgmt
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer lg(11,11,6)
c*LOC lg : coefficients entre 1970 et 1995
c
      integer kmax
c*LOC kmax : nombre de pas de calcul (= 11)
c
      integer n,m,jj,km
c*LOC n,m,jj,km : indices de boucles et de tableaux
c
      double precision cof1,cof2,delt
      double precision shmit(11,11)
      double precision g(11,11)
      double precision const(11,11),fn(11),fm(11)
      double precision p(11,11),dp(11,11),sp(11),cp(11)
      double precision tzero(6)
      double precision tmold,sph,st,temp,par
      double precision cph,ct,aor,ar,rrmag
c*LOC Variables de travail intermediaires
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     --------------------------
c*FON Affectation des constantes
c     --------------------------
c
c*   coefficients pour 1970,1975,1980,1985,1990,1995
c
c    Revision octobre 2006, nouveaux coefficients pour le champ dgrf95
c
      data lg/
     >     1,-30220,-1781, 1287,  952, -216,   43,   72,  14,    8, -3,
     >  5737, -2068, 3000,-2091,  800,  359,   64,  -57,   6,   10, -3,
     > -2047,    25, 1611, 1278,  461,  262,   15,    1,  -2,    2,  2,
     >  -366,   251, -196,  838, -395,  -42, -212,   14, -13,  -12, -5,
     >   167,  -266,   26, -279,  234, -160,    2,  -22,  -3,   10, -1,
     >    26,   139, -139,  -91,   83,  -56,    3,   -2,   5,   -1,  6,
     >   -12,   100,   72,  -37,   -6,    1, -112,   13,   0,    0,  4,
     >   -70,   -27,   -4,    8,   23,  -23,  -11,   -2,  11,    3,  1,
     >     7,   -15,    6,  -17,    6,   21,   -6,  -16,   3,    1,  0,
     >   -21,    16,    6,   -4,   -5,   10,   11,   -2,   1,   -1,  3,
     >     1,     1,    3,    4,   -4,    0,   -1,    3,   1,   -4, -1,
     >
     >    1,-30100,-1902, 1276, 946, -218,   45,  71,  14,   7, -3,
     >  5675,-2013, 3010,-2144, 791,  356,   66, -56,   6,  10, -3,
     > -2067,  -68, 1632, 1260, 438,  264,   28,   1,  -1,   2,  2,
     >  -333,  262, -223,  830,-405,  -59, -198,  16, -12, -12, -5,
     >   191, -265,   39, -288, 216, -159,    1, -14,  -8,  10, -2,
     >    31,  148, -152,  -83,  88,  -49,    6,   0,   4,  -1,  5,
     >   -13,   99,   75,  -41,  -4,   11, -111,  12,   0,  -1,  4,
     >   -77,  -26,   -5,   10,  22,  -23,  -12,  -5,  10,   4,  1,
     >     6,  -16,    4,  -19,   6,   18,  -10, -17,   1,   1,  0,
     >   -21,   16,    7,   -4,  -5,   10,   11,  -3,   1,  -2,  3,
     >     1,    1,    3,    4,  -4,   -1,   -1,   3,   1,  -5, -1,
     >
     >     1,-29992,-1997, 1281,  938, -218,   48,  72, 18,    5, -4,
     >  5604, -1956, 3027,-2180,  782,  357,   66, -59,  6,   10, -4,
     > -2129,  -200, 1663, 1251,  398,  261,   42,   2,  0,    1,  2,
     >  -336,   271, -252,  833, -419,  -74, -192,  21, -11, -12, -5,
     >   212,  -257,   53, -297,  199, -162,    4, -12,  -7,   9, -2,
     >    46,   150, -151,  -78,   92,  -48,   14,   1,   4,  -3,  5,
     >   -15,    93,   71,  -43,   -2,   17, -108,  11,   3,  -1,  3,
     >   -82,   -27,   -5,   16,   18,  -23,  -10,  -2,   6,   7,  1,
     >     7,   -18,    4,  -22,    9,   16,  -13, -15,  -1,   2,  2,
     >   -21,    16,    9,   -5,   -6,    9,   10,  -6,   2,  -5,  3,
     >     1,     0,    3,    6,   -4,    0,   -1,   4,   0,  -6,  0,
     >
     >     1,-29873,-2072, 1296,  936, -214,   53,  74,  21,   5, -4,
     >  5500, -1905, 3044,-2208,  780,  355,   65, -62,   6,  10, -4,
     > -2197,  -306, 1687, 1247,  361,  253,   51,   3,   0,   1,  3,
     >  -310,   284, -297,  829, -424,  -93, -185,  24, -11, -12, -5,
     >   232,  -249,   69, -297,  170, -164,    4,  -6,  -9,   9, -2,
     >    47,   150, -154,  -75,   95,  -46,   16,   4,   4,  -3,  5,
     >   -16,    88,   69,  -48,   -1,   21, -102,  10,   4,  -1,  3,
     >   -83,   -27,   -2,   20,   17,  -23,   -7,   0,   4,   7,  1,
     >     8,   -19,    5,  -23,   11,   14,  -15, -11,  -4,   1,  2,
     >   -21,    15,    9,   -6,   -6,    9,    9,  -7,   2,  -5,  3,
     >     1,     0,    3,    6,   -4,    0,   -1,   4,   0,  -6,  0,
     >
     >     1,-29775,-2131, 1314,  939, -214,   61,  77,  23,   4, -3,
     >  5406, -1848, 3059,-2239,  780,  353,   65, -64,   5,   9, -4,
     > -2279,  -373, 1686, 1248,  325,  245,   59,   2,  -1,   1,  2,
     >  -284,   293, -352,  802, -423, -109, -178,  26, -10, -12, -5,
     >   247,  -240,   84, -299,  141, -165,    3,  -1, -12,   9, -2,
     >    46,   154, -153,  -69,   97,  -36,   18,   5,   3,  -4,  4,
     >   -16,    82,   69,  -52,    1,   24,  -96,   9,   4,  -2,  3,
     >   -80,   -26,    0,   21,   17,  -23,   -4,   0,   2,   7,  1,
     >    10,   -19,    6,  -22,   12,   12,  -16, -10,  -6,   1,  3,
     >   -20,    15,   11,   -7,   -7,    9,    8,  -7,   2,  -6,  3,
     >     2,     1,    3,    6,   -4,    0,   -2,   3,  -1,  -6,  0,
     >
     s 1,-29692,-2200,1335,940,-214,68,77,25,4,-3,
     s 5306,-1784,3070,-2267,780,352,67,-72,6,9,-6,
     s -2366,-413,1681,1249,290,235,68,1,-6,3,2,
     s -262,302,-427,759,-418,-118,-170,28,-9,-10,-4,
     s 262,-236,97,-306,122,-166,-1,5,-14,8,-1,
     s 46,165,-143,-55,107,-17,19,4,9,-8,4,
     s -17,72,67,-58,1,36,-93,8,6,-1,2,
     s -69,-25,4,24,17,-24,-6,-2,-5,10,2,
     s 11,-21,8,-23,15,11,-16,-4,-7,-2,5,
     s -20,15,12,-6,-8,8,5,-8,3,-8,1,
     s  1,0,4,5,-5,-1,-2,1,-2,-7,0/
c
      data kmax /11/
c
      data sp(1)   /0.d0/
      data p(1,1)  /1.d0/
      data cp(1)   /1.d0/
      data dp(1,1) /0.d0/
c
      data shmit(1,1) /0.d0/
c
      data tzero /1970.d0, 1975.d0, 1980.d0, 1985.d0,
     >            1990.d0, 1995.d0/
c
      data tmold /0.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
c     -------------------------------------------
c*FON Initialisation a 0 des composantes du champ
c*FON ainsi que du module
c     -------------------------------------------
c
      br = 0.0d0
      bt = 0.0d0
      bp = 0.0d0
      b  = 0.0d0
c
c     -------------------------------------------------
c*FON Initialisations faites au premier appel du module
c     -------------------------------------------------
c
      if (shmit(1,1) .ne. -1.d0) then
c
c     -----------------------------------
c*FON  Calcul des tableaux fn,fm et const
c     -----------------------------------
c
         shmit(1,1) = -1.d0
         do 10 n = 1, 11
            fn(n) = n
            do 20 m = 1, 11
               fm(m)      = m - 1
               const(n,m) = dble((n - 2)**2 - (m - 1)**2) /
     >                      dble((2 * n - 3) * (2 * n - 5))
20         continue
10       continue
c
c     ----------------------------
c*FON Calcul du tableau shmit(n,m)
c     ----------------------------
c
         do 30 n = 2, 11
            shmit(n,1) = dble(2 * n - 3) * shmit(n-1,1) 
     >                   / dble(n - 1)
            jj = 2
            do 40 m = 2, n
               shmit(n,m)   = shmit(n,m-1)
     >                        * sqrt(dble((n - m + 1) * jj) /
     >                        (n + m - 2))
               shmit(m-1,n) = shmit(n,m)
               jj = 1
40          continue
30       continue
c
c     -------------------------------------------
c*FON Traitements effectues si c'est le 1er appel
c     -------------------------------------------
c
      endif
c
c     --------------------------------------------------
c*FON Calcul du g(n,m), interpolation entre deux jeux de
c*FON coefficients avec un delai de 5 ans entre eux
c     --------------------------------------------------
c
      if (year .ne. tmold) then
c
         tmold = year
c
         km = idint((year - 1970.d0) / 5.d0) + 1
c
         delt = (year - tzero(km)) / 5.d0
c
         do 70 n = 1, kmax
            do 80 m = 1, kmax
            cof1   = dble(lg(n,m,km))
            cof2   = dble(lg(n,m,km+1))
                g(n,m) = cof1 + (cof2 - cof1) * delt
80          continue
70       continue
c
         do 77 n = 1, kmax
            do 78 m= 1, kmax
               g(n,m) = g(n,m) * shmit(n,m)
78          continue
77       continue
c
      endif
c
c     ----------------------------------------------------
c*FON Calcul des cosinus directeurs de la colatitude et de
c*FON la longitude. On test si la colatitude est comprise
c*FON entre 0 et PI
c     ----------------------------------------------------
c
      cph = cos(phi)
      sph = sin(phi)
c
      if (thet .ne. 0.0d0 .and. thet .ne. acos(-1.0d0)) then
c
         ct = cos(thet)
         st = sin(thet)
c
      else
c
         if (thet .ne. acos(-1.0d0)) then
            ct = 1.0d0
            st = 1.d-15
         elseif (thet .ne. 0.0d0) then
            ct = -1.0d0
            st = 1.d-15 
         endif
c
      endif
c
      sp(2) = sph
      cp(2) = cph
c
c     ----------------------------
c*FON Calcul des tableaux cp et sp
c     ----------------------------
c
      do 90 m = 3, kmax
         sp(m) = sp(2) * cp(m-1) + cp(2) * sp(m-1)
         cp(m) = cp(2) * cp(m-1) - sp(2) * sp(m-1)
90    continue
c
c     --------------------------------
c*FON Calcul des trois premiers termes
c     --------------------------------
c
      rrmag   = rre * rayt / rgmt
      aor     = 1.d0 / rrmag
      ar      = aor * aor * aor
      p(2,1)  = ct
      dp(2,1) = -st
      p(2,2)  = st
      dp(2,2) = ct
      br      = -(ar + ar) * (g(2,1) * p(2,1) + p(2,2)
     >          * (g(2,2) * cp(2) + g(1,2) * sp(2)))
      bt      = ar * (g(2,1) * dp(2,1) + dp(2,2)
     >          * (g(2,2) * cp(2) + g(1,2) * sp(2)))
      bp      = ar * (g(1,2) * cp(2) - g(2,2) * sp(2)) * p(2,2)
c
c     ------------------------------------------------
c*FON Calcul final des composantes du champ magnetique
c     ------------------------------------------------
c
      do 100 n = 3, kmax
         ar = aor * ar
         do 110 m = 1, n
c
            if (m .eq. n) then
               p(n,n)  = st * p(n-1,n-1)
               dp(n,n) = st * dp(n-1,n-1) + ct * p(n-1,n-1)
            else
               if (const(n,m) .eq. 0) then
                  p(n,m)  = ct * p(n-1,m)
                  dp(n,m) = ct * dp(n-1,m) - st * p(n-1,m)
               else
                  p(n,m)  = ct * p(n-1,m) - const(n,m) * p(n-2,m)
                  dp(n,m) = ct * dp(n-1,m) - st * p(n-1,m)
     >                      - const(n,m) * dp(n-2,m)
               endif
            endif
c
            par = p(n,m) * ar
c
            if (m .eq. 1) then
               temp = g(n,m)
            else
               temp  = g(n,m) * cp(m) + g(m-1,n) * sp(m)
               bp    = bp - (g(n,m) * sp(m) - g(m-1,n) * cp(m))
     >                 * fm(m) * par
            endif
c
            br = br - temp * fn(n) * par
            bt = bt + temp * dp(n,m) * ar
c
110      continue
c
100   continue
c
      bp = bp / st
c
c     ---------------------------------------
c*FON Expression du champ magnetique en gauss
c     ---------------------------------------
c
       b  = sqrt(bp * bp + bt * bt + br * br)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end