subroutine integ (arc,beg,bend,b,jep,eco,fi,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.06.05 - V 2.0
c*VER 03.01.06 - V 2.1
c*
c*AUT spec. Mac Ilwain - HAWITT
c*AUT adap. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Calcul du deuxieme invariant adiabatique I.
c*
c*PAR arc  (I) : tableau (dim = 200) des arcs calcules
c*
c*PAR beg  (I) : tableau (dim = 200) des rapports des champs aux
c*PAR          : extremites d'un segment
c*PAR bend (I) : tableau (dim = 200) des rapports des champs aux
c*PAR          : extremites d'un segment
c*
c*PAR b    (I) : tableau (dim = 200) des valeurs du champ
c*
c*PAR jep  (I) : nombre de points
c*
c*PAR eco  (I) : tableau (dim = 200) des valeurs de la fonction
c*PAR          : logarithmique du rapport des champs
c*
c*PAR fi   (O) : valeur de l'invariant I
c*
c*PAR ier  (O) : code de retour
c*
c*NOT jep      : pour jep points dans le calcul, il faut un nombre
c*NOT          : minimum de (jep + 1) valeurs dans les tableaux
c*
c*NOT ier      : sans objet
c*
c*NOT common   : util
c*
c*INF utilise  : sans objet
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.05 - correction de commentaires de code
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 arc(200), beg(200), bend(200), b(200)
      integer jep
      double precision eco(200)
      double precision fi
      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 kk
      integer i
c*LOC i : indice de boucles
c
      double precision a,x2,x3,asum,dn,bb,c,t,te,tb
      double precision arg1, arg2
c*LOC Variables de travail intermediaires
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
c     ---------------------------------
c*FON Calcul de l'invariant adiabatique
c     ---------------------------------
c
      kk = jep
      if ((kk - 4) .le. 0) then
c
         if ((kk - 4) .eq. 0) then
            kk = kk - 1
         else
            a    = b(kk-1) / b(2)
            x2   = b(kk) / b(2)
            x3   = b(kk+1) / b(2)
            asum = arc(kk) + arc(kk+1)
            dn   = arc(kk) * arc(kk+1) * asum
            bb   = (-a * arc(kk+1) * (arc(kk) + asum) + x2 * asum**2
     >             - x3 * arc(kk)**2) / dn
            c    = (a * arc(kk+1) - x2 * asum + x3 * arc(kk)) / dn
            fi   = pid * (1.d0 - a + bb * bb / (4.d0 * c))
     >             / sqrt (abs(c))
         endif
c
      else
c
         t  = sqrt(1.d0 - bend(2) / b(2))
         fi = (2.d0 * t - log((1.d0 + t) / (1.d0 - t))) / eco(2)
c
         if ((b(2) - bend(kk)) .gt. 0.0d0) then
            kk = kk + 1
         endif
c
         t  = sqrt(abs(1.d0 - beg(kk) / b(2)))
         fi = fi - (2.d0 * t - log((1.d0 + t) / (1.d0 - t))) / eco(kk)
         kk = kk - 1
c
         do 10 i = 3 ,kk
c
            arg2 = bend(i) / b(2)
            arg1 = 1.d0 - arg2
c
            if (arg1 .gt. 0.0d0) then
               te = sqrt(arg1)
            else
               te = 1.d-5
            endif
c
            arg2 = beg(i) / b(2)
            arg1 = 1.d0 - arg2
c
            if (arg1 .gt. 0.0d0) then
               tb = sqrt(arg1)
            else
               tb = 1.d-5
            endif
c
            if ((abs(eco(i)) - 2.d-5) .gt. 0.0d0) then
               fi = fi + (2.d0 * (te - tb) -
     >              log((1.d0 + te) * (1.d0 - tb) / ((1.d0 - te) *
     >              (1.d0 + tb)))) / eco(i)
            else
               fi = fi + ((te + tb) * (arc(i) + arc(i+1))) / 4.d0
            endif
c
10       continue
c
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end