CDPP / AMDA_COTS

Commit 510e044c14a20ad7b361323119d7a512584e9e6c

Authored by Hacene SI HADJ MOHAND
1 parent bacb77ca
Exists in master

extension validite 1945-2025

Showing 8 changed files with 1759 additions and 9 deletions   Show diff stats
Inline Side-by-side
maglib/src/dconjrv.f 0 → 100755
  Diff comments   View file @510e044
... ... @@ -0,0 +1,420 @@
  1 + subroutine dconjrv (magin,year,magout,indgm,indval,tilt,rb,
  2 + > rggsm,rgsmg,dir,rre,theto,
  3 + > phio,rfin,rmax,np,tr,tthet,tphi,ier)
  4 +c*
  5 +c***********************************************************************
  6 +c*
  7 +c* "Copyright [c] CNES 98 - tous droits reserves"
  8 +c* **********************************************
  9 +c*
  10 +c*PRO MAGLIB
  11 +c*
  12 +c*VER 99.03.31 - V 1.0
  13 +c*VER 01.05.30 - V 2.0
  14 +c*VER 03.01.06 - V 2.1
  15 +c*VER 01.01.07 - V 3.0
  16 +c*VER 13.10.10 - V 4.0
  17 +c*VER 17.02.23 - V 5.0
  18 +c*VER 20.07.15 - V 6.0
  19 +c*
  20 +c*AUT spec. CNES - JC KOSIK - juin 1995
  21 +c*AUT port. CISI
  22 +c*AUT adapt. AKKA
  23 +c*
  24 +c*ROL Theme : Calculs de geophysique
  25 +c*ROL Calcul du point conjugue d'un point donne en tenant
  26 +c*ROL compte d'une combinaison de champs interns et externe.
  27 +c*
  28 +c*PAR magin (I) : type de champ magnetique interne
  29 +c*
  30 +c*PAR year (I) : annee fractionnaire
  31 +c*
  32 +c*PAR magout (I) : type de champ magnetique externe
  33 +c*
  34 +c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae
  35 +c*PAR indval (I) : indice geomagnetique de 1 a 6
  36 +c*
  37 +c*PAR tilt (I) : angle de tilt (radians)
  38 +c*
  39 +c*PAR rb (I) : distance subsolaire normalisee de la magnetopause
  40 +c*PAR : (rayons terrestres) rb = 10. * re
  41 +c*
  42 +c*PAR rggsm (I) : matrice (3,3) de passage du repere geographique au
  43 +c*PAR : repere solaire magnetospherique
  44 +c*PAR rgsmg (I) : matrice (3,3) de passage du repere solaire
  45 +c*PAR : magnetospherique au repere geographique
  46 +c*
  47 +c*PAR dir (I) : direction du trace des lignes de champ
  48 +c*PAR : dir = -1 trace vers le meme hemisphere
  49 +c*PAR : dir = +1 trace vers l'hemisphere oppose
  50 +c*
  51 +c*PAR rre (I) : distance radiale geocentrique (rayons terrestres)
  52 +c*PAR theto (I) : colatitude geocentrique (radians)
  53 +c*PAR phio (I) : longitude geocentrique (radians)
  54 +c*
  55 +c*PAR rfin (I) : distance geocentrique d'arret (rayons terrestres)
  56 +c*
  57 +c*PAR rmax (I) : distance geocentrique maximum autorisee
  58 +c*PAR : (rayons terrestres)
  59 +c*
  60 +c*PAR np (O) : nombre de points calcules
  61 +c*
  62 +c*PAR tr (O) : tableau des distances geocentriques des points de la
  63 +c*PAR : ligne de champ (rayons terrestres)
  64 +c*PAR tthet (O) : tableau des colatitudes des points de la
  65 +c*PAR : ligne de champ (radians)
  66 +c*PAR tphi (O) : tableau des longitudes des points de la
  67 +c*PAR : ligne de champ (radians)
  68 +c*
  69 +c*PAR ier (O) : code de retour
  70 +c*
  71 +c*NOT magin : 1 = dipole + g11 + h11
  72 +c*NOT magin : 2 = champ DGRF 1945-1970
  73 +c*NOT magin : 3 = champ DGRF 1970-1995
  74 +c*NOT magin : 3 = champ DGRF 1995-2000
  75 +c*
  76 +c*NOT year : year doit etre >= 1945.
  77 +c*
  78 +c*NOT magout : 0 = pas de champ externe
  79 +c*NOT magout : 1 = Tsyganenko 87
  80 +c*NOT magout : 2 = Tsyganenko 89
  81 +c*NOT magout : 3 = Kosik 97
  82 +c*
  83 +c*NOT indgm : 1 ---> indices geomagnetique Kp
  84 +c*NOT indval : 1 ---> Kp = 0 , 0+
  85 +c*NOT indval : 2 ---> Kp = 1- , 1 , 1+
  86 +c*NOT indval : 3 ---> Kp = 2- , 2 , 2+
  87 +c*NOT indval : 4 ---> Kp = 3- , 3 , 3+
  88 +c*NOT indval : 5 ---> Kp = 4- , 4 , 4+
  89 +c*NOT indval : 6 ---> Kp > 5-
  90 +c*
  91 +c*NOT indgm : 2 ---> geomagnetique Ae
  92 +c*NOT indval : 1 ---> Ae = 0 - 50
  93 +c*NOT indval : 2 ---> Ae = 50 - 100
  94 +c*NOT indval : 3 ---> Ae = 100 - 150
  95 +c*NOT indval : 4 ---> Ae = 150 - 250
  96 +c*NOT indval : 5 ---> Ae = 250 - 400
  97 +c*NOT indval : 6 ---> Ae >= 400
  98 +c*
  99 +c*NOT dir : +1 = trace vers les altitudes plus elevees
  100 +c*NOT : (depuis la surface vers l hemisphere oppose)
  101 +c*NOT : -1 = trace vers les altitudes plus basses
  102 +c*NOT : (vers la surface du meme hemisphere
  103 +c*
  104 +c*NOT rfin : 1 = dernier point sur Terre
  105 +c*
  106 +c*NOT rmax : point d arret des lignes de champ ouvertes
  107 +c*
  108 +c*NOT ier : 1 = l > 500 ou r > rmax ou (rre-1) <= 0.01
  109 +c*NOT : 2 = sortie de magnetopause
  110 +c*
  111 +c*NOT common : util
  112 +c*
  113 +c*INF utilise : magtot, spcar, mpause, step, carsp
  114 +c*
  115 +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
  116 +c*HST version 2.0 - 01.05.30 - Ajout indice Kosik 97
  117 +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
  118 +c*HST version 3.0 - 01.01.07 - mise a jour des modeles de champ interne
  119 +c*HST version 4.0 - 13.10.10 - mise a jour des modeles de champ interne
  120 +c*HST (DGRF05 et IGRF10)
  121 +c*HST version 5.0 - 17.02.23 - mise a jour des modeles de champ interne
  122 +c*HST (DGRF10 et IGRF15)
  123 +c*HST version 5.0 - 20.07.15 - mise a jour des modeles de champ interne
  124 +c*HST (DGRF15 et IGRF20)
  125 +c*
  126 +c***********************************************************************
  127 +c*
  128 + implicit none
  129 +c
  130 +c ---------------------------------
  131 +c*FON Declaration identificateur rcs_id
  132 +c ---------------------------------
  133 +c
  134 + character rcs_id*100
  135 +c
  136 +c --------------------------
  137 +c*FON Declaration des parametres
  138 +c --------------------------
  139 +c
  140 + integer magin
  141 + double precision year
  142 + integer magout, indgm, indval
  143 + double precision tilt, rb
  144 + double precision rggsm(3,3), rgsmg(3,3)
  145 + double precision dir
  146 + double precision rre, theto, phio
  147 + double precision rfin, rmax
  148 + integer np
  149 + double precision tr(500), tthet(500), tphi(500)
  150 + integer ier
  151 +c
  152 +c ----------------------------------
  153 +c*FON Declaration des variables communes
  154 +c ----------------------------------
  155 +c
  156 + double precision pi,dpi,rad,deg,pid,xmu,rayt
  157 +c
  158 +c*COM pi : constante pi (obtenue a partir de acos(-1.))
  159 +c*COM dpi : constante 2 * pi
  160 +c*COM pid : constante pi / 2
  161 +c*COM rad : facteur de conversion degres ----> radians
  162 +c*COM deg : facteur de conversion radians ----> degres
  163 +c*COM xmu : constante de gravitation terrestre (km**3/sec**2)
  164 +c*COM rayt : rayon equatorial terrestre (km)
  165 +c
  166 + common/util/pi,dpi,rad,deg,pid,xmu,rayt
  167 +c
  168 +c ---------------------------------
  169 +c*FON Declaration des variables locales
  170 +c ---------------------------------
  171 +c
  172 + integer ipop
  173 +c*LOC ipop : indicateur de positionnement par rapport a la magnetosphere
  174 +c*
  175 + integer ier1,ier2,ier3,ier4,ier5,ier6,ier7
  176 +c*LOC ier1,ier2,ier3,ier4,ier5,ier6,ier7 : codes retour des modules appeles
  177 +c
  178 + double precision br,bt,bp
  179 +c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
  180 +c
  181 + double precision x1,x2,x3
  182 + double precision xxpp,xxp,xx,yypp,yyp,yy,zzpp
  183 + double precision zzp,zz
  184 +c*LOC x1,x2,x3 : coordonnees cartesiennes x, y et z du point courant
  185 +c
  186 + double precision rr,thetr,phir
  187 +c*LOC rr,thetr,phir : coordonnees spheriques
  188 +c
  189 + double precision errin
  190 +c*LOC errin: erreur maximale permise lors de l'integration
  191 +c
  192 + double precision tx(500),ty(500),tz(500)
  193 +c*LOC tx,ty,tz : tableau des coordonnees cartesiennes x, y et z
  194 +c
  195 + double precision a1,a2,a3,a4,a5,a6,a7
  196 +c*LOC a1,a2,a3,a4,a5,a6,a7 : termes des equations
  197 +c
  198 + double precision altgc
  199 +c*LOC altgc : altitude en km du 1er point de la ligne de champ
  200 +c
  201 + double precision epsr
  202 +c*LOC epsr : epsilon pour validite de l'altitude
  203 +c
  204 + double precision ds0,rpp,r,b,sgn,ds,rp
  205 +c*LOC ds0,rpp,r,b,sgn,ds,rp : variables pour la determination du ds
  206 +c
  207 + double precision px,py,pz
  208 +c*LOC px,py,pz : coordonees cartesiennes
  209 +c
  210 +c -----------------------------------------
  211 +c*FON Declaration de la fonction interne pinter
  212 +c -----------------------------------------
  213 +c
  214 + double precision pinter
  215 +c
  216 + SAVE
  217 +c
  218 +c ----------------------------------------
  219 +c*FON Definition de la fonction interne pinter
  220 +c ----------------------------------------
  221 +c
  222 + pinter(a1,a2,a3,a4,a5,a6,a7)=
  223 + > ((a2 - a3) * (a7 - a2) * (a7 - a3) * a4
  224 + > - (a1 - a3) * (a7 - a1) * (a7 - a3) * a5
  225 + > + (a1 - a2) * (a7 - a2) * (a7 - a1) * a6)
  226 + > / ((a1 - a2) * (a2 - a3) * (a1 - a3))
  227 +c
  228 +c ---------------------------------
  229 +c*FON Affectation identificateur rcs_id
  230 +c ---------------------------------
  231 +c
  232 + data rcs_id /"
  233 + >$Id$"/
  234 +c
  235 +c ******************
  236 +c Debut de programme
  237 +c ******************
  238 +c
  239 +c ---------------
  240 +c*FON Initialisations
  241 +c ---------------
  242 +c
  243 + ier = 0
  244 + ier1 = 0
  245 + ier2 = 0
  246 + ier3 = 0
  247 + ier4 = 0
  248 + ier5 = 0
  249 + ier6 = 0
  250 + ier7 = 0
  251 + ipop = 0
  252 +c
  253 + rr = rre
  254 + thetr = theto
  255 + phir = phio
  256 +c
  257 +c -------------------
  258 +c*FON Determination du ds
  259 +c -------------------
  260 +c
  261 + altgc = (rr - 1.0d0)
  262 + epsr = 0.001d0
  263 + ds0 = 0.5d0
  264 + if (dir .lt. 0.0d0) then
  265 + if (altgc .gt. epsr .and. altgc .lt. 2.0d0) then
  266 + ds0 = altgc / 20.0d0
  267 + else if (altgc .le. epsr) then
  268 + ier = 0
  269 + np = 1
  270 + tr(1) = rr
  271 + tthet(1) = thetr
  272 + tphi(1) = phir
  273 + endif
  274 + endif
  275 + ds = ds0
  276 +c
  277 +c ---------------------------------------------------------------
  278 +c*FON Hemisphere Nord Br < 0, dir < 0, sgn > 0, trace dans le sens du
  279 +c*FON champ, donc vers l'hemisphere Nord. Si dir > 0, sgn <0, trace
  280 +c*FON dans le sens oppose au champ, donc vers l'hemisphere Sud.
  281 +c*FON Hemisphere Sud : Br > 0, dir < 0, sgn < 0, trace dans le sens
  282 +c*FON oppose au champ, donc vers l'hemisphere Sud. Si dir > 0, sgn > 0,
  283 +c*FON trace dans le sens du champ, donc vers l'hemisphere Nord.
  284 +c ---------------------------------------------------------------
  285 +c
  286 + call magtotv(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,rr,
  287 + > thetr,phir,br,bt,bp,b,ier1)
  288 +c
  289 + errin = 0.0005d0
  290 + sgn = sign(1.d0, br * dir)
  291 + ds = ds0 * sgn
  292 + np = 1
  293 +c
  294 +c ----------------------------------------------------
  295 +c*FON Calcul de la position cartesienne du point de depart
  296 +c ----------------------------------------------------
  297 +c
  298 + call spcar(rr,thetr,phir,x1,x2,x3,ier2)
  299 +c
  300 + tr(1) = rr
  301 + tthet(1) = thetr
  302 + tphi(1) = phir
  303 +c
  304 + tx(np) = x1
  305 + ty(np) = x2
  306 + tz(np) = x3
  307 +c
  308 +c -------------------------------------------------------------------
  309 +c*FON Calcul et test d'appartenance a la magnetosphere de type Shabansky,
  310 +c*FON si ipop = 2 point en dehors ipop = 0, point a l'interieur
  311 +c*FON Si le point de depart est a l'exterieur on ne fait pas les calculs
  312 +c*FON et on met 999 dans les valeurs en 25
  313 +c -----------------------------------------------
  314 +c
  315 + call mpause(rb,rggsm,rr,thetr,phir,ipop,ier3)
  316 +c
  317 + if (ipop .eq. 2) go to 10
  318 +c
  319 +40 continue
  320 + np = np + 1
  321 +c
  322 +c -------------------------------------------------------
  323 +c*FON Test de sortie de la magnetosphere, si ipop=2 on
  324 +c*FON arrete les calculs sans mettre 999 aux valeurs
  325 +c*FON et on se debranche sur le return, on a trace un bout de
  326 +c*FON ligne de champ c'est tout
  327 +c -------------------------------------------------------
  328 +c
  329 + call mpause(rb,rggsm,rr,thetr,phir,ipop,ier4)
  330 +c
  331 + if (ipop .eq. 2) go to 30
  332 +c
  333 +c --------------------------------------------------
  334 +c*FON Integration pas a pas par appels successifs a step
  335 +c --------------------------------------------------
  336 +c
  337 + call stepv(magin,magout,indgm,indval,year,tilt,
  338 + > rggsm,rgsmg,x1,x2,x3,ds,errin,ier5)
  339 +c
  340 + tx(np) = x1
  341 + ty(np) = x2
  342 + tz(np) = x3
  343 +c
  344 +c ------------------------------------------------
  345 +c*FON Conversion des coordonnees cartesiennes du point
  346 +c*FON en coordonnees spheriques
  347 +c ------------------------------------------------
  348 +c
  349 + call carsp(x1,x2,x3,rr,thetr,phir,ier6)
  350 +c
  351 + tr(np) = rr
  352 + tthet(np) = thetr
  353 + tphi(np) = phir
  354 +c
  355 + if (rr .gt. rfin .and. np .lt. 500 .and. rr .lt. rmax) go to 40
  356 +c
  357 + if (np .ge. 500 .or. rr .ge. rmax) then
  358 + ier = 1
  359 + tr(np) = 999.0d0
  360 + tthet(np) = 999.0d0 * rad
  361 + tphi(np) = 999.0d0 * rad
  362 + endif
  363 +c
  364 +c ----------------------------------------------------------
  365 +c*FON La limite rfin est depassee, on interpole pour obtenir les
  366 +c*FON coordonnees exactes du point en cartesiennes converties
  367 +c*FON en coordonnees spheriques apres calcul
  368 +c ----------------------------------------------------------
  369 +c
  370 + if (rr .lt. rfin) then
  371 +c
  372 + rpp = tr(np-2)
  373 + rp = tr(np-1)
  374 + r = tr(np)
  375 + xxpp = tx(np-2)
  376 + xxp = tx(np-1)
  377 + xx = tx(np)
  378 + yypp = ty(np-2)
  379 + yyp = ty(np-1)
  380 + yy = ty(np)
  381 + zzpp = tz(np-2)
  382 + zzp = tz(np-1)
  383 + zz = tz(np)
  384 + px = pinter(rpp,rp,r,xxpp,xxp,xx,rfin)
  385 + py = pinter(rpp,rp,r,yypp,yyp,yy,rfin)
  386 + pz = pinter(rpp,rp,r,zzpp,zzp,zz,rfin)
  387 +c
  388 + call carsp (px,py,pz,rr,thetr,phir,ier7)
  389 +c
  390 + tr(np) = rr
  391 + tthet(np) = thetr
  392 + tphi(np) = phir
  393 +c
  394 + endif
  395 +c
  396 +c --------------------------------------------------------------
  397 +c*FON Si la distance limite rmax ou si le nbr de points depasse 500
  398 +c*FON ou si on arrive en dehors de la magnetosphere de Shabansky on
  399 +c*FON remplit les coordonnes par des 999 ou 999*rad pour obtenir des
  400 +c*FON 999 en degres apres la conversion radians degres
  401 +c --------------------------------------------------------------
  402 +c
  403 +10 continue
  404 +c
  405 + if (ipop .eq. 2) then
  406 +c
  407 + tr(np) = 999.0d0
  408 + tthet(np) = 999.0d0 * rad
  409 + tphi(np) = 999.0d0 * rad
  410 + ier = 2
  411 +c
  412 + endif
  413 +30 continue
  414 +c
  415 +c ****************
  416 +c Fin de programme
  417 +c ****************
  418 +c
  419 + return
  420 + end
... ...
maglib/src/derivv.f 0 → 100755
  Diff comments   View file @510e044
... ... @@ -0,0 +1,228 @@
  1 + subroutine derivv (magin,magout,indgm,indval,year,tilt,
  2 + > rggsm,rgsmg,ds3,xgc,ygc,zgc,rxp1,rxp2,
  3 + > rxp3,ier)
  4 +c*
  5 +c***********************************************************************
  6 +c*
  7 +c* "Copyright [c] CNES 98 - tous droits reserves"
  8 +c* **********************************************
  9 +c*
  10 +c*pro MAGLIB
  11 +c*
  12 +c*VER 99.03.31 - V 1.0
  13 +c*VER 01.05.30 - V 2.0
  14 +c*VER 03.01.06 - V 2.1
  15 +c*VER 01.01.07 - V 3.0
  16 +c*VER 13.10.10 - V 4.0
  17 +c*VER 17.02.23 - V 5.0
  18 +c*VER 20.07.15 - V 6.0
  19 +c*
  20 +c*AUT spec. CNES - JC KOSIK - avril 1996
  21 +c*AUT port. CISI
  22 +c*AUT adpat. AKKA
  23 +c*
  24 +c*ROL Theme : Calculs de geophysique
  25 +c*ROL Calcul des pas elementaires rxp1,rxp2,rxp3 au point
  26 +c*ROL courant xgc,ygc,zgc.
  27 +c*
  28 +c*PAR magin (I) : type de champ magnetique interne
  29 +c*PAR magout (I) : type de champ magnetique externe
  30 +c*
  31 +c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae
  32 +c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ
  33 +c*
  34 +c*PAR year (I) : annee fractionnaire (>= 2000.)
  35 +c*
  36 +c*PAR tilt (I) : angle de tilt (radians)
  37 +c*
  38 +c*PAR rggsm (I) : matrice de passage du repere geographique au
  39 +c*PAR : repere solaire magnetospherique
  40 +c*PAR rgsmg (I) : matrice de passage du repere solaire
  41 +c*PAR : magnetospherique au repere geographique
  42 +c*
  43 +c*PAR ds3 (I) : pas elementaire d'integration
  44 +c*
  45 +c*PAR xgc (I) : coordonnee geocentrique en x
  46 +c*PAR ygc (I) : coordonnee geocentrique en y
  47 +c*PAR zgc (I) : coordonnee geocentrique en z
  48 +c*
  49 +c*PAR rxp1 (O) : pas elementaire sur x (au point xgc)
  50 +c*PAR rxp2 (O) : pas elementaire sur y (au point ygc)
  51 +c*PAR rxp3 (O) : pas elementaire sur z (au point zgc)
  52 +c*
  53 +c*PAR ier (O) : code de retour
  54 +c*
  55 +c*NOT magin : 1 = dipole + g11 + h11
  56 +c*NOT magin : 2 = champ DGRF 1945-1970
  57 +c*NOT magin : 3 = champ DGRF 1970-1995
  58 +c*
  59 +c*NOT magout : 0 = pas de champ externe
  60 +c*NOT magout : 1 = Tsyganenko 87
  61 +c*NOT magout : 2 = Tsyganenko 89
  62 +c*NOT magout : 3 = Kosik 97
  63 +c*
  64 +c*NOT indgm : 1 ---> indice geomagnetique Kp
  65 +c*NOT indval : 1 ---> Kp = 0 , 0+
  66 +c*NOT indval : 2 ---> Kp = 1- , 1 , 1+
  67 +c*NOT indval : 3 ---> Kp = 2- , 2 , 2+
  68 +c*NOT indval : 4 ---> Kp = 3- , 3 , 3+
  69 +c*NOT indval : 5 ---> Kp = 4- , 4 , 4+
  70 +c*NOT indval : 6 ---> Kp > 5-
  71 +c*
  72 +c*NOT indgm : 2 ---> indice geomagnetique Ae
  73 +c*NOT indval : 1 ---> Ae = 0 - 50
  74 +c*NOT indval : 2 ---> Ae = 50 - 100
  75 +c*NOT indval : 3 ---> Ae = 100 - 150
  76 +c*NOT indval : 4 ---> Ae = 150 - 250
  77 +c*NOT indval : 5 ---> Ae = 250 - 400
  78 +c*NOT indval : 6 ---> Ae >= 400
  79 +c*
  80 +c*NOT On utilise la formule generale dx/ds = bx/b
  81 +c*NOT ou ds est le pas elementaire le long de la ligne de champ, b le
  82 +c*NOT champ total, dx est la composante de ds suivant l'axe x et bx
  83 +c*NOT la composante du champ suivant ce meme axe.
  84 +c*
  85 +c*INF utilise : carsp, inmag, outma1, vspvcar
  86 +c*
  87 +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
  88 +c*HST version 2.0 - 01.05.30 - correction de commentaires de code
  89 +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
  90 +c*HST version 3.0 - 01.01.07 - mise a jour des modeles de champ interne
  91 +c*HST version 4.0 - 13.10.10 - mise a jour des modeles de champ interne
  92 +c*HST (DGRF05 et IGRF10)
  93 +c*HST version 5.0 - 17.02.23 - mise a jour des modeles de champ interne
  94 +c*HST (DGRF10 et IGRF15)
  95 +c*HST version 6.0 - 20.07.15 - mise a jour des modeles de champ interne
  96 +c*HST (DGRF15 et IGRF20)
  97 +c*
  98 +c***********************************************************************
  99 +c*
  100 + implicit none
  101 +c
  102 +c ---------------------------------
  103 +c*FON Declaration identificateur rcs_id
  104 +c ---------------------------------
  105 +c
  106 + character rcs_id*100
  107 +c
  108 +c --------------------------
  109 +c*FON Declaration des parametres
  110 +c --------------------------
  111 +c
  112 + integer indgm,indval,magin,magout
  113 + double precision year
  114 + double precision tilt
  115 + double precision rggsm(3,3), rgsmg(3,3)
  116 + double precision ds3
  117 + double precision xgc, ygc, zgc
  118 + double precision rxp1, rxp2, rxp3
  119 + integer ier
  120 +c
  121 +c ---------------------------------
  122 +c*FON Declaration des variables locales
  123 +c ---------------------------------
  124 +c
  125 + integer ier1,ier2,ier3,ier4
  126 +c*LOC ier1,ier2,ier3,ier4 : codes retour des modules appeles
  127 +c
  128 + double precision rre,thetr,phir
  129 +c*LOC rre,thetr,phir : coordonnees spheriques
  130 +c
  131 + double precision bre,bte,bpe
  132 +c*LOC bre,bte,bpe : composantes radiale, tangentielle et azimuthale du champ
  133 +c*LOC : d'origine externe
  134 +c
  135 + double precision bri,bti,bpi
  136 +c*LOC bri,bti,bpi : composantes radiale, tangentielle et azimuthale du champ
  137 +c*LOC : d'origine interne
  138 +c
  139 + double precision br,bt,bp
  140 +c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ total
  141 +c
  142 + double precision bi
  143 +c*LOC bi : module du champ
  144 +c
  145 + double precision bxt,byt,bzt
  146 +c*LOC bxt,byt,bzt : coordonnees solaires magnetospheriques x,y et z du champ
  147 +c
  148 + SAVE
  149 +c
  150 +c ---------------------------------
  151 +c*FON Affectation identificateur rcs_id
  152 +c ---------------------------------
  153 +c
  154 + data rcs_id /"
  155 + >$Id$"/
  156 +c
  157 +c ******************
  158 +c Debut de programme
  159 +c ******************
  160 +c
  161 +c -----------------------------------
  162 +c*FON Initialisations des codes de retour
  163 +c -----------------------------------
  164 +c
  165 + ier = 0
  166 + ier1 = 0
  167 + ier2 = 0
  168 + ier3 = 0
  169 + ier4 = 0
  170 +c
  171 +c ---------------------------------------------------
  172 +c*FON Initialisations des composantes du champ magnetique
  173 +c ---------------------------------------------------
  174 +c
  175 + bre = 0.d0
  176 + bte = 0.d0
  177 + bpe = 0.d0
  178 +c
  179 + bri = 0.d0
  180 + bti = 0.d0
  181 + bpi = 0.d0
  182 +c
  183 +c ----------------------------------------------------------------------
  184 +c*FON Transformation des coordonnees geocentriques en coordonnees spheriques
  185 +c ----------------------------------------------------------------------
  186 +c
  187 + call carsp (xgc,ygc,zgc,rre,thetr,phir,ier1)
  188 +c
  189 +c --------------------------------------------
  190 +c*FON Calcul du champ magnetique d'origine interne
  191 +c --------------------------------------------
  192 +c
  193 + call inmagv(magin,year,rre,thetr,phir,bri,bti,bpi,bi,ier2)
  194 +c
  195 +c --------------------------------------------
  196 +c*FON Calcul du champ magnetique d'origine externe
  197 +c --------------------------------------------
  198 +c
  199 + call outma1(magout,indgm,indval,tilt,rggsm,rgsmg,
  200 + > rre,thetr,phir,bre,bte,bpe,ier3)
  201 +c
  202 +c ---------------------------------------------------------------
  203 +c*FON Calcul des composantes du champ magnetique total en coordonnees
  204 +c*FON spheriques puis conversion en coordonnees cartesiennes.
  205 +c ---------------------------------------------------------------
  206 +c
  207 + br = bri + bre
  208 + bt = bti + bte
  209 + bp = bpi + bpe
  210 +c
  211 + call vspvcar(thetr,phir,br,bt,bp,bxt,byt,bzt,ier4)
  212 +c
  213 + bt = sqrt(bxt*bxt + byt*byt + bzt*bzt)
  214 +c
  215 +c -------------------------------------
  216 +c*FON Calcul des pas elementaires sur x,y,z
  217 +c -------------------------------------
  218 +c
  219 + rxp1 = ds3 * bxt / bt
  220 + rxp2 = ds3 * byt / bt
  221 + rxp3 = ds3 * bzt / bt
  222 +c
  223 +c ****************
  224 +c Fin de programme
  225 +c ****************
  226 +c
  227 + return
  228 + end
... ...
maglib/src/getmag.f
  Diff comments   View file @510e044
... ... @@ -100,21 +100,34 @@ c*FON by the geomagnetic calculations for the given year
100 100 c --------------------------------------------------
101 101  
102 102 c
103   - call inigeom(iyear,imonth,iday,ihour,imin,isec,year,alfag,
  103 + if (iyear .ge. 2000) then
  104 + call inigeom(iyear,imonth,iday,ihour,imin,isec,year,alfag,
104 105 > tetdip,phidip,alfas,deltas,rig,rgi,rgdip,rdipg,
105 106 > rgsm,rsmg,tilt,rggsm,rgsmg,rgse,rseg,rigsm,rgsmi,
106 107 > ifail)
107 108  
108   -c -------------------------------------------------------------
109   -c*FON Calculation of the position of the point in the magnetosphere
110   -c -------------------------------------------------------------
111   -c
112   -c
113 109 call posmag(magout,isatex,year,rrmag,thetr,phir,alfag,
114 110 > alfas,deltas,tilt,rgsm,rggsm,rgsmg,rgdip,rgse,
115 111 > tetdip,phidip,xgsm,ygsm,zgsm,xgse,ygse,zgse,tgl,
116 112 > flg,xlamb,tglc,hsl,clatgmr,clongmr,iposmg,ifail)
117 113  
  114 + else
  115 + call inigeomv(iyear,imonth,iday,ihour,imin,isec,year,alfag,
  116 + > tetdip,phidip,alfas,deltas,rig,rgi,rgdip,rdipg,
  117 + > rgsm,rsmg,tilt,rggsm,rgsmg,rgse,rseg,rigsm,rgsmi,
  118 + > ifail)
  119 +
  120 + call posmagv(magout,isatex,year,rrmag,thetr,phir,alfag,
  121 + > alfas,deltas,tilt,rgsm,rggsm,rgsmg,rgdip,rgse,
  122 + > tetdip,phidip,xgsm,ygsm,zgsm,xgse,ygse,zgse,tgl,
  123 + > flg,xlamb,tglc,hsl,clatgmr,clongmr,iposmg,ifail)
  124 +
  125 + endif
  126 +
  127 +c -------------------------------------------------------------
  128 +c*FON Calculation of the position of the point in the magnetosphere
  129 +c -------------------------------------------------------------
  130 +c
118 131  
119 132 c return
120 133  
... ...
maglib/src/inmagv.f 0 → 100755
  Diff comments   View file @510e044
... ... @@ -0,0 +1,144 @@
  1 + subroutine inmagv (magin,year,rre,thet,phi,bri,bti,bpi,bi,ier)
  2 +c*
  3 +c***********************************************************************
  4 +c*
  5 +c* "Copyright [c] CNES 98 - tous droits reserves"
  6 +c* **********************************************
  7 +c*
  8 +c*PRO MAGLIB
  9 +c*
  10 +c*VER 99.03.31 - V 1.0
  11 +c*VER 01.05.30 - V 2.0
  12 +c*VER 03.01.06 - V 2.1
  13 +c*VER 01.01.07 - V 3.0
  14 +c*VER 13.10.10 - V 4.0
  15 +c*VER 17.02.24 - V 5.0
  16 +c*VER 20.07.15 - V 6.0
  17 +c*
  18 +c*AUT spec. CNES - JC KOSIK - juin 1995
  19 +c*AUT port. CISI
  20 +c*AUT adapt. AKKA
  21 +c*
  22 +c*ROL Theme : Calculs de geophysique
  23 +c*ROL Calcul du champ magnetique d'origine interne.
  24 +c*
  25 +c*PAR magin (I) : type de champ magnetique interne
  26 +c*
  27 +c*PAR year (I) : annee fractionnaire (>= 1945. pour magin = 2 ;
  28 +c*PAR : >= 1970. pour magin = 3)
  29 +c*
  30 +c*PAR rre (I) : distance radiale geocentrique (rayons terrestres)
  31 +c*PAR thet (I) : colatitude geocentrique (radians)
  32 +c*PAR phi (I) : longitude geocentrique (radians)
  33 +c*
  34 +c*PAR bri (O) : composante radiale du champ magnetique le long du
  35 +c* : meridien positive vers l'exterieur (gauss)
  36 +c*PAR bti (O) : composante tangentielle du champ magnetique le long
  37 +c*PAR : du meridien positive vers le sud (gauss)
  38 +c*PAR bp (O) : composante azimuthale du champ magnetique, positive
  39 +c*PAR : vers l'est (gauss)
  40 +c*
  41 +c*PAR bi (O) : module du champ (gauss)
  42 +c*
  43 +c*PAR ier (O) : code de retour
  44 +c*
  45 +c*NOT magin : 1 = dipole + g11 + h11
  46 +c*NOT magin : 2 = champ DGRF 1945-1970
  47 +c*NOT magin : 3 = champ DGRF 1970-1995
  48 +c*
  49 +c*NOT ier : sans objet
  50 +c*
  51 +c*INF utilise : dipol, dgrf45_70, dgrf70_95
  52 +c*
  53 +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
  54 +c*HST version 2.0 - 01.05.30 - correction de commentaires de code
  55 +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
  56 +c*HST version 3.0 - 01.01.07 - utilisation du champ 2005
  57 +c*HST version 4.0 - 13.10.10 - Introduction du champ igrf10 et
  58 +c*HST passage de igrf05 a dgrf05
  59 +c*HST version 5.0 - 17.02.24 - Introduction du champ igrf15 et
  60 +c*HST passage de igrf10 a dgrf10
  61 +c*HST version 6.0 - 20.07.15 - Introduction du champ igrf20 et
  62 +c*HST passage de igrf15 a dgrf15
  63 +c*
  64 +c***********************************************************************
  65 +c*
  66 + implicit none
  67 +c
  68 +c ---------------------------------
  69 +c*FON Declaration identificateur rcs_id
  70 +c ---------------------------------
  71 +c
  72 + character rcs_id*100
  73 +c
  74 +c --------------------------
  75 +c*FON Declaration des parametres
  76 +c --------------------------
  77 +c
  78 + integer magin
  79 + double precision year
  80 + double precision rre, thet, phi
  81 + double precision bri, bti, bpi, bi
  82 + integer ier
  83 +c
  84 +c ---------------------------------
  85 +c*FON Declaration des variables locales
  86 +c ---------------------------------
  87 +c
  88 + integer ier1,ier2,ier3
  89 +c*LOC ier1,ier2,ier3 : codes retour des modules appeles
  90 +c
  91 + SAVE
  92 +c
  93 +c ---------------------------------
  94 +c*FON Affectation identificateur rcs_id
  95 +c ---------------------------------
  96 +c
  97 + data rcs_id /"
  98 + >$Id$"/
  99 +c
  100 +c ******************
  101 +c Debut de programme
  102 +c ******************
  103 +c
  104 + ier = 0
  105 + ier1 = 0
  106 + ier2 = 0
  107 + ier3 = 0
  108 +c
  109 +c ------------------------------------
  110 +c*FON Calcul du champ (dipole + g11 + h11)
  111 +c ------------------------------------
  112 +c
  113 + if (magin .eq. 1) then
  114 +c
  115 + call dipol(year,rre,thet,phi,bri,bti,bpi,bi,ier1)
  116 +c
  117 +c -------------------------
  118 +c*FON Calcul du champ DGRF 2000
  119 +c -------------------------
  120 +c
  121 + else if (magin .eq. 2) then
  122 +c
  123 + call dgrf45_70(year,rre,thet,phi,bri,bti,bpi,bi,ier2)
  124 +c
  125 +c -------------------------
  126 +c*FON Calcul du champ DGRF 2005
  127 +c -------------------------
  128 +c
  129 + else if (magin .eq. 3) then
  130 +c
  131 + call dgrf70_95(year,rre,thet,phi,bri,bti,bpi,bi,ier2)
  132 +c
  133 + else if(magin .eq. 4) then
  134 +c
  135 + call dgrf95(year,rre,thet,phi,bri,bti,bpi,bi,ier2)
  136 +c
  137 + endif
  138 +c
  139 +c ****************
  140 +c Fin de programme
  141 +c ****************
  142 +c
  143 + return
  144 + end
... ...
maglib/src/magtotv.f 0 → 100755
  Diff comments   View file @510e044
... ... @@ -0,0 +1,199 @@
  1 + subroutine magtotv (magin,year,magout,indgm,indval,tilt,rggsm,
  2 + > rgsmg,rre,thet,phi,br,bt,bp,bb,ier)
  3 +c*
  4 +c***********************************************************************
  5 +c*
  6 +c* "Copyright [c] CNES 98 - tous droits reserves"
  7 +c* **********************************************
  8 +c*
  9 +c*PRO MAGLIB
  10 +c*
  11 +c*VER 99.03.31 - V 1.0
  12 +c*VER 01.05.30 - V 2.0
  13 +c*VER 03.01.06 - V 2.1
  14 +c*VER 01.01.07 - V 3.0
  15 +c*VER 13.10.10 - V 4.0
  16 +c*VER 17.02.24 - V 5.0
  17 +c*VER 20.07.16 - V 6.0
  18 +c*
  19 +c*AUT spec. CNES - JC KOSIK - juin 1995
  20 +c*AUT port. CISI
  21 +c*AUT adapt. AKKA
  22 +c*
  23 +c*ROL Theme : Calculs de geophysique
  24 +c*ROL Calcul du champ magnetique total en un point.
  25 +c*ROL Ce champ est la combinaison possible d'un champ dipole
  26 +c*ROL ou complet et d'un champ externe, ici Tsyganenko.
  27 +c*
  28 +c*PAR magin (I) : type de champ magnetique interne
  29 +c*
  30 +c*PAR year (I) : annee fractionnaire (>= 1945.)
  31 +c*
  32 +c*PAR magout (I) : type de champ magnetique externe
  33 +c*
  34 +c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae
  35 +c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ
  36 +c*
  37 +c*PAR tilt (I) : angle de tilt (radians)
  38 +c*
  39 +c*PAR rggsm (I) : matrice (3,3) de passage du repere geocentrique
  40 +c*PAR : au repere magnetospherique
  41 +c*
  42 +c*PAR rgsmg (I) : matrice (3,3) de passage du repere solaire
  43 +c*PAR : magnetospherique au repere geocentrique
  44 +c*
  45 +c*PAR rre (I) : distance radiale geocentrique (rayons terrestres)
  46 +c*PAR thet (I) : colatitude geocentrique (radians)
  47 +c*PAR phi (I) : longitude geocentrique (radians)
  48 +c*
  49 +c*PAR br (O) : composante radiale du champ total (gauss)
  50 +c*PAR bt (O) : composante tangentielle du champ total (gauss)
  51 +c*PAR bp (O) : composante azimuthale du champ total (gauss)
  52 +c*
  53 +c*PAR bb (O) : module du champ magnetique total (gauss)
  54 +c*PAR ier (O) : code de retour
  55 +c*
  56 +c*NOT magin : 1 = dipole + g11 + h11
  57 +c*NOT magin : 2 = champ DGRF 1945-1970
  58 +c*NOT magin : 3 = champ DGRF 1970-1995
  59 +c*
  60 +c*NOT magout : 0 = pas de champ externe
  61 +c*NOT magout : 1 = Tsyganenko 87
  62 +c*NOT magout : 2 = Tsyganenko 89
  63 +c*NOT magout : 3 = Kosik 97
  64 +c*
  65 +c*NOT indgm : 1 ---> indice geomagnetique Kp
  66 +c*NOT indval : 1 ---> Kp = 0 , 0+
  67 +c*NOT indval : 2 ---> Kp = 1- , 1 , 1+
  68 +c*NOT indval : 3 ---> Kp = 2- , 2 , 2+
  69 +c*NOT indval : 4 ---> Kp = 3- , 3 , 3+
  70 +c*NOT indval : 5 ---> Kp = 4- , 4 , 4+
  71 +c*NOT indval : 6 ---> Kp > 5-
  72 +c*
  73 +c*NOT indgm : 2 ---> indice geomagnetique Ae
  74 +c*NOT indval : 1 ---> Ae = 0 - 50
  75 +c*NOT indval : 2 ---> Ae = 50 - 100
  76 +c*NOT indval : 3 ---> Ae = 100 - 150
  77 +c*NOT indval : 4 ---> Ae = 150 - 250
  78 +c*NOT indval : 5 ---> Ae = 250 - 400
  79 +c*NOT indval : 6 ---> Ae >= 400
  80 +c*
  81 +c*NOT ier : sans objet
  82 +c*
  83 +c*INF utilise : inmag, outma1
  84 +c*
  85 +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
  86 +c*HST version 2.0 - 01.05.30 - Ajout indice Kosik 97
  87 +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
  88 +c*HST version 3.0 - 01.01.07 - utilisation du champ 2005
  89 +c*HST version 4.0 - 13.10.10 - Introductiondu champ IGRF10 et
  90 +c*HST passage de IGRF05 a DGRF05
  91 +c*HST version 5.0 - 17.02.24 - Introduction du champ IGRF15 et
  92 +c*HST passage de IGRF10 a DGRF10
  93 +c*HST version 5.0 - 20.07.16 - Introduction du champ IGRF20 et
  94 +c*HST passage de IGRF15 a DGRF15
  95 +*
  96 +c***********************************************************************
  97 +c*
  98 + implicit none
  99 +c
  100 +c ---------------------------------
  101 +c*FON Declaration identificateur rcs_id
  102 +c ---------------------------------
  103 +c
  104 + character rcs_id*100
  105 +c
  106 +c --------------------------
  107 +c*FON Declaration des parametres
  108 +c --------------------------
  109 +c
  110 + integer magin
  111 + double precision year
  112 + integer magout, indgm, indval
  113 + double precision tilt
  114 + double precision rggsm(3,3), rgsmg(3,3)
  115 + double precision rre, thet, phi
  116 + double precision br, bt, bp, bb
  117 + integer ier
  118 +c
  119 +c ---------------------------------
  120 +c*FON Declaration des variables locales
  121 +c ---------------------------------
  122 +c
  123 + integer ier1,ier2
  124 +c*LOC ier1,ier2 : codes retour des modules appeles
  125 +c
  126 + double precision bpe,bre,bte
  127 +c*LOC bpe,bre,bte : composantes radiale, tangentielle et azimuthale du champ
  128 +c* : externe
  129 +c
  130 + double precision bpi,bri,bti
  131 +c*LOC bpi,bri,bti : composantes radiale, tangentielle et azimuthale du champ
  132 +c* : interne
  133 +c
  134 + double precision bi
  135 +c*LOC bi : module du champ
  136 +c
  137 + SAVE
  138 +c
  139 +c ---------------------------------
  140 +c*FON Affectation identificateur rcs_id
  141 +c ---------------------------------
  142 +c
  143 + data rcs_id /"
  144 + >$Id$"/
  145 +c
  146 +c ******************
  147 +c Debut de programme
  148 +c ******************
  149 +c
  150 + ier = 0
  151 + ier1 = 0
  152 + ier2 = 0
  153 +c
  154 +c --------------------------------------------
  155 +c*FON Initialisation a 0. des composantes du champ
  156 +c --------------------------------------------
  157 +c
  158 + bre = 0.d0
  159 + bte = 0.d0
  160 + bpe = 0.d0
  161 +c
  162 + bri = 0.d0
  163 + bti = 0.d0
  164 + bpi = 0.d0
  165 +c
  166 + if (magin .ne. 0) then
  167 +c
  168 +c --------------------------------------------
  169 +c*FON Calcul du champ magnetique d'origine interne
  170 +c --------------------------------------------
  171 +c
  172 + call inmagv(magin,year,rre,thet,phi,bri,bti,bpi,bi,ier1)
  173 +c
  174 + endif
  175 +c
  176 +c --------------------------------------------
  177 +c*FON Calcul du champ magnetique d'origine externe
  178 +c --------------------------------------------
  179 +c
  180 + if (magout .ne. 0) then
  181 + call outma1(magout,indgm,indval,tilt,rggsm,rgsmg,
  182 + > rre,thet,phi,bre,bte,bpe,ier2)
  183 + endif
  184 +c
  185 +c --------------------------------
  186 +c*FON Calcul du champ magnetique total
  187 +c --------------------------------
  188 +c
  189 + br = bri + bre
  190 + bt = bti + bte
  191 + bp = bpi + bpe
  192 + bb = sqrt(br * br + bt * bt + bp * bp)
  193 +c
  194 +c ****************
  195 +c Fin de programme
  196 +c ****************
  197 +c
  198 + return
  199 + end
... ...
maglib/src/posmag.f
  Diff comments   View file @510e044
... ... @@ -30,7 +30,7 @@ c*PAR magout (I) : type de champ magnetique externe
30 30 c*
31 31 c*PAR isatex (I) : type de satellite
32 32 c*
33   -c*PAR year (I) : annee fractionnaire >= 2005
  33 +c*PAR year (I) : annee fractionnaire >= 2000
34 34 c*
35 35 c*PAR rre (I) : distance radiale geocentrique (rayons terrestres)
36 36 c*PAR thetr (I) : colatitude geocentrique (radians)
... ... @@ -325,10 +325,39 @@ c
325 325 c
326 326 c ------------------------------------------------------
327 327 c*FON Calcul du point conjugue situe dans le meme hemisphere
328   -c*FON magin = 3 = champ IGRF 2005
  328 +c*NOT magin : 2 = champ DGRF 2000
  329 +c*NOT magin : 3 = champ DGRF 2005
  330 +c*NOT magin : 4 = champ DGRF 2010
  331 +c*NOT magin : 5 = champ DGRF 2015
  332 +c*NOT magin : 6 = champ IGRF 2020
  333 +
329 334 c ------------------------------------------------------
330 335 c
331   - magin = 3
  336 + if (year .ge. 2000 .and. year .lt. 2005) then
  337 +c
  338 + magin = 2
  339 +c
  340 + else if (year .ge. 2005 .and. year .lt. 2010) then
  341 +c
  342 + magin = 3
  343 +c
  344 + else if (year .ge. 2010 .and. year .lt. 2015) then
  345 +c
  346 + magin = 4
  347 +c
  348 + else if (year .ge. 2015 .and. year .lt. 2020) then
  349 +c
  350 + magin = 5
  351 +c
  352 + else if (year .ge. 2020) then
  353 +c
  354 + magin = 6
  355 +c
  356 + else
  357 + ier = 1
  358 + endif
  359 +c
  360 +c
332 361 indgm = 1
333 362 indval = 3
334 363 rb = 11.d0
... ...
maglib/src/posmagv.f 0 → 100755
  Diff comments   View file @510e044
... ... @@ -0,0 +1,481 @@
  1 + subroutine posmagv (magout,isatex,year,rre,thetr,phir,alfag,alfas,
  2 + > deltas,tilt,rgsm,rggsm,rgsmg,rgdip,rgse,
  3 + > tetdip,phidip,xgsm,ygsm,zgsm,xgse,ygse,zgse,
  4 + > tgl,flg,xlamb,tglc,hsl,clatgm,clongm,
  5 + > iposmg,ier)
  6 +c*
  7 +c***********************************************************************
  8 +c*
  9 +c* "Copyright [c] CNES 98 - tous droits reserves"
  10 +c* **********************************************
  11 +c*
  12 +c*PRO MAGLIB
  13 +c*
  14 +c*VER 99.03.31 - V 1.0
  15 +c*VER 01.05.30 - V 2.0
  16 +c*VER 03.01.06 - V 2.1
  17 +c*VER 01.01.07 - V 3.0
  18 +c*
  19 +c*AUT spec. CNES - JC KOSIK - janvier 1991
  20 +c*AUT port. CISI
  21 +c*
  22 +c*ROL Theme : Frontieres et regions
  23 +c*ROL Calcul de la position d'un satellite dans la magnetosphere
  24 +c*ROL par rapport aux differentres regions.
  25 +c*ROL Calcul aussi des parametres essentiels tgl, L et mlt,
  26 +c*ROL et des coordonnees dans divers systemes d axes.
  27 +c*ROL L'indicateur de positionnement iposmg est renseigne.
  28 +c*
  29 +c*PAR magout (I) : type de champ magnetique externe
  30 +c*
  31 +c*PAR isatex (I) : type de satellite
  32 +c*
  33 +c*PAR year (I) : annee fractionnaire >= 2005
  34 +c*
  35 +c*PAR rre (I) : distance radiale geocentrique (rayons terrestres)
  36 +c*PAR thetr (I) : colatitude geocentrique (radians)
  37 +c*PAR phir (I) : longitude geocentrique (radians)
  38 +c*
  39 +c*PAR alfag (I) : ascension droite de Greenwich (radians)
  40 +c*PAR alfas (I) : ascension droite du soleil (radians)
  41 +c*PAR deltas (I) : declinaison du soleil (radians)
  42 +c*
  43 +c*PAR tilt (I) : angle de tilt (radians)
  44 +c*
  45 +c*PAR rgsm (I) : matrice de passage du repere geographique au
  46 +c*PAR : repere solaire magnetique
  47 +c*
  48 +c*PAR rggsm (I) : matrice de passage du repere geographique au
  49 +c*PAR : repere solaire magnetospherique
  50 +c*PAR rgsmg (I) : matrice de passage du repere magnetospherique
  51 +c*PAR : au repere geographique
  52 +c*
  53 +c*PAR rgdip (I) : matrice (3,3) de passage du repere geographique
  54 +c*PAR : au repere dipolaire
  55 +c*
  56 +c*PAR rgse (I) : matrice de passage du repere geographique
  57 +c*PAR : au repere solaire ecliptique
  58 +c*
  59 +c*PAR tetdip (I) : colatitude geocentrique du dipole (radians)
  60 +c*PAR phidip (I) : longitude geocentrique du dipole (radians)
  61 +c*
  62 +c*PAR xgsm (O) : coordonnee solaire magnetospherique en x (rayons terrestres)
  63 +c*PAR ygsm (O) : coordonnee solaire magnetospherique en y (rayons terrestres)
  64 +c*PAR zgsm (O) : coordonnee solaire magnetospherique en z (rayons terrestres)
  65 +c*
  66 +c*PAR xgse (O) : coordonnee solaire ecliptique en x (rayons terrestres)
  67 +c*PAR ygse (O) : coordonnee solaire ecliptique en y (rayons terrestres)
  68 +c*PAR zgse (O) : coordonnee solaire ecliptique en z (rayons terrestres)
  69 +c*
  70 +c*PAR tgl (O) : temps geomagnetique local du satellite
  71 +c* : (heures fractionnaires)
  72 +c*PAR flg (O) : parametre L de Galperin
  73 +c*PAR xlamb (O) : latitude invariante (radians)
  74 +c*
  75 +c*PAR tglc (O) : temps geomagnetique local du point conjugue
  76 +c* : (heures fractionnaires)
  77 +c*
  78 +c*PAR hsl (O) : hauteur de l'ombre (kilometres)
  79 +c*
  80 +c*PAR clatgm (O) : latitude geomagnetique (radians)
  81 +c*PAR clongm (O) : longitude geomagnetique (radians)
  82 +c*PAR iposmg (O) : tableau (dim=15) des indicateurs de positions
  83 +c*
  84 +c*PAR ier (O) : code de retour
  85 +c*
  86 +c*NOT magout : 0 = pas de champ externe
  87 +c*NOT magout : 1 = Tsyganenko 87
  88 +c*NOT magout : 2 = Tsyganenko 89
  89 +c*NOT magout : 3 = kosik 97
  90 +c*
  91 +c*NOT isatex : 0 = satellite proche (apogee < 8. r. terre.)
  92 +c*NOT isatex : 1 = satellite excentrique (apogee > 8. r. terre.)
  93 +c*
  94 +c*NOT iposmg : 1 = le satellite appartient a la region
  95 +c*NOT iposmg : 0 = le satellite n'appartient pas a la region
  96 +c*
  97 +c*NOT iposmg(1) : calotte polaire
  98 +c*NOT iposmg(2) : ovale auroral
  99 +c*NOT iposmg(3) : cusp
  100 +c*NOT iposmg(4) : zone auroral diffuse
  101 +c*NOT iposmg(5) : ceinture de Van Allen pour xbelt = 6.
  102 +c*NOT iposmg(6) : plasmasphere
  103 +c*NOT iposmg(7) : couche neutre
  104 +c*NOT iposmg(8) : couche de plasma
  105 +c*NOT iposmg(9) : magnetosphere
  106 +c*NOT iposmg(10) : magnetopause
  107 +c*NOT iposmg(11) : magnetogaine
  108 +c*NOT iposmg(12) : onde de choc
  109 +c*NOT iposmg(13) : vent solaire
  110 +c*NOT iposmg(14) : ombre de la terre (satellite auroral)
  111 +c*NOT iposmg(15) : ceinture de Van Allen pour xbelt = 3.5
  112 +c*
  113 +c*NOT ier : sans objet
  114 +c*
  115 +c*NOT common : util
  116 +c*
  117 +c*INF utilise : spcar, geogsm, tgml, geose, bwshff, mpsib, msheath
  118 +c*INF utilise : dconjr, invar, invlat, ggeom, oval, calpol, cusp
  119 +c*INF utilise : gussen, rbelt, chapel, cahsl, posnsh, pospsh
  120 +c*
  121 +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
  122 +c*HST version 2.0 - 01.05.30 - correction de commentaires de code
  123 +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
  124 +c*HST version 3.0 - 01.01.07 - utilisation du champ IGRF 2005
  125 +c*
  126 +c***********************************************************************
  127 +c*
  128 + implicit none
  129 +c
  130 +c ---------------------------------
  131 +c*FON Declaration identificateur rcs_id
  132 +c ---------------------------------
  133 +c
  134 + character rcs_id*100
  135 +c
  136 +c --------------------------
  137 +c*FON Declaration des parametres
  138 +c --------------------------
  139 +c
  140 + integer magout, isatex
  141 + double precision year
  142 + double precision rre, thetr, phir
  143 + double precision alfag, alfas, deltas, tilt
  144 + double precision rgsm(3,3), rggsm(3,3), rgsmg(3,3)
  145 + double precision rgdip(3,3), rgse(3,3)
  146 + double precision tetdip, phidip
  147 + double precision tr(500), tthet(500), tphi(500)
  148 + double precision xgsm, ygsm, zgsm
  149 + double precision xgse, ygse, zgse
  150 + double precision tgl, flg, xlamb, tglc, hsl
  151 + double precision clatgm, clongm
  152 + integer iposmg(15)
  153 + double precision znsh
  154 + integer ier
  155 +c
  156 +c ----------------------------------
  157 +c*FON Declaration des variables communes
  158 +c ----------------------------------
  159 +c
  160 + double precision pi,dpi,rad,deg,pid,xmu,rayt
  161 +c
  162 +c*COM pi : constante pi (obtenue a partir de acos(-1.))
  163 +c*COM dpi : constante 2 * pi
  164 +c*COM pid : constante pi / 2
  165 +c*COM rad : facteur de conversion degres ----> radians
  166 +c*COM deg : facteur de conversion radians ----> degres
  167 +c*COM xmu : constante de gravitation terrestre (km**3/sec**2)
  168 +c*COM rayt : rayon equatorial terrestre (km)
  169 +c
  170 + common/util/pi,dpi,rad,deg,pid,xmu,rayt
  171 +c
  172 +c ---------------------------------
  173 +c*FON Declaration des variables locales
  174 +c ---------------------------------
  175 +c
  176 + integer imagn,ibwsh,imp,isheath,isolw,ioval,ipol
  177 + integer icusp,iguss,iallen,ichapp,ihsl,insh,ipsh
  178 +c*LOC Indicateurs d'appartenance aux zones
  179 +c
  180 + integer i
  181 +c*LOC i : indice de boucle
  182 +c
  183 + integer magin
  184 +c*LOC magin : type de champ magnetique interne (= 2 : DGRF45_70; 3 : DGRF70_95)
  185 +c
  186 + integer indgm
  187 +c*LOC indgm : type d'indice geomagnetique (=1 : Kp)
  188 +c
  189 + integer indval
  190 +c*LOC indval : indice geomagnetique : niveau d'amplitude du champ (1 a 6)
  191 +c
  192 + integer np
  193 +c*LOC np : nombre de points calcules
  194 +c
  195 + integer ier1,ier2,ier3,ier4,ier5,ier6,ier7,ier9
  196 + integer ier10,ier11,ier12,ier13,ier14,ier15,ier16
  197 + integer ier17,ier18,ier19
  198 +c*LOC ier1 a ier19 : codes retour des modules appeles
  199 +c
  200 + double precision rmagc,thetc,phic
  201 +c*LOC rmagc,thetc,phic : coordonnees spheriques
  202 +c*
  203 + double precision xmag,ymag,zmag
  204 +c*LOC xmag,ymag,zmag : composantes cartesiennes en x, y et z
  205 +c
  206 + double precision xbelt
  207 +c*LOC xbelt : valeur de la limite externe de la ceinture
  208 +c
  209 + double precision rfin,rmax
  210 +c*LOC rfin,rmax : distances geocentriques d'arret et maximum autorisee
  211 +c
  212 + double precision dir
  213 +c*LOC dir : direction du trace des lignes de champ
  214 +c
  215 + double precision rb
  216 +c*LOC rb : distance subsolaire de la magnetopause
  217 +c
  218 + SAVE
  219 +c
  220 +c ---------------------------------
  221 +c*FON Affectation identificateur rcs_id
  222 +c ---------------------------------
  223 +c
  224 + data rcs_id /"
  225 + >$Id$"/
  226 +c
  227 +c ******************
  228 +c Debut de programme
  229 +c ******************
  230 +c
  231 + ier = 0
  232 + ier1 = 0
  233 + ier2 = 0
  234 + ier3 = 0
  235 + ier4 = 0
  236 + ier5 = 0
  237 + ier6 = 0
  238 + ier7 = 0
  239 + ier9 = 0
  240 + ier10 = 0
  241 + ier11 = 0
  242 + ier12 = 0
  243 + ier13 = 0
  244 + ier14 = 0
  245 + ier15 = 0
  246 + ier16 = 0
  247 + ier17 = 0
  248 + ier18 = 0
  249 + ier19 = 0
  250 +c
  251 +c --------------------------------------------------
  252 +c*FON Mise a 0 du tableau d'indicateur de positionnement
  253 +c*FON Initialisations de valeurs par defaut
  254 +c --------------------------------------------------
  255 +c
  256 + do 10 i = 1, 15
  257 + iposmg(i) = 0
  258 +10 continue
  259 +c
  260 +c dans tous les cas mais surtout pour satellites auroraux :
  261 + iposmg(9) = 1
  262 +c
  263 + tgl = 999.0d0
  264 + flg = 999.0d0
  265 + xlamb = 999.0d0 * rad
  266 + tglc = 999.0d0
  267 + clatgm = 999.0d0 * rad
  268 + clongm = 999.0d0 * rad
  269 +c
  270 +c -------------------------------------------------
  271 +c*FON Calcul des coordonnees solaires magnetospheriques
  272 +c -------------------------------------------------
  273 +c
  274 + call spcar(rre,thetr,phir,xmag,ymag,zmag,ier1)
  275 +c
  276 + call geogsm(rggsm,xmag,ymag,zmag,xgsm,ygsm,zgsm,ier2)
  277 +c
  278 +c ---------------------------------------------------
  279 +c*FON Calcul du temps geomagnetique local du satellite tgl
  280 +c ---------------------------------------------------
  281 +c
  282 + call tgml(rgsm,thetr,phir,tgl,ier3)
  283 +c
  284 +c ------------------------------------------
  285 +c*FON Calcul des coordonnees solaire-ecliptiques
  286 +c ------------------------------------------
  287 +c
  288 + call geose(rgse,xmag,ymag,zmag,xgse,ygse,zgse,ier4)
  289 +c
  290 +c --------------------------------------------------------
  291 +c*FON Si nosat est superieur a 2, verification si le satellite
  292 +c*FON Appartient a la magnetosphere de sibeck
  293 +c --------------------------------------------------------
  294 +c
  295 + imagn = 1
  296 +c
  297 + if (isatex .eq. 1) then
  298 +c
  299 +c ---------------------------------------------------------------
  300 +c*FON Calcul de la position par rapport a l onde de choc de Fairfield
  301 +c*FON et la magnetopause de Sibeck
  302 +c ---------------------------------------------------------------
  303 +c
  304 + call bwshff(xgsm,ygsm,zgsm,ibwsh,ier5)
  305 +c
  306 + if (ibwsh .eq. 1) then
  307 + iposmg(12) = ibwsh
  308 + endif
  309 +c
  310 + call mpsib(xgsm,ygsm,zgsm,imp,ier6)
  311 +c
  312 + if (imp .eq. 1) then
  313 + iposmg(10) = imp
  314 + endif
  315 +c
  316 + call msheath (xgsm,ygsm,zgsm,imagn,isheath,isolw,ier7)
  317 +c
  318 + iposmg(13) = isolw
  319 + iposmg(11) = isheath
  320 + iposmg(9) = imagn
  321 +c
  322 + endif
  323 +c
  324 + if (imagn .eq. 1) then
  325 +c
  326 +c ------------------------------------------------------
  327 +c*FON Calcul du point conjugue situe dans le meme hemisphere
  328 +c*NOT magin : 1 = dipole + g11 + h11
  329 +c*NOT magin : 2 = champ DGRF 1945-1970
  330 +c*NOT magin : 3 = champ DGRF 1970-1995
  331 +c*NOT magin : 4 = champ DGRF 1995-2000
  332 +
  333 +c ------------------------------------------------------
  334 +c
  335 +c
  336 + if (year .ge. 1945 .and. year .lt. 1970) then
  337 +c
  338 + magin = 2
  339 +c
  340 + else if (year .ge. 1970 .and. year .lt. 1995) then
  341 +c
  342 + magin = 3
  343 +c
  344 + else if (year .ge. 1995 .and. year .lt. 2000) then
  345 +c
  346 + magin = 4
  347 +c
  348 + else
  349 + ier = 1
  350 + endif
  351 +
  352 +
  353 + indgm = 1
  354 + indval = 3
  355 + rb = 11.d0
  356 + rfin = 1.d0
  357 + rmax = 40.d0
  358 + dir = -1.d0
  359 +c
  360 + call dconjrv(magin,year,magout,indgm,indval,tilt,rb,rggsm,
  361 + > rgsmg,dir,rre,thetr,phir,rfin,rmax,np,tr,
  362 + > tthet,tphi,ier1)
  363 +c
  364 + if (ier1 .eq. 0) then
  365 +c
  366 + rmagc = tr(np)
  367 + thetc = tthet(np)
  368 + phic = tphi(np)
  369 +c
  370 + call invar(rmagc,thetc,phic,flg,ier2)
  371 +c
  372 + call invlat(flg,xlamb,ier3)
  373 +c
  374 +c ---------------------------------------
  375 +c*FON Calcul d appartenance a l ovale auroral
  376 +c ---------------------------------------
  377 +c
  378 +c ----------------------------------------------------------
  379 +c*FON Calcul des coordonnees geomagnetiques et mlt (appele tglc)
  380 +c*FON du pied de la ligne de champ passant par le satellite
  381 +c ----------------------------------------------------------
  382 +c
  383 + call ggeom(rgdip,thetc,phic,clatgm,clongm,ier9)
  384 + call tgml(rgsm,thetc,phic,tglc,ier10)
  385 +
  386 + call oval(tglc,clatgm,ioval,ier11)
  387 + call calpol(tglc,clatgm,ipol,ier12)
  388 +c
  389 +c -------------------------------------------
  390 +c*FON Calcul de la position par rapport a la CUSP
  391 +c -------------------------------------------
  392 +c
  393 + call cusp(clatgm,tglc,icusp,ier13)
  394 +c
  395 +c ------------------------------------------------------
  396 +c*FON Calcul de la frontiere sud de la zone aurorale diffuse
  397 +c ------------------------------------------------------
  398 +c
  399 + call gussen(tglc,xlamb,iguss,ier14)
  400 +c
  401 + iposmg(1) = ipol
  402 + iposmg(2) = ioval
  403 + iposmg(3) = icusp
  404 + iposmg(4) = iguss
  405 +c
  406 + endif
  407 +c
  408 +c -------------------------------------------------------
  409 +c*FON Si le satellite est a moins de 7 rayons terrestres.....
  410 +c -------------------------------------------------------
  411 +c
  412 + if (rre .lt. 7.d0) then
  413 +c
  414 +c ------------------------------------------------
  415 +c*FON Calcul d appartenance a la ceinture de Van Allen
  416 +c ------------------------------------------------
  417 +c
  418 +c Ceinture a 3.5 r. terr.
  419 +c -----------------------
  420 +c
  421 + xbelt = 3.5d0
  422 + call rbelt(xbelt,tetdip,phidip,rre,thetr,
  423 + + phir,iallen,ier15)
  424 + iposmg(15) = iallen
  425 +c
  426 +c Ceinture a 6.0 r. terr.
  427 +c -----------------------
  428 +c
  429 + xbelt = 6.d0
  430 + call rbelt(xbelt,tetdip,phidip,rre,thetr,
  431 + + phir,iallen,ier15)
  432 + iposmg(5) = iallen
  433 +c
  434 +c ---------------------------------------------------
  435 +c*FON Calcul d appartenance a la plasmasphere de Chappell
  436 +c ---------------------------------------------------
  437 +c
  438 + call chapel(tgl,flg,ichapp,ier16)
  439 + iposmg(6) = ichapp
  440 +c
  441 +c --------------------------------------------------------
  442 +c*FON Calcul de la distance du satellite a l ombre de la terre
  443 +c --------------------------------------------------------
  444 +c
  445 + call cahsl(alfag,alfas,deltas,rre,thetr,phir,
  446 + + hsl,ihsl,ier17)
  447 + iposmg(14) = ihsl
  448 +c
  449 + endif
  450 +c
  451 + endif
  452 +c
  453 +c -----------------------------------------------------------------
  454 +c*FON Calcul de la distance a la couche neutre et a la couche de plasma
  455 +c*FON si satellite excentrique et xgsm inferieur a -9 rt
  456 +c -----------------------------------------------------------------
  457 +c
  458 + if (isatex .eq. 1 .and. xgsm .le. -9.d0) then
  459 +c
  460 +c ----------------------------------------------------
  461 +c*FON Calcul de la position par rapport a la couche neutre
  462 +c ----------------------------------------------------
  463 +c
  464 + call posnsh(tilt,xgsm,ygsm,zgsm,insh,znsh,ier18)
  465 + iposmg(7) = insh
  466 +c
  467 +c -------------------------------------------------------
  468 +c*FON Calcul de la position par rapport a la couche de plasma
  469 +c -------------------------------------------------------
  470 +c
  471 + call pospsh(tilt,xgsm,ygsm,zgsm,ipsh,ier19)
  472 + iposmg(8) = ipsh
  473 +c
  474 + endif
  475 +c
  476 +c ****************
  477 +c Fin de programme
  478 +c ****************
  479 +c
  480 + return
  481 + end
... ...
maglib/src/stepv.f 0 → 100755
  Diff comments   View file @510e044
... ... @@ -0,0 +1,236 @@
  1 + subroutine stepv (magin,magout,indgm,indval,year,tilt,
  2 + > rggsm,rgsmg,x,y,z,ds,errin,ier)
  3 +c*
  4 +c***********************************************************************
  5 +c*
  6 +c* "Copyright [c] CNES 98 - tous droits reserves"
  7 +c* **********************************************
  8 +c*
  9 +c*PRO MAGLIB
  10 +c*
  11 +c*VER 99.03.31 - V 1.0
  12 +c*VER 01.05.30 - V 2.0
  13 +c*VER 03.01.06 - V 2.1
  14 +c*VER 01.01.07 - V 3.0
  15 +c*
  16 +c*AUT spec. CNES - JC KOSIK - avril 1996
  17 +c*AUT port. CISI
  18 +c*
  19 +c*ROL Theme : Calculs de geophysique
  20 +c*ROL Calcul d'un pas elementaire de la ligne de champ par la
  21 +c*ROL methode d'integration de Runge Kutta a pas varaiable et
  22 +c*ROL les coefficients de Merson.
  23 +c*
  24 +c*PAR magin (I) : type de champ magnetique interne
  25 +c*PAR magout (I) : type de champ magnetique externe
  26 +c*
  27 +c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae
  28 +c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ
  29 +c*
  30 +c*PAR year (I) : annee fractionnaire (>= 2000.)
  31 +c*
  32 +c*PAR tilt (I) : angle de tilt (radians)
  33 +c*
  34 +c*PAR rggsm (I) : matrice de passage du repere geographique au
  35 +c*PAR : repere solaire magnetospherique
  36 +c*PAR rgsmg (I) : matrice de passage du repere solaire
  37 +c*PAR : magnetospherique au repere geographique
  38 +c*
  39 +c*PAR x (I/O) : coordonnee geocentrique en x
  40 +c*PAR y (I/O) : coordonnee geocentrique en y
  41 +c*PAR z (I/O) : coordonnee geocentrique en z
  42 +c*
  43 +c*PAR ds (I) : pas elementaire d'integration
  44 +c*
  45 +c*PAR errin (I) : erreur maximale permise lors de l'integration
  46 +c*
  47 +c*PAR ier (O) : code retour
  48 +c*
  49 +c*NOT magin : 1 = dipole + g11 + h11
  50 +c*NOT magin : 2 = champ DGRF 1945-1970
  51 +c*NOT magin : 3 = champ DGRF 1970-1995
  52 +c*
  53 +c*NOT magout : 0 = pas de champ externe
  54 +c*NOT magout : 1 = Tsyganenko 87
  55 +c*NOT magout : 2 = Tsyganenko 89
  56 +c*NOT magout : 3 = Kosik 97
  57 +c*
  58 +c*NOT indgm : 1 ---> indice geomagnetique Kp
  59 +c*NOT indval : 1 ---> Kp = 0 , 0+
  60 +c*NOT indval : 2 ---> Kp = 1- , 1 , 1+
  61 +c*NOT indval : 3 ---> Kp = 2- , 2 , 2+
  62 +c*NOT indval : 4 ---> Kp = 3- , 3 , 3+
  63 +c*NOT indval : 5 ---> Kp = 4- , 4 , 4+
  64 +c*NOT indval : 6 ---> Kp > 5-
  65 +c*
  66 +c*NOT indgm : 2 ---> indice geomagnetique Ae
  67 +c*NOT indval : 1 ---> Ae = 0 - 50
  68 +c*NOT indval : 2 ---> Ae = 50 - 100
  69 +c*NOT indval : 3 ---> Ae = 100 - 150
  70 +c*NOT indval : 4 ---> Ae = 150 - 250
  71 +c*NOT indval : 5 ---> Ae = 250 - 400
  72 +c*NOT indval : 6 ---> Ae >= 400
  73 +c*
  74 +c*NOT La division par 3 du pas, permet des coefficients plus simples.
  75 +c*
  76 +c*INF utilise : deriv
  77 +c*
  78 +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
  79 +c*HST version 2.0 - 01.05.30 - correction de commentaires de code
  80 +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
  81 +c*HST version 3.0 - 01.01.07 - extension au champ IGRF 2005
  82 +c*
  83 +c***********************************************************************
  84 +c*
  85 + implicit none
  86 +c
  87 +c ---------------------------------
  88 +c*FON Declaration identificateur rcs_id
  89 +c ---------------------------------
  90 +c
  91 + character rcs_id*100
  92 +c
  93 +c --------------------------
  94 +c*FON Declaration des parametres
  95 +c --------------------------
  96 +c
  97 + integer magin, magout, indgm, indval
  98 + double precision year
  99 + double precision tilt
  100 + double precision rggsm(3,3), rgsmg(3,3)
  101 + double precision x, y, z
  102 + double precision ds
  103 + double precision errin
  104 + integer ier
  105 +c
  106 +c ---------------------------------
  107 +c*FON Declaration des variables locales
  108 +c ---------------------------------
  109 +c
  110 + double precision r11,r12,r13
  111 +c*LOC r11,r12,r13 : pas elementaires sur x, y et z (aux points x, y et z)
  112 +c
  113 + double precision x2,y2,z2
  114 +c*LOC x2,y2,z2 : coordonnees cartesiennes du point au 2eme pas
  115 +c
  116 + double precision r21,r22,r23
  117 +c*LOC r21,r22,r23 : pas elementaires sur x, y et z (aux points x2, y2 et z2)
  118 +c
  119 + double precision x3,y3,z3
  120 +c*LOC x3,y3,z3 : coordonnees cartesiennes du point au 3eme pas
  121 +c
  122 + double precision r31,r32,r33
  123 +c*LOC r31,r32,r33 : pas elementaires sur x, y et z (aux points x3, y3 et z3)
  124 +c
  125 + double precision x4,y4,z4
  126 +c*LOC x4,y4,z4 : coordonnees cartesiennes du point au 3eme pas
  127 +c
  128 + double precision r41,r42,r43
  129 +c*LOC r41,r42,r43 : pas elementaires sur x, y et z (aux points x4, y4 et z4)
  130 +c
  131 + double precision x5,y5,z5
  132 +c*LOC x5,y5,z5 : coordonnees cartesiennes du point au 5eme pas
  133 +c
  134 + double precision r51,r52,r53
  135 +c*LOC r51,r52,r53 : pas elementaires sur x, y et z (aux points x5, y5 et z5)
  136 +c
  137 + double precision ds3,errcur
  138 +c*LOC Variables de travail intermediaires
  139 +c
  140 + SAVE
  141 +c
  142 +c ---------------------------------
  143 +c*FON Affectation identificateur rcs_id
  144 +c ---------------------------------
  145 +c
  146 + data rcs_id /"
  147 + >$Id$"/
  148 +c
  149 +c ******************
  150 +c Debut de programme
  151 +c ******************
  152 +c
  153 + ier = 0
  154 +c
  155 +10 continue
  156 +c
  157 + ds3 = ds / 3.d0
  158 +c
  159 +c ----------------------------------------------------------
  160 +c*FON 1er appel au sous programme de calcul des pas elementaires
  161 +c ----------------------------------------------------------
  162 +c
  163 + call derivv(magin,magout,indgm,indval,year,tilt,
  164 + > rggsm,rgsmg,ds3,x,y,z,r11,r12,r13,ier)
  165 +c
  166 + x2 = x + r11
  167 + y2 = y + r12
  168 + z2 = z + r13
  169 +c
  170 +c ------------------------------------------------------------
  171 +c*FON 2ieme appel au sous programme de calcul des pas elementaires
  172 +c ------------------------------------------------------------
  173 +c
  174 + call derivv(magin,magout,indgm,indval,year,tilt,
  175 + > rggsm,rgsmg,ds3,x2,y2,z2,r21,r22,r23,ier)
  176 +c
  177 + x3 = x + .5d0 * (r11 + r21)
  178 + y3 = y + .5d0 * (r12 + r22)
  179 + z3 = z + .5d0 * (r13 + r23)
  180 +c
  181 +c ------------------------------------------------------------
  182 +c*FON 3ieme appel au sous programme de calcul des pas elementaires
  183 +c ------------------------------------------------------------
  184 +c
  185 + call derivv(magin,magout,indgm,indval,year,tilt,
  186 + > rggsm,rgsmg,ds3,x3,y3,z3,r31,r32,r33,ier)
  187 +c
  188 + x4 = x +.375d0 * (r11 + 3.d0 * r31)
  189 + y4 = y +.375d0 * (r12 + 3.d0 * r32)
  190 + z4 = z +.375d0 * (r13 + 3.d0 * r33)
  191 +c
  192 +c ------------------------------------------------------------
  193 +c*FON 4ieme appel au sous programme de calcul des pas elementaires
  194 +c ------------------------------------------------------------
  195 +c
  196 + call derivv(magin,magout,indgm,indval,year,tilt,
  197 + > rggsm,rgsmg,ds3,x4,y4,z4,r41,r42,r43,ier)
  198 +c
  199 + x5 = x + 1.5d0 * (r11 -3.d0 * r31 + 4.d0 * r41)
  200 + y5 = y + 1.5d0 * (r12 -3.d0 * r32 + 4.d0 * r42)
  201 + z5 = z + 1.5d0 * (r13 -3.d0 * r33 + 4.d0 * r43)
  202 +c
  203 +c ------------------------------------------------------------
  204 +c*FON 5ieme appel au sous programme de calcul des pas elementaires
  205 +c ------------------------------------------------------------
  206 +c
  207 + call derivv(magin,magout,indgm,indval,year,tilt,
  208 + > rggsm,rgsmg,ds3,x5,y5,z5,r51,r52,r53,ier)
  209 +c
  210 + errcur = abs(r11 - 4.5d0 * r31 + 4.d0 * r41 -.5d0 * r51) +
  211 + > abs(r12 - 4.5d0 * r32 + 4.d0 * r42 -.5d0 * r52) +
  212 + > abs(r13 - 4.5d0 * r33 + 4.d0 * r43 -.5d0 * r53)
  213 +c
  214 +c -------------------------------
  215 +c*FON Test sur la precision de calcul
  216 +c -------------------------------
  217 +c
  218 + if (errcur .lt. errin) then
  219 + x = x + .5d0 * (r11 + 4.d0 * r41 + r51)
  220 + y = y + .5d0 * (r12 + 4.d0 * r42 + r52)
  221 + z = z + .5d0 * (r13 + 4.d0 * r43 + r53)
  222 + else
  223 + ds = ds * .5d0
  224 + goto 10
  225 + endif
  226 +c
  227 + if (errcur.lt. errin * .04d0 .and. abs(ds) .lt. 1.33d0) then
  228 + ds = ds * 1.5d0
  229 + endif
  230 +c
  231 +c ****************
  232 +c Fin de programme
  233 +c ****************
  234 +c
  235 + return
  236 + end
... ...