subroutine kk97kp (indval,tilt,xgsm,ygsm,zgsm,
     >                   bxt,byt,bzt,bt,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. CNES - JC KOSIK - janvier 1998
c*AUT port. CISI
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul du champ externe dans un systeme en coordonnees 
c*ROL         solaires magnetospheriques. Les unites sont les radians,  
c*ROL         rayons terrestres et nanoteslas.
c
c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ
c*
c*PAR tilt   (I) : angle de tilt (radians)
c*
c*PAR xgsm   (I) : coordonnee solaire magnetospherique x (rayons terrestres)
c*PAR ygsm   (I) : coordonnee solaire magnetospherique y (rayons terrestres)
c*PAR zgsm   (I) : coordonnee solaire magnetospherique z (rayons terrestres)
c*
c*PAR bxt    (O) : composante solaire magnetospherique x du champ 
c*PAR            : externe (nanoteslas)
c*PAR byt    (O) : composante solaire magnetospherique y du champ 
c*PAR            : externee (nanoteslas)
c*PAR bzt    (O) : composant solaire magnetospherique z du champ 
c*PAR            : externe (nanoteslas)
c*
c*PAR bt     (O) : totalite du champ (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*INF utilise  : carsp, k899k, vspvcar, tailts1k, bdbspst2k
c*INF utilise  : cylind2k
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
      double precision tilt,xgsm,ygsm,zgsm
      double precision bxt,byt,bzt,bt
      integer indval,ier 
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer i,j,m
c*LOC i,j,m : indices de boucles
c
      double precision gg(6,6)
c*LOC gg : coefficients du champ lointain
c
      double precision rb
c*LOC rb : distance radiale geocentrique
c
      double precision xn,bn,delb,s
c*LOC  xn,bn,delb,s : coefficients du modele de la queue
c
      double precision facrc,c20
c*LOC facrc,c20 : coefficients du courant de l'anneau
c
      double precision tx(5,17)
c*LOC tx : tx(1) ..to tx(9) gg(i,j) coefficients, tx(10) = rb
c*LOC      tx(11).....tx(14) coefficients xn, bn, delb, s
c*LOC      tx(15),tx(16) facrc, c20
c*LOC      tx(17) = bb
c
      double precision bxts,byts,bzts
c*LOC bxts,byts,bzts : composantes solaires magnetospheriques en x, y et z
c
      double precision bb
c*LOC bb : coefficient du systeme de courants en retour
c
      double precision bmxe,bmye,bmze
c*LOC bmxe,bmye,bmze : composantes xgsm, ygsm et zgsm du champ lointain
c
      double precision bcylx,bcyly,bcylz
c*LOC bcylx,bcyly,bcylz : composantes du champ magnetique induit par les
c*LOC                   : courants en retour en x, y et z
c
      double precision xsm,ysm,zsm
c*LOC xsm,ysm,zsm : coordonnees solaires magnetiques cartesiennes
c
      double precision rsm,thetsm,phism
c*LOC rsm,thetsm,phism : composantes solaires magnetiques spheriques
c
      double precision brsm,btsm,bpsm
c*LOC brsm,btsm,bpsm : vecteur du courant de l'anneau
c
      double precision bxsm,bysm,bzsm
c*LOC bxsm,bysm,bzsm : coordonnees solaires magnetiques cartesiennes
c*LOC                : du courant de l'anneau
c
      double precision bxgsm,bygsm,bzgsm
c*LOC bxgsm,bygsm,bzgsm : coordonnees GSM cartesiennes
c*LOC                   : du courant de l'anneau
c
      double precision rgsm,thgsm,phgsm
c*LOC rsm,thetsm,phism : composantes GSM spheriques
c
      double precision xmax,hmax,xgsmo,rcyl,xk,delz,zgsk
c*LOC xmax,hmax,xgsmo   : variables intermediaires pour le calcul
c*LOC rcyl,xk,delz,zgsk : du champ magnetique
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 (tx(1,m),m=1,17)/
     > -10.98617d0, -8.50034d0, 1.34113d0, -2.48786d0, 0.30964d0,
     >  -0.42732d0, -0.02118d0, 0.01296d0, -0.03237d0,     12.d0,
     >       -7.d0,      50.d0,     45.d0,      70.d0,     1.d0,
     >      -0.6d0,      25.d0/

      data (tx(2,m),m=1,17)/
     > -9.91009d0, -9.34152d0,  1.68348d0, -2.85354d0,  0.43415d0,
     > -0.47281d0, -0.02807d0,  0.02365d0, -0.03232d0,      12.d0,
     >	    -7.d0,      40.d0,      30.d0,      70.d0,      1.3d0,
     >     -0.6d0,      30.d0/

      data (tx(3,m),m=1,17)/
     > -12.37934d0, -10.70398d0,  1.86899d0, -3.45258d0, 0.52744d0,
     >  -0.59398d0,  -0.02654d0,  0.03671d0, -0.04219d0,     12.d0,
     >	     -7.d0,       40.d0,      30.d0,      70.d0,     1.5d0,
     >      -0.4d0,       35.d0/

      data (tx(4,m),m=1,17)/
     > -10.32178d0, -35.92150d0, 25.41001d0, -38.68030d0, 22.99485d0,
     > -20.01264d0,  -3.43988d0,  5.45664d0,  -3.66522d0,      36.d0,
     >       -7.d0,       40.d0,      30.d0,       70.d0,      1.8d0,
     >      -0.4d0,       50.d0/

      data (tx(5,m),m=1,17)/
     > -3.03703d0, 0.35692d0, 0.12408d0, 1.13281d0, 0.10293d0,
     >  0.35461d0, 0.00433d0, 0.01610d0, 0.03338d0,     10.d0,
     >      -7.d0,     50.d0,     40.d0,     70.d0,      2.d0,
     >      0.2d0,     60.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
c     ------------------------------------------------
c*FON Les coefficients gg(i,j) sont initialises a zero
c     ------------------------------------------------
c
      do 10 i = 2,6
         do 20 j = 1,i
            gg(i,j)  = 0.d0
20       continue
10    continue
c
c     ------------------------------
c*FON Coefficients du champ lointain
c     ------------------------------
c
      gg(2,1) = tx(indval,1)
      gg(2,2) = 0.d0
      gg(3,1) = 0.d0
      gg(3,2) = tx(indval,2)
      gg(3,3) = 0.d0
      gg(4,1) = tx(indval,3)
      gg(4,2) = 0.d0
      gg(4,3) = tx(indval,4)
      gg(4,4) = 0.d0
      gg(5,1) = 0.d0
      gg(5,2) = tx(indval,5)
      gg(5,3) = 0.d0
      gg(5,4) = tx(indval,6)
      gg(5,5) = 0.d0
      gg(6,1) = tx(indval,7)
      gg(6,2) = 0.d0
      gg(6,3) = tx(indval,8)
      gg(6,4) = 0.d0
      gg(6,5) = tx(indval,9)
      gg(6,6) = 0.d0
c
      rb      = tx(indval,10)
c
c     ----------------------------------
c*FON Coefficients du modele de la queue
c     ----------------------------------
c
      xn   = tx(indval,11)
      bn   = tx(indval,12)
      delb = tx(indval,13)
      s    = tx(indval,14)
c
c     ------------------------------------
c*FON Coefficients ddu courant de l'anneau
c     ------------------------------------
c
      facrc = tx(indval,15)
      c20   = tx(indval,16)
c
c     --------------------------------------------
c*FON Coefficient du systeme de courants en retour
c     --------------------------------------------
c
      bb   = tx(indval,17)
c
c     -------------------------------------------
c*FON Tous les composants sont initialises a zero
c     -------------------------------------------
c
      bxts = 0.d0
      byts = 0.d0
      bzts = 0.d0
c
      bmxe = 0.d0
      bmye = 0.d0
      bmze = 0.d0
c
      bcylx = 0.d0
      bcyly = 0.d0
      bcylz = 0.d0
c
c     -----------------------------------------
c*FON Calcul du champ k89 (courant de anneau) :
c     -----------------------------------------
c
c     ----------------------------------------------------------
c*FON 1) Les coordonnees GSM sont transformees en coordonnees SM
c     ----------------------------------------------------------
c
      xsm = xgsm * dcos(tilt) - zgsm * dsin(tilt)
      ysm = ygsm
      zsm = xgsm * dsin(tilt) + zgsm * dcos(tilt)
c
c     -------------------------------------------------------
c*FON 2) Les coordonnees cartesiennes SM sont transformees en 
c*FON    coordonnees spheriques SM 
c      ------------------------------------------------------
c 
      call carsp(xsm,ysm,zsm,rsm,thetsm,phism,ier)
