subroutine grad95 (year,rre,thet,phi,br,bt,bp,bb,gradb,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*
c*AUT spec. CNES - JC KOSIK - mars 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 DGRF95.
c*
c*PAR year   (I) : annee fractionnaire (>= 1995.)
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 gradb  (O) : gradients du champ magnetique
c*PAR            : les quantites sont evaluees en unite de rayon terrestre
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 Attention : coefficient g7,4 et sa derivee non imprimes
c*NOT dans bulletin iaga. coefficients provisoires indiques
c*
c*NOT common     : util, util2
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*
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
      double precision gradb(3,3)
      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
      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 dgrf95
c
      integer n,m,jj
c*LOC n,m,jj : indices de boucles et de tableaux
c
      integer kmax
c*LOC kmax : nombre de pas de calcul (= 14)
c
      double precision dbrdr,dbrdt,dbrdp,dbtdr,dbtdt,dbtdp
      double precision dbpdr,dbpdt,dbpdp
c*LOC dbrdr,dbrdt,dbrdp : derivees partielles des composantes radiale, 
c*LOC dbtdr,dbtdt,dbtdp : tangentielle et azimuthale du champ total
c*LOC dbpdr,dbpdt,dbpdp :
c
      double precision gg(14,14),ggt(14,14),shmit(14,14),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 d2p(14,14)
      double precision tmold,tzero,sphi,st,t,f2
      double precision cphi,ct,f1,aor,aor3,aorn,daorn,rmag,rkm
      double precision temp1,dtemp1,d2temp1
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
      data lg/
     s 10,-296920,-22000,13350,9400,-2140,680,770,250,40,-30,3*0,
     s 53060,-17840,30700,-22670,7800,3520,670,-720,60,90,-60,3*0,
     s -23660,-4130,16810,12490,2900,2350,680,10,-60,30,20,3*0,
     s -2620,3020,-4270,7590,-4180,-1180,-1700,280,-90,-100,-40,3*0,
     s 2620,-2360,970,-3060,1220,-1660,-10,50,-140,80,-10,3*0,
     s 460,1650,-1430,-550,1070,-170,190,40,90,-80,40,3*0,
     s -170,720,670,-580,10,360,-930,80,60,-10,20,3*0,
     s -690,-250,40,240,170,-240,-60,-20,-50,100,20,3*0,
     s 110,-210,80,-230,150,110,-160,-40,-70,-20,50,3*0,
     s -200,150,120,-60,-80,80,50,-80,30,-80,10,3*0,
     s  10,0,40,50,-50,-10,-20,10,-20,-70,0,45*0/
c
      data lgt/100,
     > 1452, -1354, 92, -154, -96, 86, 40, -12, 20, 8,
     > 54, -44, -4, -2398, 1116, -32, -420, 136, -12, 24,
     > -40, 12, 8, 0, -34, -6, -18, -2312, -900, -202,
     > 62, -800, -254, 124, -20, -64, 0, -6, -38, 4,
     > 6, 688, -172, -1282, -890, 300, -248, 182, 106, 22,
     > 32, 18, 30, 18, 2, 212, 82, 456, 44, -214,
     > -52, -98, 82, -52, -34, 10, -2, -4, -8, -44,
     > 138, 198, 314, -14, 82, -42, 58, 2, -18, -6,
     > 2, 18, 26, -8, -166, -38, -64, -6, 156, 52,
     > -14, 20, -10, -20, -14, -10, -8, 88, 16, 44,
     > 0, -44, -28, 4, 16, -58, -14, 0, 14, 6,
     > 14, 18, -10, 10, 30, 10, -42, 22, 38, 0,
     > -46, -16, 34, -6, -8, 6, -32, 10, -4, -8,
     > 8, -24, -4, 36, -4, -14, 2, -8, 6, 14,
     > 0, 0, -2, -18, -4, -18, -16, -4, -8, -22,
     > 24, -2, -2, 2, 26, -18, -52, 18, -14, -56,
     > -18, -24, -38, -18, 80, -4, 8, -8, 6, 50,
     > -52, 14, 6, 0, 0, 6, -18, -8, 16, -8,
     > 0, -18, 4, 36, -8, -20, -2, 14, 6, 12,
     > 6, -4, -10, -18, 2/
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/1995.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
      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     mise a zero des constantes avant la premiere initialisation
c
         do 8 n = 1, 14
            fn(n) = 0.d0
            do 9 m = 1, 14
               fm(m)      = 0.d0
               const(n,m) = 0.d0
               shmit(n,m) = 0.d0
               dp(n,m)    = 0.d0
               d2p(n,m)   = 0.d0
9           continue
8        continue

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
      cphi = cos(phi)
      sphi = 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) = sphi
      cp(2) = cphi
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
      rkm = rre * rayt
      rmag = rkm / rgmt
