subroutine carmel (bb,xi,vl,beq,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.01 - 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 parametre L de Mac Ilwain a partir de I.
c*
c*PAR bb  (I) : module du champ magnetique (gauss)
c*
c*PAR xi  (I) : valeur de l'invariant I
c*
c*PAR vl  (O) : valeur de L
c*
c*PAR beq (O) : module du champ magnetique a l'equateur (gauss)
c*
c*PAR ier (O) : code de retour
c*
c*NOT ier     : sans objet
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.01 - 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 bb, xi
      double precision vl, beq
      integer ier
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      double precision xx,gg
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 la valeur de gg
c     -------------------------
c
      xx = log(((xi**3) * bb) / 0.311653d0)
c
      if (xx + 22.d0 .gt. 0.d0) then
c
         if (xx + 3.d0 .gt. 0.d0) then
c
            if (xx - 3.d0 .gt. 0.d0) then
c
               if (xx - 12.d0 .gt. 0.d0) then
c
                  if (xx - 23.d0 .gt. 0.d0) then
                     gg = xx - 3.0460681d0
                  else
                     gg = (((((2.8212095d-8 * xx - 3.8049276d-6)
     >                  * xx + 2.170224d-4) * xx - 6.7310339d-3)
     >                  * xx + .12038224d0) * xx - .18461796d0)
     >                  * xx + 2.0007187d0
                  endif
c
               else
                  gg = ((((((((6.3271665d-10 * xx - 3.958306d-8)
     >               * xx+9.9766148d-07) * xx -
     >               1.2531932d-5) * xx + 7.9451313d-5) * xx
     >               -3.2077032d-4) * xx + 2.1680398d-3)*
     >               xx + 1.2817956d-2) * xx + .43510529d0)
     >               * xx + .6222355d0
               endif
c
            else
               gg = ((((((((2.6047023d-10 * xx + 2.3028767d-9)
     >            * xx - 2.1997983d-8) * xx-
     >            5.3977642d-7) * xx - 3.3408822d-6) * xx +
     >            3.8379917d-5) * xx + 1.1784234d-3) *
     >            xx + 1.4492441d-2) * xx + .43352788d0)
     >            * xx + .6228644d0
            endif
c
         else
            gg = ((((((((-8.1537735d-14 * xx + 8.3232531d-13)
     >         * xx + 1.0066362d-9) * xx +
     >         8.1048663d-8) * xx + 3.2916354d-6) * xx +
     >         8.2711096d-5) * xx + 1.3714667d-3) *
     >         xx + .015017245d0) * xx + .43432642d0)
     >         * xx + .62337691d0
         endif
c
      else
         gg = .333338d0 * xx + 3.0062102d0
      endif
c
c     ----------------------------------------------
c*FON Calcul de la valeur du L et du module du champ
c     ----------------------------------------------
c
      vl  = (((1.0d0 + exp(gg)) * 0.311653d0) / bb)**(1.d0 / 3.d0)
c
      beq = 1.0d0 + exp(gg)
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end