subroutine econjr (magin,year,magout,indgm,indval,tilt,rb,
     >                   rggsm,rgsmg,dir,rre,theto,
     >                   phio,rmax,np,tr,tthet,tphi,req,
     >                   theteq,phieq,ieq,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 2017.02.23 - V 4.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 equatorial d'un point donne en 
c*ROL         tenant compte d'une combinaison de champs internes et
c*ROL         externes.
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 : niveau d'amplitude du champ
c*
c*PAR tilt   (I) : angle de tilt (radians)
c*
c*PAR rb     (I) : distance subsolaire 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            
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 rmax   (I) : distance geocentrique maximum autorisee dans le trace
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 req    (O) : distance radiale au point de champ minimum
c*               : (rayons terrestres)
c*PAR theteq (O) : colatitude geocentrique du point de champ minimum (radians)
c*PAR phieq  (O) : longitude geocentrique du point de champ minimum (radians)
c*
c*PAR ieq    (O) : code de retour indiquant le passage a l'equateur
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 >= 2015.
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 ---> indice 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 ---> indice 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*
c*NOT rmax       : point d arret des lignes de champ ouvert
c*
c*NOT ieq        : 1 = passage a l equateur
c*NOT            : 0 = pas de passage a l equate
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 - 2017.02.23 - mise a jour des modeles de champ interne
c*HST                            (DGRF10 et IGRF15)
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, rre, theto, phio
      double precision rmax, bmin
      double precision req, theteq, phieq
      integer np
      double precision tr(500), tthet(500), tphi(500)
      integer ieq
      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
c*LOC ier1,ier2,ier3,ier4,ier5,ier6 : codes retour des modules appeles
c
      double precision br,bt,bp
c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ total
c
      double precision altgc
c*LOC altgc : altitude en km du 1er point de la ligne de champ
c
      double precision epsr
c*LOC espr : epsilon pour validite de l'altitude
c
      double precision rr,thetr,phir
c*LOC rr,thetr,phir : coordonnees spheriques
c
      double precision x1,x2,x3
c*LOC x1,x2,x3 : coordonnees cartesiennes du point courant
c
      double precision errin 
c*LOC errin : erreur maximale permise lors de l'integration
c
      double precision b,sgn,ds,ds0
c*LOC b,sgn,ds,ds0 : variables pour la determination du ds
c
      SAVE
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
      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 .gt. 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 dir = +1 dans l'hemisphere Nord, br < 0, sgn < 0, on trace dans
c*FON le sens oppose au champ, donc vers l'hemisphere Sud. Dans  
c*FON l'hemisphere Sud, br > 0, donc sgn > 0, on trace dans le sens du 
c*FON 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
c     --------------------------------------------------------
c*FON Calcul et test d'appartenance a la magnetosphere de type 
c*FON Shabansky, si ipop=2 point en dehors, si ipop=0, point a 
c*FON l'interieur.
c*FON Si le point de depart est a l'exterieur on ne fait pas 
c*FON les calculs 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) then
         go to 20
      endif
c
10    continue
      np = np + 1
c
c     -------------------------------------------------------------
c*FON Test de sortie de la magnetosphere, si ipop=2 on arrete les 
c*FON calculs sans mettre 999 aux valeurs et on se debranche sur le 
c*FON return, on a trace un bout de ligne de champ c'est tout.
c     -------------------------------------------------------------
c
      call mpause(rb,rggsm,rr,thetr,phir,ipop,ier4)
c
      if (ipop .eq. 2) then
         go to 30
      endif
c
c     -------------------------------------------------------------
c*FON On calcule le champ magnetique au point courant et on definit
c*FON les quantites equatoriales temporaires qui seront definitives
c*FON quand on ne passera plus dans la boucle
c     -------------------------------------------------------------
c
      call magtot(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,rr,
     >             thetr,phir,br,bt,bp,b,ier1)
c
      bmin       = b
      ieq        = 1
      req        = rr
      theteq     = thetr
      phieq      = phir
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
c     ---------------------------------------------------
c*FON Conversion des coordonnees cartesiennes du point en
c*FON coordonnees spheriques
c     ---------------------------------------------------
c
      call carsp(x1,x2,x3,rr,thetr,phir,ier6)
c
      tr(np)    = rr
      tthet(np) = thetr
      tphi(np)  = phir
c
      call magtot(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,rr,
     >            thetr,phir,br,bt,bp,b,ier1)
c
      if (b .le. bmin .and. np .lt. 500 .and. rr .lt. rmax) then
         go to 10
      endif
c
      if (np .ge. 500 .or. rr .ge. rmax .and. b .lt. bmin) then 
         ier        = 1
         ieq        = 0
         tr(np)     = 999.0d0
         tthet(np)  = 999.0d0 * rad
         tphi(np)   = 999.0d0 * rad
      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
c*FON des 999 en degres apres la conversion radians degres
c     -------------------------------------------------------------
c
20    continue
c
      if (ipop .eq. 2) then       
         tr(np)     = 999.0d0
         tthet(np)  = 999.0d0 * rad
         tphi(np)   = 999.0d0 * rad
         ier        = 2
         ieq        = 0 
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
30    continue
c
      return
      end