c
      aor   = rgmt / rkm
      aor3  = aor * aor * aor
c
      p(2,1)   =  ct
      dp(2,1)  = -st
      d2p(2,1) = -ct
      p(2,2)   =  st
      dp(2,2)  =  ct
      d2p(2,2) = -st
c
      br = -2.d0 * aor3 *
     >   (g(2,1) * ct + (g(2,2) * cphi + g(1,2) * sphi) * st)
c
      bt = aor3 *
     >   (-g(2,1) * st + (g(2,2) * cphi + g(1,2) * sphi) * ct)
c
      bp = aor3 * (g(1,2) * cphi - g(2,2) * sphi)
c
c     ------------------------------------------
c*FON Calcul des 9 derivees partielles de depart
c     ------------------------------------------
c
      dbrdr = 6.d0 * aor3 *
     >   (g(2,1) * ct + (g(2,2) * cphi + g(1,2) * sphi) * st) / rre
c
      dbrdt = -2.d0 * aor3 *
     >   (-g(2,1) * st + (g(2,2) * cphi + g(1,2) * sphi) * ct)
c
      dbrdp = -2.d0 * aor3 *
     >   (g(1,2) * cphi - g(2,2) * sphi) * st
c
      dbtdr = -3.d0 * aor3 *
     >   (-g(2,1) * st + (g(2,2) * cphi + g(1,2) * sphi) * ct) / rre
c
      dbtdt = aor3 *
     >   (-g(2,1) * ct - (g(2,2) * cphi + g(1,2) * sphi) * st)
c
      dbtdp = aor3 *
     >   (-g(2,2) * sphi + g(1,2) * cphi) * ct
c
      dbpdr = -3.d0 * aor3 *
     >   (g(1,2) * cphi - g(2,2) * sphi) / rre
c
      dbpdt = 0.d0
c
      dbpdp = aor3 * (-g(2,2) * cphi - g(1,2) * sphi)
c
      aorn  = aor3
c
c     -----------------------------------------------------
c*FON Evaluation finale des composantes du champ magnetique
c     -----------------------------------------------------
c
      do 100 n = 3, kmax
c
         aorn  = aor * aorn
         daorn = -dble(n + 1) * aorn / rre
c
         do 110 m = 1, n
c
            if (m .eq. n) then
c
               p(n,n)   = st * p(n-1,n-1)
               dp(n,n)  = st * dp(n-1,n-1) + ct * p(n-1,n-1)
               d2p(n,n) = ct * dp(n-1,n-1) + st * d2p(n-1,n-1)
     >                    - st * p(n-1,n-1) + ct * dp(n-1,n-1)
c
            else
c
               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)
               d2p(n,m) = - st * dp(n-1,m) + ct * d2p(n-1,m)
     >                    - ct * p(n-1,m) - st * dp(n-1,m)
     >                    - const(n,m) * d2p(n-2,m)
c
            endif
c
            if (m .eq. 1) then
c
               temp1   = g(n,m)
               dtemp1  = 0.d0
               d2temp1 = 0.d0
c
            else
c
               temp1   =  g(n,m) * cp(m) + g(m-1,n) * sp(m)
               dtemp1  =  fm(m) * (-g(n,m) * sp(m) + g(m-1,n) * cp(m))
               d2temp1 = -fm(m) * fm(m) * temp1
c
            endif
c
            br = br - temp1 * p(n,m) * aorn * fn(n)
            bt = bt + temp1 * dp(n,m) * aorn
            bp = bp + dtemp1 * p(n,m) * aorn / st
c
            dbrdr =  dbrdr - temp1 * p(n,m) * daorn * fn(n)
            dbrdt =  dbrdt - temp1 * dp(n,m) * aorn * fn(n)
            dbrdp =  dbrdp - dtemp1 * p(n,m) * aorn * fn(n)
c
            dbtdr =  dbtdr + temp1 * dp(n,m) * daorn
            dbtdt =  dbtdt + temp1 * d2p(n,m) * aorn
            dbtdp =  dbtdp + dtemp1 * dp(n,m) * aorn
c
            dbpdr =  dbpdr + dtemp1 * p(n,m) * daorn / st
            dbpdt =  dbpdt + dtemp1 * aorn *
     >               (dp(n,m) - ct * p(n,m) / st) / st
            dbpdp =  dbpdp + d2temp1 * p(n,m) * aorn / st
c
110      continue
c
100   continue
c
c      mise en memoire des gradients de champ
c
      gradb(1,1) = dbrdr
      gradb(1,2) = dbrdt
      gradb(1,3) = dbrdp
      gradb(2,1) = dbtdr
      gradb(2,2) = dbtdt
      gradb(2,3) = dbtdp
      gradb(3,1) = dbpdr
      gradb(3,2) = dbpdt
      gradb(3,3) = dbpdp
c
      bb = sqrt(bp * bp + bt * bt + br * br)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end