c
c     ---------------------------------------------------
c*FON 3) Le champ courant de l'anneau est calcule dans le  
c*FON    systeme de coordonnees speriques SM
c     ---------------------------------------------------
c
      call k899k (tilt,facrc,c20,rsm,thetsm,phism,brsm,btsm,bpsm,ier)
c
c     -------------------------------------------------
c*FON 4) Les composants spheriques SM sont transformees
c*FON    en composants cartesiens SM  
c     -------------------------------------------------
c
      call vspvcar(thetsm,phism,brsm,btsm,bpsm,bxsm,bysm,bzsm,ier)
c
c     -----------------------------------------------------
c*FON 5) Cette contribution est transformee dans le systeme
c*FON    de coordonnees GSM
c     -----------------------------------------------------
c
      bxgsm =  bxsm * dcos(tilt) + bzsm * dsin(tilt)
      bygsm =  bysm
      bzgsm = -bxsm * dsin(tilt) + bzsm * dcos(tilt)
c
c     ---------------------------------------------------------------
c*FON Le champ queue est calcule en coordonnees systeme GSM. C'est le 
c*FON modele de queue 1982 Tsyganenko avec des parametres modifies
c     ---------------------------------------------------------------
c
      call tailts1k(tilt,xn,bn,delb,s,xgsm,ygsm,zgsm,bxts,byts,bzts,ier)
