subroutine dlgalp (magin,magout,year,tilt,rggsm,rgsmg,
     >                   rre,theto,phio,fl,xlamb,rmagc,thetc,phic,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 17.02.23 - V 5.0
c*VER 20.07.15 - V 6.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 L de Y. Galperin.
c*ROL         On calcule la ligne de champ qui passe par le satellite
c*ROL         jusqu'au point de l'hemisphere nord avec le champ complet
c*ROL         puis, on calcule L avec le champ interne seulement et
c*ROL         les programmes de Mac Ilwain associes.
c*ROL         Est calculee aussi la latitude invariante.
C*
C*PAR magin  (I) : type de champ magnetique interne
c*
c*PAR magout (I) : type de champ magnetique externe
c*
c*PAR year   (I) : annee fractionnaire (>= 2015.)
c*
c*PAR tilt   (I) : angle de tilt (radians)
c*
c*PAR rggsm  (I) : matrice de passage du repere geographique au
c*PAR            : repere solaire magnetospherique
c*PAR rgsmg  (I) : matrice de passage du repere solaire
c*PAR            : magnetospherique au repere geographique
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 fl     (O) : parametre L de Y. Galperin.
c*
c*PAR xlamb  (O) : latitude invariante (radians)
c*
c*PAR rmagc  (O) : distance radiale du point conjugue (rayons terrestres)
c*PAR thetc  (O) : colatitude geocentrique du point conjugue (radians)
c*PAR phic   (O) : longitude geocentrique du point conjugue (radians)
c*
c*PAR ier    (O) : code de retour
c*
C*NOT magin      : 1 = dipole + g11 + h11
C*NOT magin      : 2 = dgrf2000
C*NOT magin      : 3 = dgrf2005
C*NOT magin      : 4 = dgrf2010
C*NOT magin      : 5 = dgrf2015
C*NOT magin      : 6 = igrf2020
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 xlamb      : si thetc > PI/2, xlamb est pris negatif (trace)
c*                  
c*NOT ier        : 1 = le calcul du point conjugue a echoue
c*NOT ier        : 2 = la valeur de L vaut zero et xlamb non 
c*NOT            : calculee
c*
c*NOT common     : util
c*
c*INF utilise    : dconjr, invar, invlat
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.05.30 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - adoption du champ interne igrf 2005
c*HST version 4.0 - 13.10.10 - adoption du champ interne igrf 2010
c*HST version 5.0 - 17.02.23 - adoption du champ interne igrf 2015
c*HST version 6.0 - 20.07.15 - adoption du champ interne igrf 2020
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, magout
      double precision year
      double precision tilt
      double precision rggsm(3,3), rgsmg(3,3)
      double precision rre, theto, phio
      double precision fl, xlamb
      double precision rmagc, thetc, phic
      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 indgm
c*LOC indgm : type d'indice geomagnetique (= 1 : Kp)
c
      integer indval
c*LOC indval : indice geomagnetique : niveau d'amplitude du champ (1 a 6)
c
      integer np
c*LOC np : nombre de points calcules
c
      integer ier1,ier2,ier3
c*LOC ier1,ier2,ier3 : codes retour des modules appeles
c
      double precision tr(500),tthet(500),tphi(500)
c*LOC tr,tthet,tphi : tableaux des coordonnees spheriques
c
      double precision dir
c*LOC dir : direction du trace des lignes de champ
c
      double precision rb
c*LOC rb : distance subsolaire de la magnetopause
c
      double precision rfin,rmax
c*LOC rfin,rmax : distances geocentriques d'arret et maximum autorisee
c
      double precision xlambr
c*PAR xlambr : latitude invariante en radians
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
c
      rmagc  = 999.d0
      thetc  = 999.d0
      phic   = 999.d0
      fl     = 999.d0
      xlamb  = 999.d0
      dir    = -1.d0
      rfin   = 1.0172d0
      rb     = 10.d0
      indgm  = 1
      indval = 1
      rmax   = 60.d0
c
c     ---------------------------------------------------------------
c*FON On choisit de tracer systematiquement vers le meme hemisphere
c*FON que celui ou est le point soit dir = -1.d0 , decision "Sauvaud"
c*FON Du 10 octobre 1991.
c     ---------------------------------------------------------------

c
      call dconjr(magin,year,magout,indgm,indval,tilt,rb,rggsm,rgsmg,
     >            dir,rre,theto,phio,rfin,rmax,np,tr,tthet,tphi,ier1)
c
c     ------------------------------------------------------------
c*FON Calcul du L de Galperin. On utilise invar
c*FON avec le champ interne seulement
c     ------------------------------------------------------------
c
      if (ier1 .eq. 0) then
c
         rmagc = tr(np)
         thetc = tthet(np)
         phic  = tphi(np)
c
         call invar(rmagc,thetc,phic,fl,ier2)
 
         call invlat(fl,xlambr,ier3)
c
         if (ier3 .ne. 0) then    
            ier = 2
         else
            xlamb  = xlambr * deg
            if (thetc .gt. pid) then
               xlamb = - xlamb
            endif
         endif
c
      else
         ier = 1
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end