subroutine stepv (magin,magout,indgm,indval,year,tilt,
     >                 rggsm,rgsmg,x,y,z,ds,errin,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.30 - V 2.0
c*VER 03.01.06 - V 2.1
c*VER 01.01.07 - V 3.0
c*
c*AUT spec. CNES - JC KOSIK - avril 1996
c*AUT port. CISI 
c*
c*ROL Theme : Calculs de geophysique
c*ROL         Calcul d'un pas elementaire de la ligne de champ par la 
c*ROL         methode d'integration de Runge Kutta a pas varaiable et
c*ROL         les coefficients de Merson.
c*
c*PAR magin  (I) : type de champ magnetique interne
c*PAR magout (I) : type de champ magnetique externe
c*
c*PAR indgm  (I) : type d'indice geomagnetique Kp ou Ae
c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ
c*
c*PAR year   (I) : annee fractionnaire (>= 2000.)
c*
c*PAR tilt   (I) : angle de tilt (radians)
c*
c*PAR rggsm  (I) : matrice de passage du repere geographique au
c*PAR            : repere solaire magnetospherique
c*PAR rgsmg  (I) : matrice de passage du repere solaire
c*PAR            : magnetospherique au repere geographique
c*
c*PAR x    (I/O) : coordonnee geocentrique en x
c*PAR y    (I/O) : coordonnee geocentrique en y
c*PAR z    (I/O) : coordonnee geocentrique en z
c*
c*PAR ds     (I) : pas elementaire d'integration
c*
c*PAR errin  (I) : erreur maximale permise lors de l'integration
c*
c*PAR ier    (O) : code retour
c*
c*NOT magin      : 1 = dipole + g11 + h11
c*NOT magin      : 2 = champ DGRF 1945-1970
c*NOT magin      : 3 = champ DGRF 1970-1995
c*
c*NOT magout     : 0 = pas de champ externe
c*NOT magout     : 1 = Tsyganenko 87
c*NOT magout     : 2 = Tsyganenko 89  
c*NOT magout     : 3 = Kosik 97  
c*
c*NOT indgm      : 1 ---> indice geomagnetique Kp
c*NOT indval     : 1 ---> Kp =  0  , 0+
c*NOT indval     : 2 ---> Kp =  1- , 1 , 1+
c*NOT indval     : 3 ---> Kp =  2- , 2 , 2+
c*NOT indval     : 4 ---> Kp =  3- , 3 , 3+
c*NOT indval     : 5 ---> Kp =  4- , 4 , 4+
c*NOT indval     : 6 ---> Kp >  5-
c*
c*NOT indgm      : 2 ---> indice geomagnetique Ae
c*NOT indval     : 1 ---> Ae =    0 - 50
c*NOT indval     : 2 ---> Ae =   50 - 100
c*NOT indval     : 3 ---> Ae =  100 - 150
c*NOT indval     : 4 ---> Ae =  150 - 250
c*NOT indval     : 5 ---> Ae =  250 - 400
c*NOT indval     : 6 ---> Ae >= 400
c*
c*NOT La division par 3 du pas, permet des coefficients plus simples.
c*       
c*INF utilise   : deriv
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.05.30 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - extension au champ IGRF 2005
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 magin, magout, indgm, indval
      double precision year
      double precision tilt
      double precision rggsm(3,3), rgsmg(3,3)
      double precision x, y, z
      double precision ds
      double precision errin
      integer ier
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      double precision r11,r12,r13
c*LOC r11,r12,r13 : pas elementaires sur x, y et z (aux points x, y et z)
c
      double precision x2,y2,z2
c*LOC x2,y2,z2 : coordonnees cartesiennes du point au 2eme pas
c
      double precision r21,r22,r23
c*LOC r21,r22,r23 : pas elementaires sur x, y et z (aux points x2, y2 et z2)
c
      double precision x3,y3,z3
c*LOC x3,y3,z3 : coordonnees cartesiennes du point au 3eme pas
c
      double precision r31,r32,r33
c*LOC r31,r32,r33 : pas elementaires sur x, y et z (aux points x3, y3 et z3)
c
      double precision x4,y4,z4
c*LOC x4,y4,z4 : coordonnees cartesiennes du point au 3eme pas
c
      double precision r41,r42,r43
c*LOC r41,r42,r43 : pas elementaires sur x, y et z (aux points x4, y4 et z4)
c
      double precision x5,y5,z5
c*LOC x5,y5,z5 : coordonnees cartesiennes du point au 5eme pas
c
      double precision r51,r52,r53
c*LOC r51,r52,r53 : pas elementaires sur x, y et z (aux points x5, y5 et z5)
c
      double precision ds3,errcur
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
10    continue
c
      ds3 = ds / 3.d0
c
c     ----------------------------------------------------------
c*FON 1er appel au sous programme de calcul des pas elementaires
c     ----------------------------------------------------------
c
      call derivv(magin,magout,indgm,indval,year,tilt,
     >      rggsm,rgsmg,ds3,x,y,z,r11,r12,r13,ier)
c
      x2 = x + r11
      y2 = y + r12
      z2 = z + r13
c
c     ------------------------------------------------------------
c*FON 2ieme appel au sous programme de calcul des pas elementaires
c     ------------------------------------------------------------
c
      call derivv(magin,magout,indgm,indval,year,tilt,
     >               rggsm,rgsmg,ds3,x2,y2,z2,r21,r22,r23,ier)
c
      x3 = x + .5d0 * (r11 + r21)
      y3 = y + .5d0 * (r12 + r22)
      z3 = z + .5d0 * (r13 + r23)
c
c     ------------------------------------------------------------
c*FON 3ieme appel au sous programme de calcul des pas elementaires
c     ------------------------------------------------------------
c
      call derivv(magin,magout,indgm,indval,year,tilt,
     >               rggsm,rgsmg,ds3,x3,y3,z3,r31,r32,r33,ier)
c
      x4 = x +.375d0 * (r11 + 3.d0 * r31)
      y4 = y +.375d0 * (r12 + 3.d0 * r32)
      z4 = z +.375d0 * (r13 + 3.d0 * r33)
c
c     ------------------------------------------------------------
c*FON 4ieme appel au sous programme de calcul des pas elementaires
c     ------------------------------------------------------------
c
      call derivv(magin,magout,indgm,indval,year,tilt,
     >               rggsm,rgsmg,ds3,x4,y4,z4,r41,r42,r43,ier)
c
      x5 = x + 1.5d0 * (r11 -3.d0 * r31 + 4.d0 * r41)
      y5 = y + 1.5d0 * (r12 -3.d0 * r32 + 4.d0 * r42)
      z5 = z + 1.5d0 * (r13 -3.d0 * r33 + 4.d0 * r43)
c
c     ------------------------------------------------------------
c*FON 5ieme appel au sous programme de calcul des pas elementaires
c     ------------------------------------------------------------
c
      call derivv(magin,magout,indgm,indval,year,tilt,
     >           rggsm,rgsmg,ds3,x5,y5,z5,r51,r52,r53,ier)
c
      errcur = abs(r11 - 4.5d0 * r31 + 4.d0 * r41 -.5d0 * r51) +
     >         abs(r12 - 4.5d0 * r32 + 4.d0 * r42 -.5d0 * r52) +
     >         abs(r13 - 4.5d0 * r33 + 4.d0 * r43 -.5d0 * r53)
c
c     -------------------------------
c*FON Test sur la precision de calcul
c     -------------------------------
c
      if (errcur .lt. errin) then
         x = x + .5d0 * (r11 + 4.d0 * r41 + r51)
         y = y + .5d0 * (r12 + 4.d0 * r42 + r52)
         z = z + .5d0 * (r13 + 4.d0 * r43 + r53)
      else
         ds = ds * .5d0
         goto 10
      endif
c
      if (errcur.lt. errin * .04d0 .and. abs(ds) .lt. 1.33d0) then
         ds = ds * 1.5d0
      endif
c
c     **************** 
c     Fin de programme
c     **************** 
c
      return
      end