subroutine tsyg87 (indval,ps,x,y,z,bx,by,bz,ier)
c*
c***********************************************************************
c*
c*          "Copyright [c] CNES 98 - tous droits reserves"
c*          **********************************************
c*
c*PRO MAGLIB
c*
c*VER 99.03.31 - V 1.0
c*VER 01.05.30 - V 2.0
c*VER 03.01.06 - V 2.1
c*
c*AUT spec. Nikolai a. TSYGANENKO
c*AUT       Institute of Physics, Leningrad University
c*AUT       Stary Petergof 1989 04 Leningrad USSR
c*AUT port. CISI
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul du champ externe d'origine magnetospherique en
c*ROL         fonction de l'angle de tilt (ps) et de l'indice
c*ROL         geomagnetique (indval).
c*
c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ
c*
c*PAR ps     (I) : angle de tilt (radians)
c*
c*PAR x      (I) : coordonnee solaire magnetique en x (rayons terrestres)
c*PAR y      (I) : coordonnee solaire magnetique en y (rayons terrestres)
c*PAR z      (I) : coordonnee solaire magnetique en z (rayons terrestres)
c*
c*PAR bx     (O) : composante solaire magnetospherique x du champ 
c*PAR            : externe (nanoteslas)
c*PAR by     (O) : composante solaire magnetospherique y du champ 
c*PAR            : externee (nanoteslas)
c*PAR bz     (O) : composant solaire magnetospherique z du champ 
c*PAR            : externe (nanoteslas)
c*
c*PAR ier    (O) : code de retour
c*
c*NOT indval     : 1 ---> Kp =  0  , 0+
c*NOT indval     : 2 ---> Kp =  1- , 1 , 1+
c*NOT indval     : 3 ---> Kp =  2- , 2 , 2+
c*NOT indval     : 4 ---> Kp =  3- , 3 , 3+
c*NOT indval     : 5 ---> Kp =  4- , 4 , 4+
c*NOT indval     : 6 ---> Kp >  5-
c*
c*NOT ier      : sans objet
c*
c*NOT common   : util
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.05.30 - correction de commentaires de code
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
      integer indval
      double precision ps
      double precision x, y, z
      double precision bx, by, bz
      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
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer ip
c*LOC ip : valeur initiale de indval (= 100)
c
      integer j
c*LOC j : indice de boucles
c
      double precision tga(6,32),pa(32)
c
      double precision c1,rrc2,xn,rh,dy,delx,b0,b1,xn2,xn22,xn21,psi
      double precision sps,cps,zs,rps,zp,rt,fy,s2m,g2m,fc,xnx2,xnx,xc1
      double precision x1,xc2,x2,xc22,xr2,xc12,b,b20,bp,b2p,bm,b2m
      double precision xa1,xap1,xam1,xa2,xap2,xam2,xna,xnap,xnam,f,fp
      double precision fm,s0,s0p,dstr,b2,xnr,adln,zm,xln1,xlnp1,xlnm1
      double precision s1p,s2p,g1p,g2,g2p,ex1,ex2,y2,z2,s0m,xln2,xlnp2
      double precision yz,zsm,rr2,xsm,zn,bzsm,by1,brsm,bxsm,xlnm2
      double precision aln,s1,s2,s1m,g1,g1m
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 (tga(1,j),j = 1,32)/
     >  -0.09673d0,    -10.63d0,    1.210d0,    34.57d0,
     >  -0.04502d0,    -0.06553, -0.02952d0,   0.3852d0,
     >  -0.03665d0,    -2.084d0, 0.001795d0,  0.00638d0,
     >    -23.49d0,   0.06082d0,  0.01642d0, -0.02137d0,
     >     32.21d0,  -0.04373d0, -0.02311d0,  -0.2832d0,
     > -0.002303d0, -0.000631d0,   -6.397d0,    -967.d0,
     >    -8650.d0,    -20.55d0,     5.18d0,   -2.796d0,
     >     2.715d0,     13.58d0,    8.038d0,    29.21d0/
c
      data (tga(2,j),j = 1,32)/
     >   -0.4850d0,    -12.84d0,     1.856d0,    40.06d0,
     >   -0.0294d0,  -0.09071d0,  -0.02993d0,   0.5465d0,
     >  -0.04928d0,    -2.453d0,  0.001587d0, 0.007402d0,
     >    -29.41d0,   0.08101d0,   0.02322d0,  -0.1091d0,
     >     40.75d0,  -0.07995d0,  -0.03859d0,  -0.2755d0,
     > -0.002759d0, -0.000408d0,    -6.189d0,   -957.8d0,
     >    -7246.d0,    -25.51d0,     5.207d0,   -4.184d0,
     >     2.641d0,     16.56d0,     7.795d0,    29.36d0/
