subroutine iterat (l,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.31 - V 2.0
c*VER 03.01.06 - V 2.1
c*
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Integration d'une ligne de force.
c*
c*PAR l   (I) : indicateur de traitement (nombre de pas)
c*
c*PAR ier (O) : code de retour
c*
c*NOT Integration d'une ligne de forces, similaire a la version FAP de
c*NOT "gl aide share pgm no. 413"). Utilise au depart une procedure
c*NOT "bootstrap", puis la formule Adams, avec lct = 1 au depart, 
c*NOT puis incrementation de lct de 1 chaque fois que phi contient les
c*NOT cooerdonnees initiales - distance geocentrique, latitude and l
c*NOT en radians. br, etc, sont les composants du champs a phi. 
c*NOT st = sin(theta) latitude a phi. Chaque pas genere un nouveau phi 
c*NOT et conserve le point phip precedant.
c*NOT Si sgn = +1, l'integration se deroule dans le sens direct
c*NOT c'est a dire , du nord vers le sud dans le champ terrestres.
c*NOT Si sgn = -1 la direction est inversee.
c*NOT d = valeur du pas le long de la ligne de forces r et dont la
c*NOT meme unite. 3 pas ont lieu en 7 iterations, puis 1 pas pour chaque  
c*NOT iteration supplentaire, yp(3,4) different de phi et 3 points 
c*NOT precedants.
c
c*NOT ier     : sans objet
c*
c*NOT common  : iter
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.05.31 - 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
      integer l
      integer ier
c
c     ----------------------------------
c*FON Declaration des variables communes
c     ----------------------------------
c
      double precision y(3),yold(3)
      double precision br,bt,bp,b,st,sgn,ds
c*COM y        : composantes spheriques
c*COM yold     : composantes spheriques precedantes
c*COM br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
c*COM b        : module du champ
c*COM st       : sinus theta
c*COM sgn      : sens du calcul (+1, -1)
c*COM ds       : intervalle d'integration le long de la ligne de champ
c
      common/iter/y,yold,br,bt,bp,b,st,sgn,ds
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      double precision yp(3,4)
      double precision d12,d2,d24,d6,fac
      integer i,j
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
      if (l .eq. 269) go to 200
      yp(1,4) = sgn * br / b
      fac     = sgn / (b * y(1))
      yp(2,4) = bt * fac
      yp(3,4) = bp * fac / st
c
      if (l .gt. 7) go to 90
c
c     --------------------------------------------------------
c*FON Bootstrapping pour les 7 premieres iterations, sur 3 pas
c     --------------------------------------------------------
c
      do 80 i = 1, 3
c
         go to (10,20,30,40,50,60,70),l
c
10       continue
c
         d2      = ds / 2.d0
         d6      = ds / 6.d0
         d12     = ds / 12.d0
         d24     = ds / 24.d0
         yp(i,1) = yp(i,4)
         yold(i) = y(i)
         y(i)    = yold(i) + ds * yp(i,1)
         go to 80
c
20       continue
c
         yp(i,2) = yp(i,4)
         y(i)    = yold(i) + d2 * (yp(i,2) + yp(i,1))
         go to 80
c
30       continue
c
         y(i) = yold(i) + d6 * 
     >          (2.d0 * yp(i,4) + yp(i,2) + 3.d0 * yp(i,1))
         go to 80
c
40       continue
c
         yp(i,2) = yp(i,4)
         yold(i) = y(i)
         y(i)    = yold(i) + d2 * (3.d0 * yp(i,2) - yp(i,1))
         go to 80
c
50       continue
c
         y(i) = yold(i) + d12 * 
     >          (5.d0 * yp(i,4) + 8.d0 * yp(i,2) - yp(i,1))
         go to 80
c
60       continue
c
         yp(i,3) = yp(i,4)
         yold(I) = y(i)
         y(I)    = yold(i) + d12 * (23.d0 * yp(i,3) - 16.d0 * yp(i,2)
     +             + 5.d0 * yp(i,1))
         go to 80
c
70       continue
c
         y(i) = yold(i) + d24 * (9.d0 * yp(i,4) + 19.d0 * yp(i,3)
     +          - 5.d0 * yp(i,2) + yp(i,1))
c
80    continue
      go to 200
c
c     ----------------------------------------------
c*FON 8eme iteration et apres -formule Adams 4-point 
c     ----------------------------------------------
c
90    continue
c
      do 100 i = 1, 3
c
         yold(i) = y(i)
         y(i)    = yold(i) + d24 * (55.d0 * yp(i,4) - 59.d0 * yp(i,3)
     +             + 37.d0 * yp(i,2) - 9.d0 * yp(i,1))
c
         do 110 j = 1, 3
c
            yp(i,j) = yp(i,j+1)
c
110      continue
c
100   continue
c
c     ****************
c     Fin de programme
c     ****************
c
200   continue
      return
      end