subroutine dgrf00 (year,rre,thet,phi,br,bt,bp,bb,ier)
c*
c***********************************************************************
c*
c*          "Copyright [c] CNES 98 - tous droits reserves"
c*          **********************************************
c*
c*PRO MAGLIB
c*
c*VER 01.05.30 - V 2.0
c*VER 03.01.06 - V 2.1
c*VER 01.01.07 - V 3.0
c*VER 10.09.27 - V 4.0
c*
c*AUT spec. CNES - JC KOSIK - fevrier 2001
c*AUT port. CS-SI
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul du champ magnetique d'origine interne avec les
c*ROL         coefficients du modele  DGRF 2000.
c*
c*PAR year (I) : annee fractionnaire (>= 2000.)
c*
c*PAR rre  (I) : distance radiale geocentrique (rayons terrestres)
c*PAR thet (I) : colatitude geocentrique (radians)
c*PAR phi  (I) : longitude geocentrique (radians)
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 bb   (O) : module du champ magnetique (gauss)
c*
c*PAR ier  (O) : code de retour
c*
c*NOT ier      : sans objet
c*
c*NOT Dans le cas ou l'angle thet est nul, son sinus est
c*NOT remplacee par 1.d-15 afin d'eviter des divisions par 0 et
c*NOT son cosinus est remplace par 1.0d0.
c*NOT Dans le cas ou l'angle thet est egal a pi, son sinus est
c*NOT remplacee par 1.d-15 afin d'eviter des divisions par 0 et
c*NOT son cosinus est remplace par -1.0d0.
c*
c*NOT common   : util, util2
c*
c*INF utilise  : sans objet
c*
c*HST version 2.0 - 01.05.30 - Enrichissement de la maglib au CDPP
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - adoption des coefficients definitifs
c*HST                          du champ 2000
c*HST version 4.0 - 10.09.27 - Passage aux variations seculaires 
c*HST                          definitives
c*HST                          Modif de fn2 (100->1000) car précision
c*HST                          des coeff 2005 sur 2 décimales
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, rre, thet, phi
      double precision br, bt, bp, bb
      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(14,14),lgt(14,14)
c*LOC lg, lgt : coefficients du modele dgrf 2000
c
      integer kmax
c*LOC kmax : nombre de pas de calcul (= 14)
c
      integer n,m,jj
c*LOC n,m,jj : indices de boucles et de tableaux
c
      double precision gg(14,14),ggt(14,14)
      double precision shmit(14,14)
      double precision g(14,14)
      double precision const(14,14),fn(14),fm(14)
      double precision p(14,14),dp(14,14),sp(14),cp(14)
      double precision tmold,tzero,sph,st,t,temp,par,f2
      double precision cph,ct,f1,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     Revision octobre 2006 - coefficients lg 
c     Revision septembre 2010 - coefficients lgt
c
      data lg/10,
     >  -296194,-22677,13396,9323,-2188,723,790,244,50,-26,27,-22,-2,
     >  51861,-17282,30684,-22880,7868,3514,682,-740,66,94,-60,-17,-3,
     > -9,
     > -24816,-4580,16709,12521,2500,2223,742,0,-92,30,17,-19,2,3,
     > -2276,2934,-4911,7145,-4030,-1304,-1609,333,-79,-84,-31,15,9,1,
     >  2726,-2319,1198,-3038,1113,-1686,-59,91,-166,63,-5,-1,-2,-4,
     >  438,1719,-1331,-393,1063,-129,169,69,91,-89,37,1,9,13,
     > -174,637,651,-612,7,438,-904,73,70,-15,10,-7,-5,-4,
     > -646,-242,62,240,148,-254,-58,-12,-79,93,20,7,3,7,
     >  119,-215,85,-215,155,89,-149,-21,-70,-43,42,17,-3,-4,
     > -197,134,125,-62,-84,84,38,-82,48,-82,3,1,-4,3,
     >  17,0,40,49,-59,-12,-29,2,-22,-74,-11,12,-1,-1,
     >  1,13,-9,-26,9,-7,-28,-9,-12,-19,-9,40,-2,4,
     > -4,3,25,-26,7,3,0,0,3,-9,-4,8,-4,0,
     > -9,2,18,-4,-10,-1,7,3,6,3,-2,-5,-9,1
     > /
