subroutine invar (rre,thet,phi,flg,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*
c*AUT spec. Mac Ilwain
c*AUT adap. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Calcul du L de Mac Ilwain. Le calcul de L est
c*ROL         effectue en utilisant le champ total dans ses
c*ROL         differentes possibilites : champ interne seulement,
c*ROL         champ interne + champ externe.
c*
c*PAR rre  (I) : distance radiale geocentrique (rayons terrestres)
c*PAR thet (I) : colatitude geocentrique (radians)
c*PAR phi  (I) : longitude geocentrique (radians)
c*
c*PAR flg  (O) : parametre L de Mac Ilwain
c*
c*PAR ier  (O) : code de retour
c*
c*NOT Le module invar a ete etabli par Mac Ilwain au debut des anneees
c*NOT 60. Ce module doit logiquement n'etre utilise que pour la 
c*NOT ceinture de radiations.
c*NOT Le jeu de coefficients correspond au champ interne Jenssen Cain 
c*NOT de 1965. Calcul fait pour 1965.d0
c*NOT L'interpolation presentee ici ne peut reussir qu'a des fins de 
c*NOT parametrage de lignes de champ. Elle ne peut surtout pas etre 
c*NOT utilisee pour les calculs de donnees de radiations.
c*
c*NOT invar a ete modifie (usage de magnet.) (sept 30, 62)
c*NOT f(phi) et v(1,2) ont ete modifies le 30 sept. 1962
c*NOT l'erreur sur l est de moins de 0.1 * er.
c*
c*NOT ier      : sans objet
c*
c*NOT common   : util
c*
c*INF utilise  : startlf, lines, integ, carmel
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.05.30 - colatitude superieure a pi
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
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
      double precision rre, thet, phi
      double precision flg
      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 jup,jep
      integer i,j
c*LOC i,j : indice de boucles
c
      integer ier1,ier2,ier3,ier4
c*LOC ier1,ier2,ier3,ier4 : codes retour des modules appeles
c
      double precision v(3,3),b(200),arc(200),vn(3),vp(3)
      double precision beg(200),bend(200),blog(200),eco(200)
      double precision r1(3),r2(3),r3(3)
      double precision asum,bco,bobz,cco,dclt,dco,dn,dx
      double precision err,flint,sa,sc,epsang
      double precision pinear
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     --------------------------
c*FON Affectation des constantes
c     --------------------------
c
      data err /0.01d0/
c
      data epsang /0.0174532d0/
c
      data pinear /3.124139454d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier  = 0
      ier1 = 0
      ier2 = 0
      ier3 = 0
      ier4 = 0
c
c     -----------------------------------
c*FON Traitement angle inferieur a epsang
c     -----------------------------------
c
      if (thet .lt. epsang) then
         thet = epsang
      endif
      if (thet .gt. pinear) then
         thet = pinear
      endif
c
c     ---------------------------
c*FON Calcul de la ligne de champ
c     ---------------------------
c
      v(1,2) = rre
      v(2,2) = thet
      v(3,2) = phi
      arc(1) = 0.d0
      arc(2) = v(1,2) * sqrt(err) * 0.3d0
      dclt   = pid - 0.2007d0 * cos(v(3,2) + 1.239d0)
c
      if (v(2,2) - dclt) 10,10,11
c
11    continue
      arc(2) = -arc(2)
10    continue
      call startlf(r1,r2,r3,b,arc,v,ier1)
c
      do 12 i = 1, 3
         vp(i) = v(i,2)
         vn(i) = v(i,3)
12    continue
c
      call lines(r1,r2,r3,b,arc,err,j,vp,vn,ier2)
c
      if (j - 200) 16,17,17
c
17    flg = -1.0d0
      go to 18
16    jup = j
c
      do 40 j = 1, jup
         arc(j)  = abs(arc(j))
         blog(j) = log(b(j))
40    continue
c
      jep = jup - 1
c
      do 21 j = 2, jep
         asum    = arc(j) + arc(j+1)
         dx      = blog(j-1) - blog(j)
         dn      = asum * arc(j) * arc(j+1)
         bco     = ((blog(j-1) - blog(j+1)) * arc(j)**2 - dx * asum**2)
     >             / dn
         cco     = (dx * arc(j+1) - (blog(j) - blog(j+1)) * arc(j)) / dn
         sa      = .75d0 * arc(j)
         sc      = sa + .25d0 * asum
         dco     = blog(j-1) - cco * sa * sc
         eco(j)  = bco + cco * (sa + sc)
         beg(j)  = exp(dco + eco(j) * .5d0 * arc(j))
         bend(j) = exp(dco + eco(j) * .5d0 * (asum + arc(j)))
21    continue
c
      beg(jup)  = bend(jep)
      bend(jup) = b(jup)
      eco(jup)  = (2.0d0 / arc(jup)) * log(bend(jup) / beg(jup))
c
      call integ(arc,beg,bend,b,jep,eco,flint,ier3)
c
27    continue
c
      call carmel(b(2),flint,flg,bobz,ier4)
c
18    continue
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end