subroutine ddparab (rb,x1,y1,z1,dd,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.05 - 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 : Frontieres et regions
c*ROL         Calcul de la distance d'un satellite a la parabole de
c*ROL         Shabansky.
c*
c*PAR rb  (I) : distance subsolaire de la magnetopause
c*PAR         : (rayons terrestres) rb = 11. * re
c*
c*PAR x1  (I) : coordonnee solaire magnetospherique ou ecliptique en x
c*PAR y1  (I) : coordonnee solaire magnetospherique ou ecliptique en y
c*PAR z1  (I) : coordonnee solaire magnetospherique ou ecliptique en z
c*
c*PAR dd  (O) : plus courte distance du satellite a la magnetopause
c*PAR         : parabolique de Shabansky
c*
c*PAR ier (O) : code de retour
c*
c*NOT ier     : sans objet
c*
c*INF utilise : kardan, distpar
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.05 - 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 rb
      double precision x1, y1, z1
      double precision dd
      integer ier
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer i
c*LOC i : indice de boucle
c
      integer icode
c*LOC icode : nombre de racines reelles trouvees
c
      integer ier1,ier2
c*LOC ier1,ier2 : codes retour des modules appeles
c
      double precision ddist
c*LOC ddist : distance a la parapole
c
      double precision a(4)
c*LOC a : tableau a des coefficients de l'equation
c
      double precision rr(3),ri(3)
c*LOC rr,ri : tableaux des racines reelles et imaginaires
c
      double precision td(3)
c*LOC td : tableau des distances a la parabole
c
      double precision x
c*LOC x : abscisse xgsm du satellite
c
      double precision ro1
c*LOC ro1 : distance radiale du satellite
c
      double precision ro
c*LOC ro : distance radiale du satellite de la magnetopause
c
      double precision calc,signe
c*LOC calc,signe : variables de travail intermediaires
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier  = 0
      ier1 = 0
      ier2 = 0
c
c     ----------------------------------------
c*FON Test de la valeur de rb par rapport a 0.
c     ----------------------------------------
c
      if (rb .ne. 0.0d0) then
c
c     ---------------------------------------------
c*FON Preparation des donnees pour l'appel a kardan
c     ---------------------------------------------
c
         ro1  = sqrt(y1**2 + z1**2)
         a(1) = 1.d0
         a(2) = 0.d0
         a(3) = +2.d0 * rb * x1
         a(4) = -2.d0 * rb * rb * ro1
         calc = 2.d0 * rb * x1 + ro1 * ro1 - 2.d0 * rb * rb
         signe = +1.d0
         if (calc .lt. 0.d0) then
            signe = -1.d0
         endif
c
c     ----------------------------------------------
c*FON Calcul des racines de l'equation du 3eme degre
c     ----------------------------------------------
c
         call kardan(a,rr,ri,icode,ier1)
c
c     -----------------------------------
c*FON Calcul de la distance a la parabole
c     -----------------------------------
c
         do 10 i = 1, icode
            ro = rr(i)
            x  = rb - (ro * ro / 2.d0 / rb)
            call distpar(x1,ro1,x,ro,ddist,ier2)
            td(i) = ddist
10       continue
c
         if (icode .eq. 1) then
            dd = td(1) * signe
         else if (icode .eq. 3) then
            dd = min(td(1),td(2),td(3)) * signe
         endif
c
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end