c
c     -------------------------------------
c*FON Calcul des coordonnees spheriques GSM  
c     -------------------------------------
c
      call carsp(xgsm,ygsm,zgsm,rgsm,thgsm,phgsm,ier)
c
c     ------------------------------------------
c*FON Calcul du champ magnetique lointain en GSM 
c     ------------------------------------------
c
      call bdbspst2k(tilt,rb,gg,rgsm,thgsm,phgsm,bmxe,bmye,bmze,ier)
c
c     -----------------------------------------------------------------
c*FON Calcul du champ magnetique a partir des courants en retour en GSM
c     -----------------------------------------------------------------
c
      xmax  = xn * cos(tilt)
      hmax  = dabs(xn) * sin(tilt)
      xgsmo = 0.d0
      rcyl  = 30.d0
      xk    = 0.05d0
c
      if (xgsm .ge. xmax) then
c
         delz = xgsm * hmax / xmax
         zgsk = zgsm - delz
         call cylind2k(xk,bb,xgsmo,rcyl,xgsm,ygsm,zgsk,
     >                 bcylx,bcyly,bcylz,ier)
c
      else if (xgsm .lt. xmax) then
c
         delz = hmax
         zgsk = zgsm - delz            
         call cylind2k(xk,bb,xgsmo,rcyl,xgsm,ygsm,zgsk,
     >                 bcylx,bcyly,bcylz,ier)
c
      endif
c
c     -----------------------------------------------------------------
c*FON Somme des contributions au champ magnetique a partir de l'anneau, 
c*FON de la queue, du champ lointain et des courants en retour
c     -----------------------------------------------------------------
c
      bxt = bxgsm + bxts + bmxe + bcylx
      byt = bygsm + byts + bmye + bcyly
      bzt = bzgsm + bzts + bmze + bcylz
c
c     -------------------------------
c*FON Calcul du module du champ total
c     -------------------------------
c
      bt = sqrt(bxt * bxt + byt * byt + bzt * bzt)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end