c
      data (tga(3,j),j = 1,32)/
     >     -1.132d0,    -18.05d0,     2.625d0,    48.55d0,
     >  -0.004868d0,   -0.1087d0,  -0.03824d0,   0.8514d0,
     >    -0.0522d0,    -2.881d0, -0.000295d0, 0.009055d0,
     >     -29.48d0,   0.06394d0,   0.03864d0,  -0.2288d0,
     >      41.77d0,  -0.05849d0,  -0.06443d0,  -0.4683d0,
     >   0.001222d0, -0.000519d0,    -3.696d0,   -991.1d0,
     >     -6955.d0,    -31.43d0,     4.878d0,   -3.151d0,
     >      3.277d0,     19.19d0,     7.248d0,    28.99d0/
c
      data (tga(4,j),j = 1,32)/
     >    -1.003d0,    -16.98d0,    3.140d0,   52.81d0,
     >  -0.08625d0,   -0.1478d0, -0.03501d0,  0.5500d0,
     >  -0.07778d0,    -2.970d0, 0.002086d0, 0.01275d0,
     >    -26.79d0,   0.06328d0,  0.03622d0, 0.08345d0,
     >     39.72d0,  -0.06009d0, -0.07825d0, -0.9698d0,
     >  0.000178d0, -0.000573d0,  -0.9328d0,  -872.5d0,
     >    -5851.d0,    -39.68d0,    4.902d0,  -3.848d0,
     >     2.790d0,     20.91d0,    6.193d0,   26.81d0/
c
      data (tga(5,j),j = 1,32)/
     >     -1.539d0,    -14.29d0,    3.479d0,   53.36d0,
     >  -0.004201d0,   -0.2043d0, -0.03932d0,   0.6409d0,
     >    -0.1058d0,    -3.221d0, -0.00114d0,  0.02166d0,
     >     -30.43d0,   0.04049d0,  0.05464d0, 0.008884d0,
     >      42.00d0,  -0.01035d0,  -0.1053d0,   -1.630d0,
     >   0.003802d0, -0.001029d0,    4.204d0,   -665.6d0,
     >    -1011.0d0,    -43.49d0,    4.514d0,   -2.948d0,
     >       2.99d0,     21.59d0,    6.005d0,    22.00d0/
c
      data (tga(6,j),j = 1,32)/
     >   -2.581d0,  -7.726d0,    5.045d0,   53.31d0,
     >  0.02262d0, -0.1972d0, -0.01981d0,  0.4280d0,
     >  -0.1055d0,  -5.075d0, 0.002762d0, 0.03277d0,
     >   -27.35d0, 0.04986d0,  0.06119d0, -0.1211d0,
     >    47.48d0, -0.0502d0,  -0.1477d0,   0.838d0,
     > -0.01008d0, -0.0057d0,    9.231d0,  -674.3d0,
     >   -900.0d0,  -74.43d0,    4.658d0,  -3.245d0,
     >     3.39d0,   21.80d0,    5.620d0,   25.17d0/
c
      data ip /100/
c
      data fc /0.3183099031d0/
c
      data rt /30.d0/
c
      data x1 /4.d0/
      data x2 /5.d0/
c
      data psi /10.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
      if (indval .eq. ip) go to 2
c
      ip = indval
c
      do 10 j = 1, 32
         pa(j) = tga(ip,j)
10    continue
c
      c1   = pa(29)**2
      rrc2 = pa(27)**2
      dstr = pa(26) / rrc2 * 4.d0
      xn   = pa(28)
      rh   = pa(31)
      dy   = pa(30)
      delx = pa(32)
      b0   = pa(23)
      b1   = pa(24)
      b2   = pa(25)
      xn21 = (xn-x1)**2
      xn2  = xn - x2
      xnr  = 1.d0 / xn2
      xn22 = xn2**2
      adln = log(xn22 / xn21)
c
2     continue
c
      if (abs(ps - psi) .lt. 1.d-6) go to 3
c
      psi = ps
      sps = sin(ps)
      cps = cos(ps)
      rps = rh * sps
