subroutine chp05 (year,rre,thetr,phir,br,bt,bp,bb,ier)
c*
c***********************************************************************
c*
c*PRO MAGLIB
c*
c*VER 01.05.30 - V 2.0
c*VER 03.01.06 - V 2.1
c*VER 10.09.27 - V 3.0
c*VER 17.02.22 - 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 des composantes du champ magnetique apres 2000.
c*
c*PAR year  (I) : annee fractionnaire (>= 2005.)
c*
c*PAR rre   (I) : distance radiale geocentrique (rayons terrestres)
c*PAR thetr (I) : colatitude geocentrique (radians)
c*PAR phir  (I) : longitude geocentrique (radians)
c*
c*PAR br    (O) : composante radiale du champ magnetique (gauss)
c*PAR bt    (O) : composante tangentielle du champ magnetique le long
c*PAR           : du meridien positive vers le bas (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 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 - 10.09.27 - Passage aux coefficiens de Schmidt
c*HST                          définitifs et aux variations seculaires 
c*HST                          interpolees
c*HST                          Modif de fn1 (10->100) car précision
c*HST                          des coeff 2005 sur 2 décimales
c*HST                          Modif de fn2 (100->1000) car précision
c*HST                          des coeff 2005 sur 2 décimales
c*HST version 5.0 - 22.02.17 - Passage aux variations séculaires 
c*                             definitives
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,thetr,phir,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(195),lgt(195)
c*LOC lg, lgt : coefficients apres l'an 2005
c
      integer kmax
c*LOC kmax : nombre de pas de calcul (= 14)
c
      integer l,k,jj,nm,ntot,inc
c*LOC l,k,jj,nm,ntot,inc : indices de boucles et de tableaux
c
      double precision gg(14,14),ggt(14,14),hh(14,14),hht(14,14)
      double precision g(14,14),h(14,14),sh(14,14),fn(14),fm(14)
      double precision cp(14),sp(14),p(14,14),dp(14,14),const(14,14)
      double precision f1, f2, tzero, tmold,t,cph,sph,ct,st
      double precision rrmag,ar,aor,temp,par
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 septembre 2010 - coefficients lg
c    Revision fevrier 2017 - coefficients lgt
c
      data lg/
     > -2955463,-166905,507799,-233724,304769,-259450,165776,-51543,
     >  133630,-230583,-19886,124639,26972,67251,
     > -52472,92055,79796,28207,21065,-22523,-37986,14515,10000,
     > -30536,-22700,35441,4272,20895,
     >  18025,-13654,-12345,-16805,-1957,-1355,10385,7360,6956,-2033,
     >  7674,5475,-15134,6363,
     > -1458,-6353,1458,24,-8636,5094,7988,-7446,-6114,-165,-2257,
     >  3873,682,1230,
     >  2535,937,1093,542,-2632,194,-464,2480,762,1120,-1173,-2088,
     > -688,983,
     > -1811,-1971,1017,1622,936,761,-1125,-1276,-487,-6,558,976,
     > -2011,358,
     >  1269,-694,1267,501,-672,-1076,-816,-125,810,876,292,-666,-773,
     > -922,
     >  601,-217,-612,219,142,10,-235,446,-15,476,306,-658,29,-101,
     >  206,-347,377,-86,-21,-231,-209,-793,295,-160,26,-188,144,144,
     > -77,-31,-227,29,90,-79,-58,53,-269,180,-108,16,-158,96,
     > -190,399,-139,-215,-29,-55,21,23,89,238,-38,-263,96,61,
     > -30,40,46,1,-35,2,-36,28,8,-87,-49,-34,-8,88,
     > -16,-88,-76,30,33,28,172,-43,-54,118,-107,-37,-4,75,
     >  63,-26,21,35,53,-5,38,41,-22,-10,-57,-18,-82
     > /
c
       data lgt/
     >  11612,16526,-26746,-11764,-4270,-22808,2082,-12060,710,-4142,
     >  7692,-2858,-3594,-7756,
     > -2462,-1578,2202,882,-8814,2840,4606,3862,-2120,-872,-774,576,
     >	372,-1738,
     >  1752,-902,1078,976,3912,1104,-562,-164,-174,-114,-164,-2114,
     >  1988,-418,
     > -1650,-546,-296,556,1654,892,112,-108,668,-580,274,1302,-56,340,
     > -78,218,-780,-756,-258,596,272,-78,118,-72,-554,170,258,400,
     > -246,460,288,98,298,-130,-560,404,266,340,-16,-62,-86,-26,
     > -236,334,16,-376,-84,-324,148,98,-26,-66,-156,-352,330,-172,
     >  200,46,-24,108,-106,-40,256,50,-2,-64,-122,-128,-124,10,
     >  14,-96,-136,-226,-164,68,-142,-76,20,24,-26,-30,46,42,
     >  22,-40,102,50,-10,0,38,-32,36,-2,-38,-8,-106,-42,
     > -8,-48,-94,6,16,-64,18,8,30,-50,-50,28,-2,-24,
     >  38,38,12,-2,-8,22,-2,-2,26,2,-56,22,24,-2,
     >  14,-2,-22,2,-6,28,-12,-4,-10,-20,-14,12,-6,6,
     > -18,16,-22,6,-8,14,12,2,-6,-32,8,-16,6
     > /
c
      data kmax /14/
c
      data f1 /100.d0/, f2 /1000.d0/
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 tzero /2005.d0/
c
      data tmold /0.d0/
c
      data sh(1,1) /0.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
c     ---------------
c*FON Initialisations
c     ---------------
c
      ier  = 0
c
      if (sh(1,1) .ne. -1.d0) then
         sh(1,1) = -1.d0
c
c     ---------------------------------------------------
c*FON Calcul des coefficients de Shmidt et des constantes
c     ---------------------------------------------------
c
         do 20 k = 1, kmax
            fn(k) = k
            do 10 l = 1, k
               fm(l)      = l - 1
               const(k,l) = dble((k - 2)**2 - (l - 1)**2) /
     >                      dble((2 * k - 3) * (2 * k - 5))
10          continue
20       continue
c
         do 40 k = 2, kmax
            sh(k,1) = dble(2 * k - 3) * sh(k-1,1) / dble(k - 1)
            jj      = 2
            do 30 l = 2,k
               sh(k,l) = sh(k,l-1) * dsqrt(dble((k - l + 1) * jj)
     >                   / (k + l - 2))
               jj      = 1
30          continue
40       continue
c
         ntot = 0
         do 60 k = 2, kmax
c
            gg(k,1)  = dble(lg(ntot+1))
            ggt(k,1) = dble(lgt(ntot+1))
c
            nm = 2 * k - 1
            do 50 l = 2, k
               inc      = 2 * l - 2
               gg(k,l)  = dble(lg(ntot+inc))
               ggt(k,l) = dble(lgt(ntot+inc))
               hh(k,l)  = dble(lg(ntot+inc+1))
               hht(k,l) = dble(lgt(ntot+inc+1))
50          continue
c
            ntot = ntot + nm
60       continue
c
         do 80 k = 2, kmax
            gg(k,1)  = gg(k,1) * sh(k,1) / f1
            ggt(k,1) = ggt(k,1) * sh(k,1) / f2
            do 70 l = 2, k
               gg(k,l)  = gg(k,l)  * sh(k,l) / f1
               ggt(k,l) = ggt(k,l) * sh(k,l) / f2
               hh(k,l)  = hh(k,l)  * sh(k,l) / f1
               hht(k,l) = hht(k,l) * sh(k,l) / f2
70          continue
80       continue
      endif
c
      if (year .ne. tmold) then
         tmold = year
         t     = year - tzero
         do 100 k = 2, kmax
            g(k,1) = gg(k,1) + t * ggt(k,1)
            do 90 l = 2,k
               g(k,l) = gg(k,l) + t * ggt(k,l)
               h(k,l) = hh(k,l) + t * hht(k,l)
90          continue
c
100      continue
c
      endif
c
      cph = dcos(phir)
      sph = dsin(phir)
c
      if (thetr .ne. 0.d0) then
         ct = dcos(thetr)
         st = dsin(thetr)
      else
         ct = 1.0d0
         st = 1.d-15
      endif
      sp(2) = sph
      cp(2) = cph
c
      do 110 k = 3, kmax
c
         sp(k) = sp(2) * cp(k-1) + cp(2) * sp(k-1)
         cp(k) = cp(2) * cp(k-1) - sp(2) * sp(k-1)
c
110   continue
      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
c
      br = -(2.d0 * ar) * (g(2,1) * p(2,1) + p(2,2)
     >     * (g(2,2) * cp(2) + h(2,2) * sp(2)))
      bt = ar * (g(2,1) * dp(2,1) + dp(2,2) * (g(2,2)
     >     * cp(2) + h(2,2) * sp(2)))
      bp = ar * (h(2,2) * cp(2) - g(2,2) * sp(2)) * p(2,2)
c
      do 130 k = 3, kmax
c
         ar = aor * ar
c
         do 120 l = 1, k
c
            if (k .eq. l) then
               p(k,k)  = st * p(k-1,k-1)
               dp(k,k) = st * dp(k-1,k-1) + ct * p(k-1,k-1)
            else
               p(k,l)  = ct * p(k-1,l) - const(k,l) * p(k-2,l)
               dp(k,l) = ct * dp(k-1,l) - st * p(k-1,l) - const(k,l)
     >                * dp(k-2,l)
            endif
c
            par = p(k,l) * ar
c
            if (l .eq. 1) then
               temp = g(k,l)
            else
               temp = g(k,l) * cp(l) + h(k,l) * sp(l)
               bp   = bp - (g(k,l) * sp(l) - h(k,l) *
     >                cp(l)) * fm(l) * par
            endif
            br = br - temp * fn(k) * par
            bt = bt + temp * dp(k,l) * ar
c
120      continue
c
130   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 = dsqrt(br * br + bt * bt + bp * bp)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end