subroutine startlf (r1,r2,r3,b,arc,v,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 - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Initialisation du calcul des lignes de champ.
c*ROL         Le calcul est modifie pour pouvoir utiliser
c*ROL         un champ complet.
c*
c*PAR r1  (I/O) : champ b du 1er point calcule
c*PAR r2  (I/O) : champ b du 2eme point calcule
c*PAR r3    (O) : coordonnees du 3eme point calcule
c*
c*PAR b     (O) : module du champ magnetique (gauss)
c*
c*PAR arc   (O) : tableau des arcs calcules 
c*
c*PAR v     (O) : v(i,j) = i eme coordonnee du j eme point
c*
c*PAR ier   (O) : code de retour
c*
c*NOT arc       : en sortie arc(3) est change de signe de facon a
c*NOT           : ce que le point 1 soit plus pres de la terre que
c*NOT           : le point 2 sur la meme ligne de champ.
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 v(3,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 is,i
      integer ier1,ier2
c*LOC ier1,ier2 : codes 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 dn,r,sit
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
      ier1 = 0
      ier2 = 0
c
      r   = v(1,2)
      sit = abs(sin(v(2,2)))
c
10    continue
c
      if (v(3,2)) 20,30,30
c
20    continue
c
      v(3,2) = v(3,2) + dpi
c
      go to 10
c
30    continue
c
      call gsfc65(year,r,v(2,2),v(3,2),br,bt,bp,b(2),ier1)
c
      r2(1) = br / b(2)
      dn    = b(2) * v(1,2)
      r2(2) = bt / dn
      r2(3) = bp / (dn * sit)
      is    = 0
c
40    continue
c
      do 50 i = 1 , 3
         v(i,1) = v(i,2) - arc(2) * r2(i)
50    continue
c
60    continue
c
      r   = v(1,1)
      sit = abs(sin(v(2,1)))
c
      call gsfc65(year,r,v(2,1),v(3,1),br,bt,bp,b(1),ier2)
c
      if (b(1) - b(2)) 70,80,80
c
70    continue
c
      arc(2) = -arc(2)
c
      go to 40
c
80    continue
c
      r1(1)  = br / b(1)
      arc(3) = arc(2)
      dn     = b(1) * v(1,1)
      r1(2)  = bt / dn
      r1(3)  = bp / (dn * sit)
c
      do 90 i = 1, 3
         v(i,1) = v(i,2) - arc(2) * (r1(i) + r2(i)) / 2.d0
90    continue
c
      sit = abs(sin(v(2,1)))
      is  = is + 1
      if (is.eq.1) go to 60
c
c
      do 100 i = 1, 3
         v(i,3) = v(i,2) + arc(3) * ((1.5d0) 
     >            * r2(i) - .5d0 * r1(i))
100   continue
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end