c
c     ----------------------------------------------------------
c*FON Le traitement commence ici si les parametress indval et ps 
c*FON restent inchanges apres l'appel precedant
c     ----------------------------------------------------------
c
3     continue
c
      zs    = z - rps
      zp    = z - rt
      zm    = z + rt
      fy    = fc / (1.d0 + (y / dy)**2)
      xnx   = xn - x
      xnx2  = xnx**2
      xc1   = x - x1
      xc2   = x - x2
      xc22  = xc2**2
      xr2   = xc2 * xnr
      xc12  = xc1**2
      b20   = zs**2 + c1
      b2p   = zp**2 + c1
      b2m   = zm**2 + c1
      b     = sqrt(b20)
      bp    = sqrt(b2p)
      bm    = sqrt(b2m)
      xa1   = xc12 + b20
      xap1  = xc12 + b2p
      xam1  = xc12 + b2m
      xa2   = 1.d0 / (xc22 + b20)
      xap2  = 1.d0 / (xc22 + b2p)
      xam2  = 1.d0 / (xc22 + b2m)
      xna   = xnx2 + b20
      xnap  = xnx2 + b2p
      xnam  = xnx2 + b2m
      f     = b20 - xc22
      fp    = b2p - xc22
      fm    = b2m - xc22
      xln1  = log(xn21 / xna)
      xlnp1 = log(xn21 / xnap)
      xlnm1 = log(xn21 / xnam)
      xln2  = xln1 + adln
      xlnp2 = xlnp1 + adln
      xlnm2 = xlnm1 + adln
      aln   = 0.25d0 * (xlnp1 + xlnm1 -2.d0 * xln1)
      s0    = (atan(xnx / b) + pid)/ b
      s0p   = (atan(xnx / bp) + pid) / bp
      s0m   = (atan(xnx / bm) + pid) / bm
      s1    = (xln1 * .5d0 + xc1 * s0) / xa1
      s1p   = (xlnp1 * .5d0 + xc1 * s0p) / xap1
      s1m   = (xlnm1 * .5d0 + xc1 * s0m) / xam1
      s2    = (xc2 * xa2 * xln2 - xnr - f * xa2 * s0) * xa2
      s2p   = (xc2 * xap2 * xlnp2 - xnr - fp * xap2 * s0p) * xap2
      s2m   = (xc2 * xam2 * xlnm2 - xnr - fm * xam2 * s0m) * xam2
      g1    = (b20 * s0 - 0.5d0 * xc1 * xln1) / xa1
      g1p   = (b2p * s0p - 0.5d0 * xc1 * xlnp1) / xap1
      g1m   = (b2m * s0m - 0.5d0 * xc1 * xlnm1) / xam1
      g2    = ((0.5d0 * f * xln2+ 2.d0 * s0 * b20 * xc2) 
     >        * xa2 + xr2) * xa2
      g2p   = ((0.5d0 * fp * xlnp2 + 2.d0 * s0p * b2p * xc2)
     >        * xap2 + xr2) * xap2
      g2m   = ((0.5d0 * fm * xlnm2 + 2.d0 * s0m *b2m * xc2)
     >        * xam2 + xr2) * xam2
      bx    = fy * (b0 * (zs * s0 - 0.5d0 * (zp * s0p + zm * s0m))
     >        + b1 * (zs * s1 - 0.5d0 * (zp * s1p + zm * s1m))
     >        + b2 * (zs * s2 - 0.5d0 * (zp * s2p + zm * s2m)))
      by    = 0.d0
      bz    = fy * (b0 * aln + b1 *(g1 - 0.5d0 * (g1p + g1m))
     >        + b2 * (g2 - 0.5d0 * (g2p + g2m)))
c
c     ---------------------------------------------------------
c*FON La contribution courante au calcul de la queue magnetique 
c*FON est terminee
c     ---------------------------------------------------------
c
      ex1 = exp(x / delx)
      ex2 = ex1**2
      y2  = y**2
      z2  = z**2
      yz  = y * z
c
      bx = bx + (ex1 * pa(1) + ex2 * pa(3)) * z * cps
     >     + (ex1 * pa(2) + ex2 * (pa(4) + pa(5) * y2
     >     + pa(6) * z2)) * sps
      by = (ex1 * pa(7) + ex2 * pa(9)) * yz * cps 
     >     + (ex1 * pa(8) + ex2 * (pa(10) + pa(11) * y2
     >     + pa(12) * z2)) * y * sps
      bz = bz+(ex1 *(pa(13) + pa(14) * y2 + pa(15) * z2)
     >     + ex2 * (pa(17) + pa(18) * y2 + pa(19) * z2)) 
     >     * cps + (ex1 * pa(16) + ex2 * (pa(20) + pa(21)
     >     * y2 + pa(22) * z2)) * z * sps
c
c     ------------------------------------------------------------
c*FON Les contributions dcf et fac ont ete ajoutees a bx, by et bz
c     ------------------------------------------------------------
c
      xsm  = x * cps - z * sps
      zsm  = x * sps + z * cps
      z2   = zsm**2
      rr2  = xsm**2 + y**2
      zn   = sqrt((rr2 + z2) / rrc2 + 4.d0)**5
      brsm = dstr * 3.d0*zsm/zn
      bzsm = dstr * (2.d0 * z2-rr2+8.d0*rrc2)/zn
      by1  = brsm * y
      bxsm = brsm * xsm
c
c     --------------------------------------------------
c*FON Reconstitution des composantes du champ magnetique
c     --------------------------------------------------
c
      bx = bx + bxsm * cps + bzsm * sps 
      bz = bz - bxsm * sps + bzsm * cps 
      by = by + by1 
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end