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