subroutine pconjr (magin,year,magout,indgm,indval,tilt,rb,rggsm,
     >                   rgsmg,tetdip,phidip,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.24- V 5.0
c*
c*AUT spec. CNES - JC KOSIK - janvier 1991
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 internes et externes.
c*
c*PAR magin  (I) : type de champ interne
c*
c*PAR year   (I) : annee fractionnaire de calcul
c*
c*PAR magout (I) : type de champ externe
c*
c*PAR indgm  (I) : type d'indice geomagnetique
c*PAR indval (I) : indice geomagnetique pour Kp ou Ae
c*
c*PAR tilt   (I) : angle de tilt (r adians)
c*
c*PAR rb     (I) : distance subsolaire de la magnetopause
c*PAR            : (rayons terrestres) rb = 11. * 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 tetdip (I) : colatitude geocentrique du dipole (radians)
c*PAR phidip (I) : longitude geocentrique du dipole (radians)
c*
c*PAR dir    (I) : direction du trace des lignes de force
c*
c*PAR rre    (I) : distance radiale geocentrique du point de depart
c*PAR            : (rayons terrestres)
c*
c*PAR theto  (I) : colatitude geocentrique du point de depart (radians)
c*PAR phio   (I) : longitude geocentrique du point de depart (radians)
c*
c*PAR rfin   (I) : distance geocentrique d'arret (rayons terrestres)
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 point de la
c*PAR            : ligne de force
c*PAR tthet  (O) : tableau des colatitudes des points calcules de la
c*PAR            : ligne de force
c*PAR tphi   (O) : tableau des longitudes des points calcules de la
c*PAR            : ligne de force
c*
c*PAR ier    (O) : code de retour
c*
c*NOT magin      : 1 = dipole
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 ---> 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 ou 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       : 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*
c*NOT ieq        : 1 = passage a l'equateur
c*NOT            : 0 = non passage a l'equateur
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     : iter
c*NOT common     : util
c*
c*INF utilise    : magtot, iterat, mpause
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 - utilisation du champ 2005
c*HST version 4.0 - 13.10.10 - introduction du champ igrf10 et passage 
c*HST                          de igrf05 à dgrf05 
c*HST version 5.0 - 2017.02.24 - introduction du champ igrf15 et passage 
c*HST                            de igrf10 à dgrf10 
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
      double precision rb
      double precision rggsm(3,3), rgsmg(3,3)
      double precision tetdip, phidip, dir
      double precision rre, theto, phio
      double precision rfin
      double precision rmax
      integer np
      double precision tr(500)
      double precision tthet(500)
      double precision tphi(500)
      integer ier
c
c     ----------------------------------
c*FON Declaration des variables communes
c     ----------------------------------
c
      double precision y(3),yold(3)
      double precision br,bt,bp,b,st,sgn,ds
c*COM y        : composantes spheriques
c*COM yold     : composantes spheriques precedantes
c*COM br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
c*COM b        : module du champ
c*COM st       : sinus theta
c*COM sgn      : sens du calcul (+1, -1)
c*COM ds       : intervalle d'integration le long de la ligne de champ
c
      common/iter/y,yold,br,bt,bp,b,st,sgn,ds
c
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 l
c*LOC l : pas d'integration
c
      integer ier1,ier2,ier3
c*LOC ier1,ier2,ier : codes retour des modules appeles
c
      double precision r,thet,phi,rp,thetp,phip
c*LOC r,thet,phi,rp,thetp,phip : coordonnees spheriques
c
      double precision phipp,pphi,pphid,pthet,rpp,thetpp
c*LOC phipp,pphi,pphid,pthet,rpp,thetpp : coordonnees spheriques
c
      double precision cglond,singml
c*LOC cglond,singml : variables intermediaires pour la determination du ds
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,epsr,epsang,ds0     
c*LOC Variables de travail intermediaires
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*FON Affectation des constantes
c     --------------------------
c
      data epsr/0.001d0/
      data epsang/0.00174532d0/
c
c     ******************
c     Debut de programme
c     ******************
c
c     ---------------
c*FON Initialisations
c     ---------------
c
      ier  = 0
      ier1 = 0
      ier2 = 0
      ier3 = 0
      ipop = 0
c
      sgn    = dir
c
      altgc  = (rre - 1.0d0)
      if (theto .lt. epsang) then
         theto = epsang
      endif
c
c     -------------------
c*FON Determination du ds
c     -------------------
c
      cglond = cos(phio - phidip)
      singml = cos(tetdip) * cos(theto) + sin(tetdip) 
     >         * sin(theto) * cglond
      ds0    = .06d0 * rre / (1.d0 - singml**2) - 0.058d0
      if (ds0 .gt. 0.3d0) ds0 = 0.3d0
      if (dir .lt.0. 0d0) then
         if (altgc .gt. epsr .and. altgc .lt. 2.0d0) then
            ds0 = altgc / 20.0d0
         elseif (altgc .le. epsr) then
            ier      = 0
            np       = 1
            tr(1)    = rre
            tthet(1) = theto
            tphi(1)  = phio
            goto 30
         endif
      endif
      ds   = ds0
c
      r    = rre
      phi  = phio
      thet = theto
      l    = 0
c
c     --------------------------------------
c*FON Boucle de traitement tant que r > rfin
c*FON et que     = 0
c     --------------------------------------
c
10    continue
c
c        ---------------------
c*FON    Calcul du champ total
c        ---------------------
c
         st = sin(thet)
         call magtot(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,r,
     >               thet,phi,br,bt,bp,b,ier1)
c
c        ----------------------------------
c*FON    Memorisation du champ a l'equateur
c        ----------------------------------
c
         if (l .eq. 0) then
            sgn = sign(1.d0,br * dir)
         endif
c
         l = l + 1
c
         thetpp  = thetp
         phipp   = phip
         rpp     = rp
c
c        ----------------------------------------
c*FON    Integration le long de la ligne de champ
c        ----------------------------------------
c
         if (l .gt. 30 .and. tr(l) .lt. 2.d0) then
            ds = ds0 / 2.d0
         endif
c
         y(1) = r
         y(2) = thet
         y(3) = phi
         yold(1) = rp
         yold(2) = thetp
         yold(3) = phip
         call iterat(l,ier2) 
         r    = y(1)
         thet = y(2)
         phi  = y(3)
         rp    = yold(1)
         thetp = yold(2)
         phip  = yold(3)
c   
         if (l .ge. 500 .or. rp .gt. rmax) then
            ier  = 1
            ipop = 1
         endif
c
         if (ier .eq. 1) go to 20
c
c        -----------------------------------------------
c*FON    Calcul et test d'appartenance a la magnetopause
c*FON    De type Shabansky
c        -----------------------------------------------
c
         call mpause(rb,rggsm,rp,thetp,phip,ipop,ier3)
c
         if (ipop .eq. 2) then
            ier = 2
            go to 20
         endif
c
         tr(l) = rp
c
         if (thet .gt. pi) then
            thet = dpi - thet
            phi  = pi  + phi
         endif
c 
         if (thetp .gt. pi) then
            thetp = dpi - thetp
            phip  = pi  + phip
         endif
c
         tthet(l) = thetp
         tphi(l)  = phip
c
      if (r .gt. rfin .and. ipop .eq. 0) go to 10
c
c     -------------
c*FON Fin de boucle
c     -------------
c
      pthet     = pinter(rpp,rp,r,thetpp,thetp,thet,rfin)
      pphi      = pinter(rpp,rp,r,phipp,phip,phi,rfin)
      pphi      = pphi + dpi
      pphi      = mod(pphi,dpi)
      np        = l + 1
      tr(np)    = rfin
      tthet(np) = pthet
      tphi(np)  = pphi
c
20    continue
c
c*    -----------------------------------------------------
c*FON Initialisation des parametes en sortie a des valeurs
c*FON arbitraires en cas de non calcul de la ligne de force
c*    -----------------------------------------------------
c
      if (ipop .eq. 1 .or. ipop .eq. 2) then
c
         np        = l
         pthet     = 999.0d0 * rad
         pphid     = 999.0d0 * rad
c
         tr(np)    = 999.0d0
         tthet(np) = pthet
         tphi(np)  = pphid
c
      endif
c
30    continue
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end