subroutine lines (r1,r2,r3,b,arc,err,j,vp,vn,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
c*AUT adap. CNES - JC KOSIK - juin 1995
c*AUT port. CISI
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Calcul de la ligne de champ passant par les points.
c*ROL         Logiciel modifie pour tenir compte du champ total.
c*
c*PAR r1  (I/O) : valeur du champ magnetique B pour le 1 er  point
c*PAR r2  (I/O) : valeur du champ magnetique B pour le 2 eme point
c*PAR r3  (I/O) : coordonnees du 3eme point
c*
c*PAR b     (O) : champ magnetique
c*
c*PAR arc (I/O) : tableau des arcs
c*
c*PAR err   (I) : parametre d'erreurs
c*
c*PAR j     (O) : nombre de points de la ligne de force
c*
c*PAR vp    (O) : coordonnees du 2 eme point
c*PAR vn    (O) : coordonnees du 3 eme point
c*
c*PAR ier   (O) : code de retour
c*
c*NOT ier       : sans objet
c*
c*NOT common    : util
c*
c*INF utilise   : gsfc65
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 year
      double precision r1(3), r2(3), r3(3)
      double precision b(200)
      double precision arc(200)
      double precision err
      integer j
      double precision vn(3), vp(3)
      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 ilp,is,i
c*LOC ilp,is,i : indices de boucles
c
      integer ier1
c*LOC ier1 : code retour des modules appeles
c
      double precision br,bt,bp
c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
c
      double precision ra(3)
      double precision a1,a2,a3,aab,aa,ad,am,ao6,arcj,asum
      double precision bb,bd,cc,cd,cre,dd,dn
      double precision pre1,pre2,pre3,qrt,r,rbar,rt,sit,sna,x
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*FON Affectation des constantes
c     --------------------------
c
      data year /1965.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
      cre = 0.25d0
c
      if (err - 0.15625d0) 10,20,20
c
10    continue
c
      cre = (err**0.333333333d0)
c
20    continue
c
      a3  = arc(3)
      aab = abs (a3)
      sna = a3 / aab
      a1  = arc(1)
      a2  = arc(2)
      ao6 = a3 * a3 / 6.0d0
      j   = 3
      ilp = 1
      is  = 1
c
      go to 60
c
30    continue
c
      is   = 1
      j    = j + 1
      ao6  = a3 * a3 / 6.0d0
      arcj = a1 + a2 + a3
      ad   = (asum + a1) / aa
      bd   = asum / bb
      cd   = a1 / cc
c
40    continue
c
      do 50 i = 1 ,3
c
         dd = r1(i) / aa - r2(i) / bb + r3(i) / cc
         if (is .ne. 2) then
c
            rt    = r1(i) - (ad * r1(i) - bd * r2(i) + cd * r3(i) -
     >              dd * arcj) * arcj
            ra(i) = r1(i)
            r1(i) = r2(i)
            r2(i) = r3(i)
            r3(i) = rt
            vp(i) = vn(i)
c
         endif
c
         rbar  = (r2(i) + r3(i)) / 2.d0 - dd * ao6
         vn(i) = vp(i) + a3 * rbar
c
50    continue
c
60    continue
c
      if (vn(2)) 70,80,80
c
70    continue
c
      vn(2) = -vn(2)
c
80    continue
c
      if (vn(2) - pi) 100,100,90
c
90    continue
c
      vn(2) = dpi - vn(2)
c
      go to 80
c
100   continue
c
      if (vn(3)) 110,120,120
c
110   continue
c
      vn(3) = vn(3) + dpi
c
      go to 100
c
120   continue
c
      if (vn(3) - dpi) 140,140,130
c
130   continue
c
      vn(3) = vn(3) - dpi
c
      go to 120
c
140   continue
c
      if (is .eq. 2) go to 150
c
      sit  = abs(sin(vn(2)))
      pre1 = vn(1)
      pre2 = pre1 * vn(2)
      pre3 = pre1 * sit * vn(3)
      r    = vn(1)
c
c     --------------------------------------------
c*FON Calcul du champ magnetique d'origine interne
c     --------------------------------------------
c
      call gsfc65(year,r,vn(2),vn(3),br,bt,bp,b(j),ier1)
c
      r3(1) = br / b(j)
      dn    = b(j) * vn(1)
      r3(2) = bt / dn
      r3(3) = bp / (dn * sit)
      asum  = a3 + a2
      aa    = asum * a2
      bb    = a3 * a2
      cc    = asum * a3
      is    = 2
c
      go to 40
c
150   continue
c
      sit  = abs(sin(vn(2)))
      b(j) = b(j) * ((pre1 / vn(1))**3)
c
      qrt = .5d0 * abs(r3(1)) / (.1d0 + abs(r3(2) * vn(1)))
      x   = (abs(vn(1) - pre1) + qrt * abs(vn(1) * vn(2) - pre2) +
     >      abs(vn(1) * sit * vn(3) - pre3)) / 
     >      (aab * err * sqrt(1.d0 + qrt * qrt))
c
      if (ilp .ne. 2) go to 180
c
      if (x - 3.3d0) 180,160,160
c
160   continue
c
      a3   = a3 * 0.2d0 * (8.0d0 + x) / (0.8d0 + x)
      j    = j - 1
      ilp  = 3
      asum = a2 + a1
      aa   = asum * a1
      bb   = a2 * a1
      cc   = asum * a2
c
      do 170 i = 1, 3
c
         vn(i) = vp(i)
         r3(i) = r2(i)
         r2(i) = r1(i)
         r1(i) = ra(i)
c
170   continue
c
      go to 250
c
180   continue
c
      if (j - 200) 190,260,260
c
190   continue
c
      a1 = a2
      if (b(j) - b(2)) 200,200,260
c
200   continue
c
      ilp = 2
      a2  = a3
      a3  = a3 * .2d0 * (8.d0 + x) / (.8d0 + x)
      am  = (2.d0 - r3(2) * vn(1)) * vn(1) * cre
      if (abs(a3) - am) 220,220,210
c
210   continue
c
      a3 = sna * am
c
220   continue
c
      if (sna * r3(1) + .5d0) 230,230,250
c
230   continue
c
      am = -.5d0 * sna * vn(1) / r3(1)
      if (abs(a3) - am) 250,250,240
c
240   continue
c
      a3 = sna * am
c
250   continue
c
      arc(j+1) = a3
      aab      = abs(a3)
c
      go to 30
c
260    continue
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end