subroutine ex89kp (indval,tilt,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 interne d'origine magnetospherique en
c*ROL         fonction de l'angle de tilt et de l'indice geomagnetique Kp.
c*
c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ
c*
c*PAR tilt   (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 en x dans le systeme solaire magnetospherique
c*PAR            : (gauss)
c*PAR by     (O) : composante en y dans le systeme solaire magnetospherique
c*PAR            : (gauss)
c*PAR bz     (O) : composante en z dans le systeme solaire magnetospherique
c*PAR            : (gauss)
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 Traite les composants GSM du champ magnetique produits par les 
c*NOT systemes extraterrestres courantsdans la geomagnetosphere. 
c*NOT Le modele est valide jusqu'a des distances geocentriques de 
c*NOT 70 rayons terrestres et bases sur l'union des jeux de donnees
c*NOT satellite imp-a,c,d,e,f,g,h,i,j (1966-1974) and heos-1,-2 
c*NOT (1969-1974).
c*NOT Reference: n.a. Tsyganenko, A magnetospheric
c*NOT magnetic field model with a warped tail current sheet: planet.
c*NOT space. sci., v.37, pp.5-20, 1989.
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*
      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 tilt
      double precision x, y, z
      double precision bx, by, bz
      integer ier
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer i
c*LOC i : indice de boucles
c
      integer ip
c*LOC ip : valeur initiale de indval (= 100)
c
      double precision ga1(28),ga2(28),ga3(28),ga4(28),ga5(28)
      double precision ga6(28),pa(28)
      double precision a,a02,aa4sps,adr,adrt,adrt2,adt2r2,at,brzr1
      double precision brzr2,bxc,bxcf,bxdr,bxt,byc,bycf,byt
      double precision bzc,bzcf,bzt,cps,d,d0,dadd,dadd05,dadd18
      double precision dd,ddel,ddop,ddopdx,ddr,ddx,ddy,del,delx
      double precision dfdrx,drdy2m,drdycm,dwcx,dwcy,dwx,dwy
      double precision dy,dyc,dzsx,dzsy,ex,f1,f3,f5,f7,f9,fcy,fdr
      double precision fk1,fk2,fr,fs,fxm,fxp,fy,fym,fyp,fzm,fzp
      double precision g,gam,gspm,h,ha02,hldxm,hlwc2m,hsx,htps
      double precision hxlw2m,hy,hys,psi,rc,rdy2,rdyc2,ro2,rq,rqc2
      double precision rsmrt,rsprt,rsqxdl,rt,rx2a2,s1,sm,sp,sps
      double precision srq,srqc2,srx2a2,sx,sxc,sxrc16,sxsix
      double precision t,tdr,w,wc,wcsm,wcsp,wt,wtfs,x2sm,xd
      double precision xld2,xlw2,xlwc2,xrc,xsixt,xsixtd,xsm
      double precision xsx,xsxc,xwcywc,xwyw,xxd, xxd2l2
      double precision y2,y4,y410,yz,z2,zm,zp,zr,zs,zsm
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 ga1/ -98.72d0,   -10014.d0,    15.03d0,    76.62d0,
     >         -10237.d0,     1.813d0,     31.1d0, -0.07464d0,
     >        -0.07764d0,  0.003303d0,   -1.129d0, 0.001663d0,
     >        0.000988d0,     18.21d0, -0.03018d0,       0.d0,
     >              0.d0,        0.d0,       0.d0,    24.74d0,
     >            8.16d0,      2.08d0,    -0.88d0,     9.08d0, 
     >            3.84d0,     13.55d0,    26.94d0,     5.75d0/
c
      data ga2/-35.65d0,   -12800.d0,    14.37d0,   124.5d0,
     >         -13543.d0,    2.316d0,    35.64d0, -0.0741d0,
     >         -0.1081d0, 0.003924d0,   -1.451d0, 0.00202d0,
     >         0.00111d0,    21.37d0, -0.04567d0,      0.d0,
     >              0.d0,       0.d0,       0.d0,   22.33d0,
     >            8.12d0,    1.664d0,    0.932d0,    9.24d0,
     >            2.43d0,    13.81d0,    28.83d0,    6.05d0/
c
      data ga3/   -77.45d0,  -14588.d0,    64.85d0,    123.9d0,
     >           -16229.d0,    2.641d0,    42.46d0, -0.07611d0,
     >           -0.1579d0, 0.004078d0,   -1.391d0,  0.00153d0,
     >          0.000727d0,    21.86d0, -0.04199d0,       0.d0,
     >                0.d0,       0.d0,       0.d0,     20.9d0,
     >              6.28d0,     1.54d0,     4.18d0,     9.61d0,
     >              6.59d0,    15.08d0,    30.57d0,     7.43d0/
c
      data ga4/  -70.12d0, -16125.d0,    90.71d0,    38.08d0,
     >          -19630.d0,   3.181d0,     47.5d0,  -0.1327d0,
     >          -0.1864d0, 0.01382d0,   -1.488d0, 0.002962d0,
     >         0.000897d0,   22.74d0, -0.04095d0,       0.d0, 
     >               0.d0,      0.d0,       0.d0,    18.64d0,
     >             6.27d0,   0.935d0,     5.39d0,     8.57d0,
     >             5.94d0,   15.63d0,    31.47d0,     8.10d0/
c
      data ga5/  -162.5d0, -15806.d0,   160.6d0,     5.888d0,
     >          -27534.d0,   3.607d0,     51.1d0,  -0.1006d0,
     >          -0.1927d0, 0.03353d0,   -1.392d0, 0.001594d0,
     >         0.002439d0,   22.41d0, -0.04925d0,       0.d0,
     >               0.d0,      0.d0,       0.d0,    18.31d0,
     >              6.2d0,   0.768d0,     5.07d0,    10.06d0,
     >             6.67d0,    16.1d0,    30.04d0,     8.26d0/
c
      data ga6/  -128.4d0, -16184.d0,    149.1d0,    215.5d0, 
     >          -36435.d0,    4.09d0,    49.09d0,  -0.0231d0,
     >          -0.1359d0, 0.01989d0,   -2.298d0, 0.004911d0,
     >         0.003421d0,   21.79d0, -0.05447d0,       0.d0,
     >               0.d0,      0.d0,       0.d0,    19.48d0,
     >             5.83d0,   0.332d0,     6.47d0,    10.47d0,
     >             9.08d0,   15.85d0,     25.27d0,    7.98d0/
c
      data del /0.01d0/
c
      data gam /4.d0/
c
      data dyc /20.d0/
c
      data xd /0.d0/
c
      data xld2 /40.d0/
c
      data xlw2 /170.d0/
c
      data a02 /25.d0/
c
      data rt /30.d0/
c
      data sxc /4.d0/
c
      data xlwc2 /50.d0/
c
      data dadd /1.d0/
c
      data ip /100/
c
      data psi /10.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
c     --------------------------------------------------------------
c*FON Calcul initial des constantes a partir d'un jeu donne de
c*FON parametres d'un modele specifie par le numero en option (indval)
c     --------------------------------------------------------------
c
c     cas ou indval est egal a 100
      if (indval.ne.ip) then
c
         ip = indval
c
         do 10 i = 1, 28
            if (ip .eq. 1) pa(i) = ga1(i)
            if (ip .eq. 2) pa(i) = ga2(i)
            if (ip .eq. 3) pa(i) = ga3(i)
            if (ip .eq. 4) pa(i) = ga4(i)
            if (ip .eq. 5) pa(i) = ga5(i)
            if (ip .eq. 6) pa(i) = ga6(i)
10       continue
c
         delx   = pa(20)
         adr    = pa(21)
         d0     = pa(22)
         dd     = pa(23)
         rc     = pa(24)
         g      = pa(25)
         a      = pa(26)
         dy     = pa(27)
         sx     = pa(28)
         ha02   = 0.5d0 * a02
         rdyc2  = 1.d0 / dyc**2
         hlwc2m = -0.5d0 * xlwc2
         drdycm = -2.d0 * rdyc2
         hldxm  = -0.5d0 * xld2
         ddel   = 2.d0 * del
         rdy2   = 1.d0 / dy**2
         drdy2m = -2.d0 * rdy2
         hxlw2m = -0.5d0 * xlw2
         dadd05 = dadd * 0.5d0
         dadd18 = -18.d0 * dadd
c
c     ---------------------------------------------------------
c*FON Coefficients pa(16)-pa(19) sont trouves dans pa(6)-pa(13)
c*FON afin que le champ magnetique soit moins divergent
c     ---------------------------------------------------------
c
         pa(16) = -0.5d0 * (pa(6) / delx + pa(10))
         pa(17) = -(pa(7) / delx + pa(11))
         pa(18) = -(pa(8) / delx + 3.d0 * pa(12))
         pa(19) = -(pa(9) / delx + pa(13)) / 3.d0
         psi    = tilt
         sps    = sin(tilt)
         cps    = cos(tilt)
         htps   = sps / (2.d0 * cps)
         gspm   = -g * sps
c
c     cas ou indval prend les autres valeurs
      else
c
         if (abs(tilt - psi) .lt. 1.d-6) then
c           aucun traitement
         else
            psi  = tilt
            sps  = sin(tilt)
            cps  = cos(tilt)
            htps = sps / (2.d0 * cps)
            gspm = -g * sps
         endif
c
      endif
c
c     ------------------------------------------------------------
c*FON Le traitement commence ici si indval et tilt sont identiques
c     ------------------------------------------------------------
c
      xsm    = x * cps - z * sps
      zsm    = z * cps + x * sps
      x2sm   = xsm**2
      y2     = y * y
      ro2    = x2sm + y2
      xxd    = xsm - xd
      xxd2l2 = 1.d0 / (xxd**2 + xld2)
      rsqxdl = sqrt(xxd2l2)
      h      = 0.5d0 * (1.d0 + xxd * rsqxdl)
      hsx    = - hldxm * xxd2l2 * rsqxdl
      xsixt  = xsm + 16.d0
      xsixtd = 1.d0 / (xsixt**2 + 36.d0)
      sxsix  = sqrt(xsixtd)
      ddop   = dadd05 * (1.d0 - xsixt * sxsix)
      ddopdx = dadd18 * xsixtd * sxsix
      d      = d0 + del * y2 + gam * h + ddop
      ddx    = gam * hsx + ddopdx
      ddy    = ddel * y
      xrc    = xsm + rc
      sxrc16 = sqrt(xrc**2 + 16.d0)
      y4     = y2 * y2
      y410   = 1.d0 / (y4 + 1.d4)
      hy     = y2 * y410 * y
      hys    = hy * y410 * 4.d4
      hy     = hy * y
      zs     = htps * (xrc - sxrc16) + gspm * hy
      dzsx   = htps * (1.d0 - xrc / sxrc16)
      dzsy   = gspm * hys
c
c     -------------------------------------------------------
c*FON zs = zs(xsm,ysm) definissent le format du champ courant
c     -------------------------------------------------------
c
      xsx   = xsm - sx
      rq    = 1.d0 / (xsx**2 + xlw2)
      srq   = sqrt(rq)
      fy    = 1.d0 /(1.d0 + y2 * rdy2)
      w     = 0.5d0 * (1.d0 - xsx * srq) * fy
      dwx   = hxlw2m * rq * srq * fy
      dwy   = drdy2m * w * y * fy
      zr    = zsm - zs
      t     = sqrt(zr**2 + d**2)
      at    = a +t
      s1    = sqrt(at**2 + ro2)
      f5    = 1.d0 / s1
      f7    = 1.d0 / (s1+ at)
      f1    = f5 * f7
      f3    = f5**3
      f9    = at *f3
      xwyw  = xsm * dwx + y * dwy
      fr    = zr * (xsm * dzsx + y * dzsy)
      fs    = fr - d * (xsm * ddx + y * ddy)
      wt    = w / t
      wtfs  = wt * fs
      brzr1 = wt * f1
      brzr2 = wt * f3
      bxt   = (pa(1) * brzr1 + pa(2) * brzr2) * zr
      byt   = bxt * y
      bxt   = bxt * xsm
      bzt   = pa(1) * (w * f5 + xwyw * f7 + wtfs * f1) +
     >        pa(2) * (w * f9 + xwyw * f1 + wtfs * f3)
c
c     ---------------------------------------------------------
c*FON La contribution de la page courante de la queue centrale 
c*FON (bxt,byt,bzt) est trouvee. Maintenant traitons les champs 
c*FON peripheriques 
c    ----------------------------------------------------------
c
      rx2a2  = 1.d0 / (x2sm + a02)
      srx2a2 = sqrt(rx2a2)
      fdr    = 0.5d0 * (1.d0 + xsm * srx2a2)
      dfdrx  = ha02 * rx2a2 * srx2a2
      ddr    = d0 + dd * fdr + ddop
      tdr    = sqrt(zr**2 + ddr**2)
      adrt   = adr + tdr
      adrt2  = adrt**2
      adt2r2 = 1.d0 / (adrt2 + ro2)
      fk1    = adt2r2**2 * sqrt(adt2r2)
      fk2    = 3.d0 * adrt * fk1 / tdr
      bxdr   = pa(5) * zr * fk2
      byt    = byt + bxdr * y
      bxt    = bxt + bxdr * xsm
      bzt    = bzt + pa(5) * ((2.d0 * adrt2 - ro2) * fk1 + 
     >         fk2 * (fr - ddr * (dd * dfdrx + ddopdx) * xsm))
c
c     ---------------------------------------------------
c*FON Calcul du champ Chapman-Ferraro et la moyenne de la 
c*FON contribution des champs alignes courants
c     ---------------------------------------------------
c
      ex     = exp(x / delx)
      z2     = z * z
      yz     = y * z
      bxcf   = ex * (cps * pa(6) * z + sps * (pa(7) + pa(8)
     >         * y2 + pa(9) * z2))
      bycf   = ex * (cps * pa(10) * yz + sps * y * 
     >         (pa(11) + pa(12) * y2 + pa(13) * z2))
      bzcf   = ex * (cps * (pa(14) + pa(15) * y2 + pa(16) * z2)
     >         + sps * z * (pa(17) + pa(18) * y2 + pa(19) * z2))
c
c     -----------------------------------------------------------
c*FON La queue magnetique renvoie les composants du champ courant
c*FON (bxc,byc,bzc)
c     -----------------------------------------------------------
c
      fcy    = 1.d0 / (1.d0 + y2 * rdyc2)
      xsxc   = x - sxc
      rqc2   = 1.d0 / (xsxc**2 + xlwc2)
      srqc2  = sqrt(rqc2)
      wc     = 0.5d0 * (1.d0 - xsxc * srqc2) * fcy
      dwcx   = hlwc2m * rqc2 * srqc2 * fcy
      dwcy   = drdycm * wc * y * fcy
      xwcywc = x * dwcx + y * dwcy
      ro2    = y2 + x**2
      zp     = z + rt
      zm     = z - rt
      sp     = sqrt(zp**2 + ro2)
      sm     = sqrt(zm**2 + ro2)
      wcsp   = wc / sp
      wcsm   = wc / sm
      rsprt  = 1.d0 / (sp + zp)
      rsmrt  = 1.d0 / (sm - zm)
      fxp    = wcsp * rsprt
      fxm    = - wcsm * rsmrt
      fyp    = fxp * y
      fym    = fxm * y
      fxp    = fxp * x
      fxm    = fxm * x
      fzp    = wcsp + xwcywc * rsprt
      fzm    = wcsm + xwcywc * rsmrt
      aa4sps = pa(4) * sps
      bxc    = pa(3) * (fxp + fxm) + (fxp - fxm) * aa4sps
      byc    = pa(3) * (fyp + fym) + (fyp - fym) * aa4sps
      bzc    = pa(3) * (fzp + fzm) + (fzp - fzm) * aa4sps
c
c     --------------------------------------------------
c*FON Somme des champs. Les composants centraux du champ 
c*FON courant sont transformes en coordonnees GSM
c     --------------------------------------------------
c
      bx = bxc + bxt * cps + bzt * sps + bxcf
      by = byc + byt + bycf
      bz = bzc + bzt * cps - bxt * sps + bzcf
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end