c      
       data lgt/1000,
     >  12954,-13908,-660,-2350,-1640,260,176,80,116,86,50,10,8,
     > -21622,11830,-4142,-3566,2232,602,272,-92,204,72,-24,20,2,4,
     > -22580,-11486,-2628,-1142,-7870,-2670,508,-330,-506,116,-56,4,2,
     >  0,
     >  5748,-4736,-6724,-8398,4628,-1228,1912,1086,204,292,150,-12,-2,
     >  36,
     >  1894,1334,5070,-312,-2260,110,-1736,640,-302,-258,70,-42,-36,-6,
     > -216,1670,1930,3946,-490,-130,-464,494,214,-372,-128,38,12,-24,
     > -586,-1790,-294,-466,-92,1428,808,-376,472,50,-142,-18,40,6,
     >  692,326,124,270,-774,-184,232,628,-670,-108,12,-34,32,10,
     > -140,124,266,358,144,-258,428,408,426,-472,-86,20,-10,28,
     > -82,-142,34,-104,48,-60,-176,94,242,-204,-102,12,8,10,
     >  98,20,92,-28,-136,38,-114,-212,-22,-106,-198,-48,36,10,
     >  32,28,26,66,0,24,22,-36,-76,0,-98,-2,-58,2,
     > -30,-14,-24,-6,-18,20,2,4,-4,6,12,16,64,-20,
     >  28,26,-16,-28,-14,12,-14,-18,-14,16,-4,-14,16,-56
     > / 
c
      data f1 /10.d0/, f2 /1000.d0/
c
      data kmax /14/
c
      data p(1,1)  /1.d0/
      data cp(1)   /1.d0/
      data dp(1,1) /0.d0/
      data sp(1)   /0.d0/
c
      data shmit(1,1) /0.d0/
c
      data tmold /0.d0/
c
      data tzero /2000.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
c     ---------------------------------------------------------------
c*FON Initialisation a 0 des composantes du champ ainsi que du module
c     ---------------------------------------------------------------
c
      br = 0.0d0
      bt = 0.0d0
      bp = 0.0d0
      bb = 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, 14
            fn(n) = n
            do 20 m = 1, 14
               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, 14
            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 Calcul des tableaux gg(n,m) et ggt(n,m)
c     ---------------------------------------
c
         f1 = dble(lg(1,1))
         f2 = dble(lgt(1,1))
         do 50 n = 1, kmax
            do 60 m = 1, kmax
               gg(n,m)  = dble(lg(n,m)) * shmit(n,m) / f1
               ggt(n,m) = dble(lgt(n,m)) * shmit(n,m) / f2
60          continue
50       continue
c
c     --------------------------------------------------
c*FON Traitements effectues si ce n'est pas le 1er appel
c     --------------------------------------------------
c
      endif
c
c     ------------------------
c*FON Calcul du tableau g(n,m)
c     ------------------------
c
      if (year .ne. tmold) then
c
         tmold = year
         t     = year - tzero
c
         do 70 n = 1, kmax
            do 80 m= 1, kmax
               g(n,m) = gg(n,m) + t * ggt(n,m)
80          continue
70       continue
c
      endif
c
c     ----------------------------------------------------
c*FON Calcul des cosinus directeurs de la colatitude et de
c*FON la longitude geocentriques
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 Premiers termes des composantes du champ magnetique
c     ---------------------------------------------------
c
      rrmag   = rre * rayt / rgmt
      aor     = 1.d0 / rrmag
      ar      = aor * aor * aor
c
      p(2,1)  = ct
      dp(2,1) = -st
      p(2,2)  = st
      dp(2,2) = ct
c
      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 Evaluation finale 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
      br = br / 100000.d0
      bt = bt / 100000.d0
      bp = bp / 100000.d0
c
      bb = sqrt(bp * bp + bt * bt + br * br)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end