subroutine bdbspst2k (tilt,rb,gg,rr,thetr,phir,bmxe,bmye,bmze,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 00.03.29 - V 1.1
c*VER 01.06.01 - 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 lointain. 
c*
c*PAR tilt   (I) : angle de tilt (radians)
c*
c*PAR rb     (I) : distance subsolaire normalisee de la magnetopause
c*PAR            : (rayons terrestres)
c*
c*PAR gg     (I) : coefficients de l'expansion spherique harmonique (6,6)
c*
c*PAR rr     (I) : distance radiale geocentrique (rayons terrestres)
c*PAR thetr  (I) : colatitude geocentrique (radians)
c*PAR phir   (I) : longitude geocentrique (radians)
c*
c*PAR bmxe   (O) : composante xgsm du champ lointain
c*PAR bmye   (O) : composante ygsm du champ lointain
c*PAR bmze   (O) : composante zgsm du champ lointain
c*
c*PAR ier    (O) : code de retour
c*
c*NOT ier        : sans objet
c*
c*NOT common     : util
c*
c*INF utilise    : vspvcar
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 1.1 - 00.03.29 - amelioration
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 tilt,rb
      double precision gg(6,6)
      double precision rr,thetr,phir
      double precision bmxe,bmye,bmze
      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 kmax
c*LOC kmax : nombre de pas de calcul (= 6)
c
      integer mtot
c*LOC mtot : nombre de pas calcules
c
      integer k,l,jj,m
c*LOC k,l,jj,m : indices de boucles
c
       double precision g(6,6)
c
       double precision br,bt,bp
c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
c
      double precision dbr(9),dbt(9),dbp(9)
c*LOC dbr,dbt,dbp : composantes radiales, tangentielles et azimuthales du champ
c
      double precision thetlim
c*LOC thetlim : limite max de la colatitude
c
      double precision sh(6,6),fn(6),fm(6)
      double precision cp(6),sp(6),p(6,6),dp(6,6),const(6,6)
      double precision cph,sph,ct,st
      double precision ar,par
      double precision temp
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 kmax /6/
c
      data thetlim /1.d-15/
c
      data sp(1)   /0.d0/
      data p(1,1)  /1.d0/
      data cp(1)   /1.d0/
      data dp(1,1) /0.d0/
c
      data sh(1,1) /0.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
      if (sh(1,1) .ne. -1.d0) then
         sh(1,1) = -1.d0
c
c     ---------------------------------------------------
c*FON Calcul des coefficients de Shmidt et des constantes
c     ---------------------------------------------------
c
         do 10 k = 1, kmax
            fn(k) = k - 1
            do 20 l = 1, k
               fm(l) = l - 1
               const(k,l) = dble((k - 2)**2 - (l - 1)**2)
     >                      / dble((2 * k - 3) * (2 * k - 5))
20          continue
10       continue
c
         do 30 k = 2, kmax
            sh(k,1) = dble(2 * k - 3) * sh(k-1,1) / dble(k - 1)
            jj = 2
            do 40 l = 2, k
               sh(k,l) = sh(k,l-1) * dsqrt(dble((k - l + 1) * jj)
     >                   / (k + l - 2))
               jj = 1
40          continue
30       continue
c
         do 50 k = 2, kmax
            g(k,1) = gg(k,1) * sh(k,1)
            do 60 l = 2, k
               g(k,l) = gg(k,l) * sh(k,l)
60          continue
50       continue
      endif
c      
      cph = dcos(phir)
      sph = dsin(phir)
c
      ct = cos(thetr)
      st = sin(thetr)
c
      if (thetr .le. thetlim) then
c
         ct  = 1.0d0
         st  = 1.d-15
         cph = 1.d0
         sph = 0.d0
c
      else if (thetr .ge. pi) then
c
         ct  = -1.d0
         st  = 1.d-15
         cph = 1.d0
         sph = 0.d0
c
      endif
c
      sp(2) = sph
      cp(2) = cph
c
      do 70 k = 3, kmax
c
         sp(k) = sp(2) * cp(k-1) + cp(2) * sp(k-1)
         cp(k) = cp(2) * cp(k-1) - sp(2) * sp(k-1)
c
70    continue
c
      p(2,1)  = ct
      dp(2,1) = -st
      p(2,2)  = st
      dp(2,2) = ct
c
      ar = 1.d0
c
      br =  ar * g(2,1) * p(2,1)  * cos(tilt)
      bt =  ar * g(2,1) * dp(2,1) * cos(tilt)
      bp =  0.d0
c
      m = 1
c
      dbr(m) = ar * p(2,1)  * cos(tilt)
      dbt(m) = ar * dp(2,1) * cos(tilt)
      dbp(m) = 0.d0
c
      do 90 k = 3, kmax
c
         ar  = ar * (rr / rb)
c
         do 80 l = 1, k
c
            mtot = k + l
c
            if (k .eq. l) then
               p(k,k)  = st * p(k-1,k-1)
               dp(k,k) = st * dp(k-1,k-1) + ct * p(k-1,k-1)
            else
               p(k,l)  = ct * p(k-1,l) - const(k,l) * p(k-2,l)
               dp(k,l) = ct * dp(k-1,l) - st * p(k-1,l)
     >                 - const(k,l) * dp(k-2,l)
            endif
c
            par = p(k,l) * ar
c
            if (mod(mtot,2) .eq. 1) then
c
               m = m + 1
               dbp(m) = 0.d0
c
               if (l .eq. 1) then
                  temp = 1.d0
               else
                  temp   = cp(l)
                  dbp(m) = -sp(l) * fm(l) * par * cos(tilt) / st
                  bp     = bp + g(k,l) * dbp(m)
               endif
               dbr(m) = temp * fn(k) * par * cos(tilt)
               br     = br + g(k,l) * dbr(m)
               dbt(m) = temp * dp(k,l) * ar * cos(tilt)
               bt     = bt + g(k,l) * dbt(m)
c
            endif
c
80       continue
c
90    continue
c
      call vspvcar(thetr,phir,br,bt,bp,bmxe,bmye,bmze,ier)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end