subroutine pconjrs (ds,dir,tilt,indval,rggsm,rgsmg,xgsm,ygsm,zgsm, 
     >                    rfin,rmax,np,tx,ty,tz,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*
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Calcul de la ligne de champ. 
c*
c*PAR ds     (I) : pas de calcul initial (rayons terrestres)
c*PAR dir    (I) : direction du trace des lignes de force
c*PAR indval (I) : valeur de l'indice geomagnetique Kp 
c*PAR rggsm  (I) : matrice (3,3) de passage du repere geographique au
c*PAR            : repere solaire magnetospherique
c*PAR rgsmg  (I) : matrice (3,3) de passage du repere solaire
c*PAR            : magnetospherique au repere geographique
c*PAR xgsm   (I) : coordonnee solaire magnetospherique x (r. terr.)
c*PAR ygsm   (I) : coordonnee solaire magnetospherique y (r. terr.)
c*PAR zgsm   (I) : coordonnee solaire magnetospherique z (r. terr.)
c*PAR rfin   (I) : point d arret de la ligne de force (kilometres)
c*PAR rmax   (I) : point d arret des lignes de force ouvertes
c*PAR np     (O) : nombre de points calcules
c*PAR tr     (O) : tableau des distances geocentriques des point de la
c*PAR            : ligne de force
c*PAR tx     (O) : tableau des coordonnees x des points de la ligne de force
c*PAR ty     (O) : tableau des coordonnees y des points de la ligne de force
c*PAR tz     (O) : tableau des coordonnees z des points de la ligne de force
c*PAR ier    (O) : code de retour
c*
c*NOT indval     : varie de 1 a 6 pour Kp 
c*NOT dir        : +1 = trace vers les altitudes plus elevees
c*NOT            :      (depuis la surface ou vers l hemisphere oppose)
c*NOT            : -1 = trace vers les altitudes plus basses
c*NOT            :      (vers la surface du meme hemisphere)
c*NOT rfin       : point d arret de la ligne de force choisi
c*NOT            : pour la surface rfin = 1
c*NOT rmax       : point d arret des lignes de force ouvertes
c*NOT rfin       : point d arret de la ligne de force choisi
c*NOT            : pour la surface rfin = 1
c*
c*NOT ier        : 0 = OK
c*NOT ier        : 1 = l > 500 ou r > rmax ou (rre - 1) <= 0.01
c*NOT ier        : 2 = sortie de magnetopause
c*
c*NOT common     : sans objet
c*
c*INF utilise    : steps
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
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 ds
      double precision dir
      double precision tilt
      integer indval
      double precision rggsm(3,3), rgsmg(3,3)
      double precision xgsm,ygsm,zgsm
      double precision rfin
      double precision rmax
      integer np
      double precision tx(900),ty(900),tz(900),tr(900)
      integer ier
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
c errin: erreur maximale permise lors de l'integration
      double precision errin
      double precision x1, x2, x3
      double precision a1,a2,a3,a4,a5,a6,a7
      double precision rr, rpp, rp, r
      double precision xxpp, xxp, xx
      double precision yypp, yyp, yy
      double precision zzpp, zzp, zz
      double precision px, py, pz
cc
c     -----------------------------------------
c*FON Declaration de la fonction interne pinter
c     -----------------------------------------
c
      double precision pinter
c
      SAVE
c
c     ----------------------------------------
c*FON Definition de la fonction interne pinter
c     ----------------------------------------
c
      pinter(a1,a2,a3,a4,a5,a6,a7) = 
     >     ((a2 - a3) * (a7 - a2) * (a7 - a3) * a4
     >    - (a1 - a3) * (a7 - a1) * (a7 - a3) * a5 
     >    + (a1 - a2) * (a7 - a2) * (a7 - a1) * a6)
     >    / ((a1 - a2) * (a2 - a3) * (a1 - a3))
c
c     --------------------------------- 
c*FON Affectation identificateur rcs_id
c     --------------------------------- 
c
      data rcs_id /"
     >$id$"/
c
c     ******************
c     Debut de programme
c     ******************
c
c     ---------------
c*FON Initialisations
c     ---------------
c
      errin = 0.000001d0
      np    = 1
      x1    = xgsm
      x2    = ygsm
      x3    = zgsm
c
      tx(np) = x1
      ty(np) = x2
      tz(np) = x3
c
c     --------------------------------------
c*FON Boucle de traitement tant que r > rfin
c*FON et que     = 0
c     --------------------------------------
c
10    continue
      np = np + 1
c
c     ----------------------
c*FON Integration sur un pas
c     ----------------------
c
      call steps(ds,dir,indval,tilt,rggsm,rgsmg,x1,x2,x3,errin,ier)
c
      tx(np) = x1
      ty(np) = x2
      tz(np) = x3
c
      rr = dsqrt(x1**2 + x2**2 + x3**2)
c
      tr(np) = rr
c
      if (rr .gt. rfin .and. np .lt. 500 .and. rr .lt. rmax) go to 10
c
c     -------------
c*FON Fin de boucle
c     -------------
c
      ier = 1
      if (rr .lt. rfin) then
         ier    = 0
         rpp    = tr(np-2)
         rp     = tr(np-1)
         r      = tr(np)
         xxpp   = tx(np-2)
         xxp    = tx(np-1)
         xx     = tx(np)
         yypp   = ty(np-2)
         yyp    = ty(np-1)
         yy     = ty(np)
         zzpp   = tz(np-2)
         zzp    = tz(np-1)
         zz     = tz(np)
         px     = pinter(rpp,rp,r,xxpp,xxp,xx,rfin)
         py     = pinter(rpp,rp,r,yypp,yyp,yy,rfin)
         pz     = pinter(rpp,rp,r,zzpp,zzp,zz,rfin)
         tx(np) = px
         ty(np) = py
         tz(np) = pz
      endif 
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end