subroutine chp95_15 (year,rre,thetr,phir,br,bt,bp,bb,ier)
c*
c***********************************************************************
c*
c*PRO MAGLIB
c*
c*VER 17.02.23 - V 1.0
c*
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. AKKA
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul du champ magnetique pour la periode
c*ROL         comprise entre 1995 et 2015.
c*
c*PAR year  (I) : annee fractionnaire (entre 1995. et 2015.)
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 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 common    : util, util2
c*
c*INF utilise   : sans objet
c*
c*HST version 1.0 - 17.02.23 - Modele 2010 definitif et ajout
c*HST                          du modele provisoire 2015
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, thetr, phir
      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 cos(-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,5)
c*LOC lg : coefficients entre 1995 et 2015
c
      integer kmax
c*LOC kmax : nombre de pas de calcul (= 14)
c
      integer i,k,l,jj,imin,imax,ntot,nm,inc
c*LOC i,k,l,jj,imin,imax,ntot,nm,inc : indices de boucles et de tableaux
c
      double precision sh(14,14),gg(14,14),hh(14,14)
      double precision g(14,14),h(14,14),fn(14),fm(14)
      double precision cp(14),sp(14),p(14,14),dp(14,14),const(14,14)
      double precision tzero(5)
      double precision tmold, delt, cph, sph, ct, st, rrmag, ar, aor
      double precision temp, par, fmin, fmax
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, nouveaux lg
c    Champ 1995  révisé octobre 2006
c
      data (lg(i,1),i=1,195)/
     > -296920, -17840, 53060, -22000, 30700, -23660, 16810, -4130,
     > 13350,-22670,
     > -2620, 12490, 3020, 7590, -4270, 9400, 7800, 2620, 2900, -2360,
     > -4180, 970, 1220, -3060, -2140, 3520, 460, 2350, 1650, -1180,
     > -1430, -1660, -550, -170, 1070, 680, 670, -170, 680, 720,
     > -1700, 670, -10, -580, 190, 10, -930, 360, 770, -720,
     > -690, 10, -250, 280, 40, 50, 240, 40, 170, 80,
     > -240, -20, -60, 250, 60, 110, -60, -210, -90, 80, 
     > -140, -230, 90, 150, 60, 110, -50, -160, -70, -40,
     >  40, 90, -200, 30, 150, -100, 120, 80, -60, -80,
     > -80, -10, 80, 100, 50, -20, -80, -80, 30, -30,
     > -60, 10, 20, 0, -40, 40, -10, 50, 40, -50,
     > 20, -10, 20, -20, 50, 10, 10, -20, 0, -70,75*0/
c
c    Champ 2000  révisé octobre 2006
c
      data (lg(i,2),i=1,195)/
     > -296194,-17282,51861,-22677,30684,-24816,16709,-4580,13396,
     > -22880,
     > -2276,12521,2934,7145,-4911,9323,7868,2726,2500,-2319,
     > -4030,1198,1113,-3038,-2188,3514,438,2223,1719,-1304,
     > -1331,-1686,-393,-129,1063,723,682,-174,742,637,
     > -1609,651,-59,-612,169,7,-904,438,790,-740,
     > -646,0,-242,333,62,91,240,69,148,73,
     > -254,-12,-58,244,66,119,-92,-215,-79,85,
     > -166,-215,91,155,70,89,-79,-149,-70,-21,
     > 50,94,-197,30,134,-84,125,63,-62,-89,
     > -84,-15,84,93,38,-43,-82,-82,48,-26,
     > -60,17,17,0,-31,40,-5,49,37,-59,
     > 10,-12,20,-29,42,2,3,-22,-11,-74,
     > 27,-17,1,-19,13,15,-9,-1,-26,1,
     > 9,-7,-7,7,-28,17,-9,1,-12,12,
     > -19,40,-9,-22,-3,-4,2,3,9,25,
     > -2,-26,9,7,-5,3,3,0,-3,0,
     > -4,3,-1,-9,-2,-4,-4,8,-2,-9,
     > -9,3,2,1,18,-4,-4,13,-10,-4,
     > -1,7,7,-4,3,3,6,-1,3,4,
     > -2,0,-5,1,-9/
c
c    Champ 2005 coeff. definitifs révisés 2010.10.05
c
      data (lg(i,3),i=1,195)/
     > -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
