subroutine cylind2k (xk,bb,xgsmo,rcyl,xgsm,ygsm,zgsm,
     >                      bcylx,bcyly,bcylz,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.06.01 - V 2.0
c*VER 03.01.06 - V 2.1
c*
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul des composants du champ magnetique induit par
c*ROL         les courants en retour.
c*
c*PAR xk    (I) : coefficient d'exponentielle
c*
c*PAR bb    (I) : coefficient du champ cylindrique
c*
c*PAR xgsmo (I) : parametre xgsm
c*
c*PAR rcyl  (I) : rayon du cylindre pour la fonction de Bessel
c*
c*PAR xgsm  (I) : coordonnee solaire magnetospherique en x
c*              : (rayons terrestres)
c*PAR ygsm  (I) : coordonnee solaire magnetospherique en y
c*              : (rayons terrestres)
c*PAR zgsm (I)  : coordonnee solaire magnetospherique en z
c*              : (rayons terrestres)
c*
c*PAR bcylx (O) : composante en x du champ magnetique induit par les
c*PAR           : courants en retour
c*PAR bcyly (O) : composante en y du champ magnetique induit par les
c*PAR           : courants en retour
c*PAR bcylz (O) : composante en z du champ magnetique induit par les
c*PAR           : courants en retour
c*
c*PAR ier   (O) : code de retour
c*
c*NOT ier       : sans objet
c*
c*INF utilise   : bessel
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.01 - 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 xk, bb, xgsmo, rcyl
      double precision xgsm, ygsm, zgsm
      double precision bcylx, bcyly, bcylz
      integer ier 
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer kk
c*LOC bessel : ordre de la fonction bessel
c
      double precision racbess0, rgsm, cphi, sphi, zz, expz
      double precision dexp, xx, xj0, xj1, xj2, bcyr, bcyphi, bcyz
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 racbess0 /2.4d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
      rgsm = dsqrt(ygsm**2 + zgsm**2)
      if (rgsm .eq. 0.d0) then
         cphi = 1.d0
         sphi = 0.d0
      else
         cphi = ygsm / rgsm
         sphi = zgsm / rgsm
      endif
c
c     -----------------------------------------------------------
c*FON Normalisation des variables suivant les zeros de J0, J1, J2
c     -----------------------------------------------------------
c
      xx  = rgsm * racbess0 / rcyl
c
c*FON calcul de l'exponentielle
c
      zz  = dabs(xgsmo - xgsm)
      expz = dexp(xk * zz)
c
c     ----------------------------------------------------------------
c*FON Calcul des valeurs des fonctions de Bessel J0, J1 et J2 (kk=1,3)
c     ----------------------------------------------------------------
c
      kk = 1
c
      call bessel (xx,kk,xj0,ier)
c
      kk = 2
c
      call bessel (xx,kk,xj1,ier)
c
      kk = 3
c
      call bessel (xx,kk,xj2,ier)
c
      bcyr   =  0.5d0 * xk * bb * (xj0 - xj2) * sphi * expz
      bcyphi = -0.5d0 * xk * bb * (xj0 + xj2) * cphi * expz
      bcyz   =  xj1 * bb * sphi * expz
c
      bcylx  =  bcyz
      bcyly  =  bcyr * cphi - bcyphi * sphi
      bcylz  =  bcyr * sphi + bcyphi * cphi
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end