subroutine steps (ds,dir,indval,tilt,rggsm,rgsmg,x,y,z,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* c*AUT spec. CNES - JC KOSIK - avril 1996 c*AUT port. CISI c* c*ROL Theme : Calculs de geophysique c*ROL Integration par une methode de Merson (methode RK de type c*ROL 3(4)). Le 4ieme pas permettant de mesurer l'erreur et c*ROL de corriger le pas d 'integration. c* c*PAR ds (I) : pas elementaire d'integration (rayons terrestres) c* c*PAR dir (I) : direction du trace des lignes de champ c* c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ 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 errin (I) : erreur maximale permise lors de l'integration c* c*PAR ier (O) : code retour c* c*NOT dir : +1 = trace vers les altitudes plus elevees c*NOT : (depuis la surface ou vers l hemisphere oppose) c*NOT : -1 = trace vers les altitudes plus basses c*NOT : (vers la surface du meme hemispher c* 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 ier : sans objet c* c*INF utilise : derivs 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* 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 ds,dir integer indval double precision tilt double precision rggsm(3,3),rgsmg(3,3) double precision x,y,z, errin integer ier c c --------------------------------- c*FON Declaration des variables locales c --------------------------------- c double precision ds3,r11,r12,r13,errcur double precision r21,r22,r23 double precision r31,r32,r33 double precision r41,r42,r43 double precision r51,r52,r53 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 ds3 = -ds / 3.d0 c c ---------------------------------------------------------- c*FON 1er appel au sous programme de calcul des pas elementaires c ---------------------------------------------------------- c call derivs(ds3,dir,indval,tilt,rggsm,rgsmg, > x,y,z,r11,r12,r13,ier) c c ------------------------------------------------------------ c*FON 2ieme appel au sous programme de calcul des pas elementaires c ------------------------------------------------------------ c call derivs(ds3,dir,indval,tilt,rggsm,rgsmg, > x + r11,y + r12,z + r13, > r21,r22,r23, > ier) c c ------------------------------------------------------------ c*FON 3ieme appel au sous programme de calcul des pas elementaires c ------------------------------------------------------------ c call derivs(ds3,dir,indval,tilt,rggsm,rgsmg, > x + .5d0 * (r11 + r21), > y + .5d0 * (r12 + r22), > z + .5d0 * (r13 + r23), > r31,r32,r33, > ier) c c ------------------------------------------------------------ c*FON 4ieme appel au sous programme de calcul des pas elementaires c ------------------------------------------------------------ c call derivs(ds3,dir,indval,tilt,rggsm,rgsmg, > x + .375d0 * (r11 + 3.d0 * r31), > y + .375d0 * (r12 + 3.d0 * r32), > z + .375d0 * (r13 + 3.d0 * r33), > r41,r42,r43, > ier) c c ------------------------------------------------------------ c*FON 5ieme appel au sous programme de calcul des pas elementaires c ------------------------------------------------------------ c call derivs(ds3,dir,indval,tilt,rggsm,rgsmg, > x + 1.5d0 * (r11 - 3.d0 * r31 + 4.d0 * r41), > y + 1.5d0 * (r12 - 3.d0 * r32 + 4.d0 * r42), > z + 1.5d0 * (r13 - 3.d0 * r33 + 4.d0 * r43), > r51,r52,r53, > ier) c errcur = dabs(r11 - 4.5d0 * r31 + 4.d0 * r41 - .5d0 * r51) + > dabs(r12 - 4.5d0 * r32 + 4.d0 * r42 - .5d0 * r52) + > dabs(r13 - 4.5d0 * r33 + 4.d0 * r43 - .5d0 * r53) c if (errcur .lt. errin) goto 20 c ds = ds * .5d0 c goto 10 c 20 continue c c ------------------------------- c*FON Test sur la precision de calcul c ------------------------------- c x = x + .5d0 * (r11 + 4.d0 * r41 + r51) y = y + .5d0 * (r12 + 4.d0 * r42 + r52) z = z + .5d0 * (r13 + 4.d0 * r43 + r53) c if (errcur .lt. errin * .04d0 .and. abs(ds) .lt. 1.d0) then ds = ds * 1.5d0 endif c c **************** c Fin de programme c **************** c return end