c    Champ 2010 coeff. definitifs révisés 2017.02.23
c
      data (lg(i,4),i=1,195)/
     > -2949657,-158642,494426,-239606,302634,-270854,166817,-57573,
     > 133985,-232654,-16040,123210,25175,
     >  63373,-53703,91266,80897,28648,16658,-21103,-35683,16446,8940,
     > -30972,-23087,35729,4458,
     >  20026,18901,-14105,-11806,-16317,-1,-803,10104,7278,6869,
     > -2090,7592,4418,-14140,
     >  6154,-2283,-6626,1310,302,-7809,5540,8044,-7500,-5780,-455,
     > -2120,4524,654,
     >  1400,2496,1046,703,164,-2761,492,-328,2441,821,1084,-1450,
     > -2003,-559,
     >  1183,-1934,-1741,1161,1671,1085,696,-1405,-1074,-354,164,550,
     >  945,-2054,
     >  345,1151,-527,1275,313,-714,-1238,-742,-76,797,843,214,-842,
     > -608,
     > -1008,701,-194,-624,273,89,-10,-107,471,-16,444,245,-722,-33,
     > -96,213,-395,309,-199,-103,-197,-280,-831,305,-148,13,-203,167,
     >  165,-66,-51,-176,54,85,-79,-39,37,-251,179,-127,12,-211,
     >  75,-194,375,-186,-212,-21,-87,30,27,104,213,-63,-249,95,
     >  49,-11,59,52,0,-39,13,-37,27,21,-86,-77,-23,4,
     >  87,-9,-89,-87,31,30,42,166,-45,-59,108,-114,-31,-7,
     >  78,54,-18,10,38,49,2,44,42,-25,-26,-53,-26,-79
     > /
c
c    Champ 2015 coeff. provisoires révisés 2017.02.23
c
      data (lg(i,5),i=1,195)/
     > -2944200,-150100,479710,-244510,301290,-284560,167670,-64190,
     >  135070,-235230,-11530,122560,24490,
     >  58200,-53840,90760,81370,28330,12040,-18870,-33490,18090,7040,
     > -32950,-23260,36010,4730,
     >  19240,19700,-14090,-11930,-15750,1600,410,10020,7000,6770,
     > -2080,7270,3320,-12990,
     >  5890,-2890,-6670,1320,730,-7090,6260,8160,-7610,-5410,-680,
     >  -1950,5180,570,
     >  1500,2440,940,340,-280,-2740,680,-220,2420,880,1010,-1690,
     > -1830,-320,
     >  1330,-2060,-1460,1340,1620,1170,570,-1590,-910,-200,210,540,
     >  880,-2160,
     >  310,1080,-330,1180,70,-680,-1330,-690,-10,780,870,100,-910,
     > -400,
     > -1050,840,-190,-630,320,10,-40,50,460,-50,440,180,-790,-70,
     > -60,210,-420,240,-280,-180,-120,-360,-870,310,-150,-10,-230,200,
     >  200,-70,-80,-110,60,80,-70,-20,20,-220,170,-140,-20,-250,
     >  40,-200,350,-240,-190,-20,-110,40,40,120,190,-80,-220,90,
     >  30,10,70,50,-10,-30,30,-40,20,20,-90,-90,-10,00,
     >  70,00,-90,-90,40,40,50,160,-50,-50,100,-120,-20,-10,
     >  80,40,-10,-10,30,40,10,50,50,-30,-40,-40,-30,-80
     > /
c
      data kmax /14/
c
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 /1995.d0,2000.d0,2005.d0,2010.d0,2015.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
      endif
c
c     ------------------------------------------------------------
c*FON Calcul des deux modeles de champ entourant la date de calcul
c*FON et calcul de l'intervalle de temps pour l'interpolation
c     ------------------------------------------------------------
c
      if (year .ne. tmold) then
         tmold = year
         imin  = idint((year - 1995.d0) / 5.d0) + 1
         imax  = imin + 1

c        ***********************************************
c        Facteurs de normalisation des modèle min et max
c        ***********************************************
         if (tzero(imin) .ge. 2005.d0) then
            fmin = 100.d0
         else
            fmin = 10.d0
         endif
         if (tzero(imax) .ge. 2005.d0) then
            fmax = 100.d0
         else
            fmax = 10.d0
         endif


         delt  = (year - tzero(imin)) / 5.d0
c
         ntot = 0
         do 60  k = 2, kmax
            gg(k,1) = dble(lg(ntot+1,imin))/fmin
     >                + ((dble(lg(ntot+1,imax))/fmax)
     >                  - (dble(lg(ntot+1,imin))/fmin))
     >                * delt
            l = 1
c
            nm = 2 * k - 1
            do 50 l = 2, k
               inc = 2 * l - 2
               gg(k,l) =  dble(lg(ntot+inc,imin))/fmin
     >                    + ((dble(lg(ntot+inc,imax))/fmax)
     >                      - (dble(lg(ntot+inc,imin))/fmin))
     >                     * delt

               hh(k,l) = dble(lg(ntot+inc+1,imin))/fmin
     >                   + ((dble(lg(ntot+inc+1,imax))/fmax)
     >                   - (dble(lg(ntot+inc+1,imin))/fmin))
     >                   * delt
c
50          continue
            ntot = ntot + nm
60       continue
         do 80 k = 2, kmax
            g(k,1) = gg(k,1) * sh(k,1)
            do 70 l = 2, k
                g(k,l) = gg(k,l) * sh(k,l)
               h(k,l) = hh(k,l) * sh(k,l)
70         continue
80      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 90 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
90    continue
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
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
      k = 2
      l = 12
c
      do 110 k = 3, kmax
c
         ar = aor * ar
c
         do 100 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
100      continue
c
110   continue
c
      bp = bp / st
c
      bb = dsqrt(br * br + bt * bt + bp * bp)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end