subroutine flgalp (magin,magout,year,tilt,rggsm,rgsmg,tetdip,
     >                   phidip,rre,theto,phio,fl,xlamb,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 les
c*ROL         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 tetdip (I) : colatitude geocentrique du dipole (radians)
c*PAR phidip (I) : longitude geocentrique du dipole (radians)
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 de Y. Galperin
c*
c*PAR xlamb  (O) : latitude invariante (radians)
c*
c*PAR ier    (O) : indicateur d'erreur
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     :  1 = Tsyganenko 87
c*NOT magout     :  2 = Tsyganenko 89  
c*NOT magout     :  3 = Kosik 97  
c*
c*NOT xlamb      : si theto < PI/2 , xlamb est pris negatif (trace)
c*                  
c*NOT ier        :  0 = OK
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    : pconjr, 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 - mise a jour du modele de champ interne
c*HST version 4.0 - 13.10.10 - mise a jour du modele de champ interne
c*HST                          (DGRF05 et IGRF10)

c*HST version 5.0 - 17.02.23 - mise a jour du modele de champ interne

c*HST                           (DGRF10 et IGRF15)

c*HST version 6.0 - 20.07.15 - mise a jour du modele de champ interne
c*HST                           (DGRF15 et IGRF20)

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 tetdip, phidip
      double precision rre, theto, phio
      double precision fl, xlamb
      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 rmagc,thetc,phic
c*LOC rmagc,thetc,phic : coordonnees spheriques 
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
      ier1   = 0
      ier2   = 0
      ier3   = 0
c
      rmagc  = 999.d0
      thetc  = 999.d0 / deg
      phic   = 999.d0 / deg
      fl     = 999.d0
      xlamb  = 999.d0 / deg
      dir    = -1.d0
      rfin   = 1.d0
      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 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*FON Calcul du l de Mac Ilwain methode 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)
c
c     --------------------------------
c*FON Calcul de la latitude invariante
c     --------------------------------
c
         call invlat(fl,xlambr,ier3)
c
         if (ier3 .ne. 0) then    
            ier = 2
         else
c
          xlamb = xlambr * deg
c
            if (thetc .gt. pid) then
               xlamb = -xlamb
            endif
         endif
      else
c
         ier = 1
c
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end