subroutine dconjr (magin,year,magout,indgm,indval,tilt,rb,
     >                   rggsm,rgsmg,dir,rre,theto,
     >                   phio,rfin,rmax,np,tr,tthet,tphi,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*VER 01.01.07 - V 3.0
c*VER 13.10.10 - V 4.0
c*VER 2017.02.23 - V 5.0
c*
c*AUT spec. CNES - JC KOSIK - juin 1995
c*AUT port. CISI 
c*AUT adapt. AKKA 
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Calcul du point conjugue d'un point donne en tenant
c*ROL         compte d'une combinaison de champs interns et externe.
c*
c*PAR magin  (I) : type de champ magnetique interne
c*
c*PAR year   (I) : annee fractionnaire
c*
c*PAR magout (I) : type de champ magnetique externe
c*
c*PAR indgm  (I) : type d'indice geomagnetique Kp ou Ae
c*PAR indval (I) : indice geomagnetique de 1 a 6
c*
c*PAR tilt   (I) : angle de tilt (radians)
c*
c*PAR rb     (I) : distance subsolaire normalisee de la magnetopause
c*PAR            : (rayons terrestres) rb = 10. * re
c*
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*
c*PAR dir    (I) : direction du trace des lignes de champ
c*PAR            : dir = -1 trace vers le meme hemisphere
c*PAR            : dir = +1 trace vers l'hemisphere oppose
c*
c*PAR rre    (I) : distance radiale geocentrique (rayons terrestres)
c*PAR theto  (I) : colatitude geocentrique (radians)
c*PAR phio   (I) : longitude geocentrique (radians)
c*
c*PAR rfin   (I) : distance geocentrique d'arret (rayons terrestres)
c*
c*PAR rmax   (I) : distance geocentrique maximum autorisee
c*PAR            : (rayons terrestres)
c*
c*PAR np     (O) : nombre de points calcules
c*
c*PAR tr     (O) : tableau des distances geocentriques des points de la
c*PAR            : ligne de champ (rayons terrestres)
c*PAR tthet  (O) : tableau des colatitudes des points de la
c*PAR            : ligne de champ (radians)
c*PAR tphi   (O) : tableau des longitudes des points de la
c*PAR            : ligne de champ (radians)
c*
c*PAR ier    (O) : code de retour
c*
c*NOT magin      : 1 = dipole + g11 + h11
c*NOT magin      : 2 = champ DGRF 2000
c*NOT magin      : 3 = champ DGRF 2005
c*NOT magin      : 4 = champ DGRF 2010
c*NOT magin      : 5 = champ IGRF 2015
c*
c*NOT year       : year doit etre >= 2000.
c*
c*NOT magout     : 0 = pas de champ externe
c*NOT magout     : 1 = Tsyganenko 87
c*NOT magout     : 2 = Tsyganenko 89
c*NOT magout     : 3 = Kosik 97
c*
c*NOT indgm      : 1 ---> indices geomagnetique Kp
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 indgm      : 2 ---> geomagnetique Ae
c*NOT indval     : 1 ---> Ae =    0 - 50
c*NOT indval     : 2 ---> Ae =   50 - 100
c*NOT indval     : 3 ---> Ae =  100 - 150
c*NOT indval     : 4 ---> Ae =  150 - 250
c*NOT indval     : 5 ---> Ae =  250 - 400
c*NOT indval     : 6 ---> Ae >= 400
c*
c*NOT dir        : +1 = trace vers les altitudes plus elevees
c*NOT            :      (depuis la surface vers l hemisphere oppose)
c*NOT            : -1 = trace vers les altitudes plus basses
c*NOT            :      (vers la surface du meme hemisphere
c*
c*NOT rfin       : 1 = dernier point sur Terre
c*
c*NOT rmax       : point d arret des lignes de champ ouvertes
c*
c*NOT ier        : 1 = l > 500 ou r > rmax ou (rre-1) <= 0.01
c*NOT            : 2 = sortie de magnetopause
c*
c*NOT common     : util
c*
c*INF utilise    : magtot, spcar, mpause, step, carsp
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.05.30 - Ajout indice Kosik 97
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - mise a jour des modeles de champ interne
c*HST version 4.0 - 13.10.10 - mise a jour des modeles de champ interne
c*HST                          (DGRF05 et IGRF10)
c*HST version 5.0 - 2017.02.23 - mise a jour des modeles de champ interne
c*HST                            (DGRF10 et IGRF15)
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
      integer magin
      double precision year
      integer magout, indgm, indval
      double precision tilt, rb
      double precision rggsm(3,3), rgsmg(3,3)
      double precision dir
      double precision rre, theto, phio
      double precision rfin, rmax
      integer np
      double precision tr(500), tthet(500), tphi(500)
      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 ipop
c*LOC ipop : indicateur de positionnement par rapport a la magnetosphere
c*
      integer ier1,ier2,ier3,ier4,ier5,ier6,ier7
c*LOC ier1,ier2,ier3,ier4,ier5,ier6,ier7 : codes retour des modules appeles
c
      double precision br,bt,bp
c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
c
      double precision x1,x2,x3
      double precision xxpp,xxp,xx,yypp,yyp,yy,zzpp
      double precision zzp,zz
c*LOC x1,x2,x3 : coordonnees cartesiennes x, y et z du point courant
c
      double precision rr,thetr,phir
c*LOC rr,thetr,phir : coordonnees spheriques
c
      double precision errin
c*LOC errin: erreur maximale permise lors de l'integration
c
      double precision tx(500),ty(500),tz(500)
c*LOC tx,ty,tz : tableau des coordonnees cartesiennes x, y et z 
c
      double precision a1,a2,a3,a4,a5,a6,a7
c*LOC a1,a2,a3,a4,a5,a6,a7 : termes des equations
c
      double precision altgc
c*LOC altgc : altitude en km du 1er point de la ligne de champ
c
      double precision epsr
c*LOC epsr : epsilon pour validite de l'altitude
c
      double precision ds0,rpp,r,b,sgn,ds,rp
c*LOC ds0,rpp,r,b,sgn,ds,rp : variables pour la determination du ds
c
      double precision px,py,pz
c*LOC px,py,pz : coordonees cartesiennes
c
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
      ier  = 0
      ier1 = 0
      ier2 = 0
      ier3 = 0
      ier4 = 0
      ier5 = 0
      ier6 = 0
      ier7 = 0
      ipop = 0
c
      rr    = rre
      thetr = theto
      phir  =  phio
c
c     -------------------
c*FON Determination du ds
c     -------------------
c
      altgc  = (rr - 1.0d0)
      epsr   = 0.001d0
      ds0    = 0.5d0
      if (dir .lt. 0.0d0) then
         if (altgc .gt. epsr .and. altgc .lt. 2.0d0) then
            ds0 = altgc / 20.0d0
         else if (altgc .le. epsr) then
            ier      = 0
            np       = 1
            tr(1)    = rr
            tthet(1) = thetr
            tphi(1)  = phir
         endif
      endif
      ds = ds0
c      
c     ---------------------------------------------------------------
c*FON Hemisphere Nord Br < 0, dir < 0, sgn > 0, trace dans le sens du 
c*FON champ, donc vers l'hemisphere Nord. Si dir > 0, sgn <0, trace
c*FON dans le sens oppose au champ, donc vers l'hemisphere Sud.
c*FON Hemisphere Sud : Br > 0, dir < 0, sgn < 0, trace dans le sens 
c*FON oppose au champ, donc vers l'hemisphere Sud. Si dir > 0, sgn > 0,
c*FON trace dans le sens du champ, donc vers l'hemisphere Nord.
c     ---------------------------------------------------------------
c
      call magtot(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,rr,
     >            thetr,phir,br,bt,bp,b,ier1)
c
      errin = 0.0005d0
      sgn   = sign(1.d0, br * dir)
      ds    = ds0 * sgn
      np    = 1
c
c     ----------------------------------------------------
c*FON Calcul de la position cartesienne du point de depart
c     ----------------------------------------------------
c
      call spcar(rr,thetr,phir,x1,x2,x3,ier2)
c
      tr(1)    = rr
      tthet(1) = thetr
      tphi(1)  = phir
c
      tx(np) = x1
      ty(np) = x2
      tz(np) = x3
c
c     -------------------------------------------------------------------
c*FON Calcul et test d'appartenance a la magnetosphere de type Shabansky,
c*FON si ipop = 2 point en dehors ipop = 0, point a l'interieur
c*FON Si le point de depart est a l'exterieur on ne fait pas les calculs 
c*FON et on met 999 dans les valeurs en 25
c     -----------------------------------------------
c
      call mpause(rb,rggsm,rr,thetr,phir,ipop,ier3)
c
      if (ipop .eq. 2) go to 10
c
40    continue
      np = np + 1
c
c     -------------------------------------------------------
c*FON Test de sortie de la magnetosphere, si ipop=2 on 
c*FON arrete les calculs sans mettre 999 aux valeurs
c*FON et on se debranche sur le return, on a trace un bout de
c*FON ligne de champ c'est tout
c     -------------------------------------------------------
c
      call mpause(rb,rggsm,rr,thetr,phir,ipop,ier4)
c
      if (ipop .eq. 2) go to 30
c
c     --------------------------------------------------
c*FON Integration pas a pas par appels successifs a step
c     --------------------------------------------------
c
      call step(magin,magout,indgm,indval,year,tilt,
     >              rggsm,rgsmg,x1,x2,x3,ds,errin,ier5)
c
      tx(np) = x1
      ty(np) = x2
      tz(np) = x3
c
c     ------------------------------------------------
c*FON Conversion des coordonnees cartesiennes du point
c*FON en coordonnees spheriques
c     ------------------------------------------------
c
      call carsp(x1,x2,x3,rr,thetr,phir,ier6)
c
      tr(np)    = rr
      tthet(np) = thetr
      tphi(np)  = phir
c
      if (rr .gt. rfin .and. np .lt. 500 .and. rr .lt. rmax) go to 40
c
      if (np .ge. 500 .or. rr .ge. rmax) then
         ier        = 1
         tr(np)     = 999.0d0
         tthet(np)  = 999.0d0 * rad
         tphi(np)   = 999.0d0 * rad
      endif
c
c     ----------------------------------------------------------
c*FON La limite rfin est depassee, on interpole pour obtenir les
c*FON coordonnees exactes du point en cartesiennes converties 
c*FON en coordonnees spheriques apres calcul
c     ----------------------------------------------------------
c
      if (rr .lt. rfin) then
c
         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)
c
         call carsp (px,py,pz,rr,thetr,phir,ier7)
c
         tr(np) = rr
         tthet(np) = thetr
         tphi(np) = phir
c
      endif
c
c     --------------------------------------------------------------
c*FON Si la distance limite rmax ou si le nbr de points depasse 500
c*FON ou si on arrive en dehors de la magnetosphere de Shabansky on
c*FON remplit les coordonnes par des 999 ou 999*rad pour obtenir des
c*FON 999 en degres apres la conversion radians degres
c     --------------------------------------------------------------
c
10    continue
c
      if (ipop .eq. 2) then
c       
         tr(np)    = 999.0d0
         tthet(np) = 999.0d0 * rad
         tphi(np)  = 999.0d0 * rad
         ier       = 2 
c       
      endif
30    continue
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end