subroutine grad00 (year,rre,thet,phi,br,bt,bp,bb,gradb,ier) c* c*********************************************************************** c* c*PRO MAGLIB c* 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*VER 14.12.10 - V 4.0 c* c*AUT spec. CNES - JC KOSIK - mars 2001 c*AUT port. CS-SI c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul du champ magnetique d'origine interne avec les c*ROL coefficients du modele DGRF 2000. c* c*PAR year (I) : annee fractionnaire (>= 2000.) c* c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) c*PAR thet (I) : colatitude geocentrique (radians) c*PAR phi (I) : longitude geocentrique (radians) c* c*PAR br (O) : composante radiale du champ magnetique le long du c* : meridien positive vers l'exterieur (gauss) c*PAR bt (O) : composante tangentielle du champ magnetique le long c*PAR : du meridien positive vers le sud (gauss) c*PAR bp (O) : composante azimuthale du champ magnetique, positive c*PAR : vers l'est (gauss) c* c*PAR bb (O) : module du champ magnetique (gauss) c* c*PAR gradb (O) : gradients du champ magnetique c*PAR : les quantites sont evaluees en unite de rayon terrestre c* c*PAR ier (O) : code de retour c* c*NOT ier : sans objet c* c*NOT Dans le cas ou l'angle thet est nul , son sinus est c*NOT remplacee par 1.d-15 afin d'eviter des divisions par 0 et c*NOT son cosinus est remplace par 1.0d0. c*NOT Dans le cas ou l'angle thet est egal a pi, son sinus est c*NOT remplacee par 1.d-15 afin d'eviter des divisions par 0 et c*NOT son cosinus est remplace par -1.0d0. c* c*NOT common : util, util2 c* c*HST version 2.0 - 01.05.30 - Enrichissement de la maglib au CDPP c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 c*HST version 3.0 - 01.01.07 - utilisation du modele DGRF 2000 c*HST version 4.0 - 14.12.10 - Mise à jour des VS 2000. 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,rre,thet,phi double precision br,bt,bp,bb double precision gradb(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 double precision rgmt c c*COM rgmt : rayon geocentrique terrestre (km) c common/util2/rgmt c c --------------------------------- c*FON Declaration des variables locales c --------------------------------- c integer lg(14,14),lgt(14,14) c*LOC lg, lgt : coefficients du modele dgrf05 c integer n,m,jj c*LOC n,m,jj : indices de boucles et de tableaux c integer kmax c*LOC kmax : nombre de pas de calcul (= 14) c double precision dbrdr,dbrdt,dbrdp,dbtdr,dbtdt,dbtdp double precision dbpdr,dbpdt,dbpdp c*LOC dbrdr,dbrdt,dbrdp : derivees partielles des composantes radiale, c*LOC dbtdr,dbtdt,dbtdp : tangentielle et azimuthale du champ total c*LOC dbpdr,dbpdt,dbpdp : c double precision gg(14,14),ggt(14,14),shmit(14,14),g(14,14) double precision const(14,14),fn(14),fm(14) double precision p(14,14),dp(14,14),sp(14),cp(14) double precision d2p(14,14) double precision tmold,tzero,sphi,st,t,f2 double precision cphi,ct,f1,aor,aor3,aorn,daorn,rmag,rkm double precision temp1,dtemp1,d2temp1 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 lg/10, > -296194,-22677,13396,9323,-2188,723,790,244,50,-26,27,-22,-2, > 51861,-17282,30684,-22880,7868,3514,682,-740,66,94,-60,-17,-3, > -9, > -24816,-4580,16709,12521,2500,2223,742,0,-92,30,17,-19,2,3, > -2276,2934,-4911,7145,-4030,-1304,-1609,333,-79,-84,-31,15,9,1, > 2726,-2319,1198,-3038,1113,-1686,-59,91,-166,63,-5,-1,-2,-4, > 438,1719,-1331,-393,1063,-129,169,69,91,-89,37,1,9,13, > -174,637,651,-612,7,438,-904,73,70,-15,10,-7,-5,-4, > -646,-242,62,240,148,-254,-58,-12,-79,93,20,7,3,7, > 119,-215,85,-215,155,89,-149,-21,-70,-43,42,17,-3,-4, > -197,134,125,-62,-84,84,38,-82,48,-82,3,1,-4,3, > 17,0,40,49,-59,-12,-29,2,-22,-74,-11,12,-1,-1, > 1,13,-9,-26,9,-7,-28,-9,-12,-19,-9,40,-2,4, > -4,3,25,-26,7,3,0,0,3,-9,-4,8,-4,0, > -9,2,18,-4,-10,-1,7,3,6,3,-2,-5,-9,1 > / c data lgt/1000, > 12954,-13908,-660,-2350,-1640,260,176,80,116,86,50,10,8, > -21622,11830,-4142,-3566,2232,602,272,-92,204,72,-24,20,2,4, > -22580,-11486,-2628,-1142,-7870,-2670,508,-330,-506,116,-56,4,2, > 0, > 5748,-4736,-6724,-8398,4628,-1228,1912,1086,204,292,150,-12,-2, > 36, > 1894,1334,5070,-312,-2260,110,-1736,640,-302,-258,70,-42,-36,-6, > -216,1670,1930,3946,-490,-130,-464,494,214,-372,-128,38,12,-24, > -586,-1790,-294,-466,-92,1428,808,-376,472,50,-142,-18,40,6, > 692,326,124,270,-774,-184,232,628,-670,-108,12,-34,32,10, > -140,124,266,358,144,-258,428,408,426,-472,-86,20,-10,28, > -82,-142,34,-104,48,-60,-176,94,242,-204,-102,12,8,10, > 98,20,92,-28,-136,38,-114,-212,-22,-106,-198,-48,36,10, > 32,28,26,66,0,24,22,-36,-76,0,-98,-2,-58,2, > -30,-14,-24,-6,-18,20,2,4,-4,6,12,16,64,-20, > 28,26,-16,-28,-14,12,-14,-18,-14,16,-4,-14,16,-56 > / c data kmax /14/ c data p(1,1) /1.d0/ data cp(1) /1.d0/ data dp(1,1) /0.d0/ data sp(1) /0.d0/ c data shmit(1,1)/0.d0/ c data tmold/0.d0/ c data tzero/2000.d0/ c c ****************** c Debut de programme c ****************** c ier = 0 c c -------------------------------------------- c*FON Initialisation a 0. des composantes du champ c*FON ainsi que du module c -------------------------------------------- c br = 0.0d0 bt = 0.0d0 bp = 0.0d0 bb = 0.0d0 c c ------------------------------------------------- c*FON Initialisations faites au premier appel du module c ------------------------------------------------- c if (shmit(1,1) .ne. -1.d0) then c c mise a zero des constantes avant la premiere initialisation c do 8 n = 1, 14 fn(n) = 0.d0 do 9 m = 1, 14 fm(m) = 0.d0 const(n,m) = 0.d0 shmit(n,m) = 0.d0 dp(n,m) = 0.d0 d2p(n,m) = 0.d0 9 continue 8 continue c c ---------------------------------- c*FON Calcul des tableaux fn,fm et const c ---------------------------------- c shmit(1,1) = -1.d0 do 10 n = 1, 14 fn(n) = n do 20 m = 1, 14 fm(m) = m - 1 const(n,m) = dble((n - 2)**2 - (m - 1)**2) / > dble((2 * n - 3) * (2 * n - 5)) 20 continue 10 continue c c ---------------------------- c*FON Calcul du tableau shmit(n,m) c ---------------------------- c do 30 n = 2, 14 shmit(n,1) = dble(2 * n - 3) * shmit(n-1,1) / dble(n - 1) jj = 2 do 40 m = 2, n shmit(n,m) = shmit(n,m-1) * > sqrt(dble((n - m + 1) * jj) / (n + m - 2)) shmit(m-1,n) = shmit(n,m) jj = 1 40 continue 30 continue c c --------------------------------------- c*FON Calcul des tableaux gg(n,m) et ggt(n,m) c --------------------------------------- c f1 = dble(lg(1,1)) f2 = dble(lgt(1,1)) do 50 n = 1, kmax do 60 m = 1, kmax gg(n,m) = dble(lg(n,m)) * shmit(n,m) / f1 ggt(n,m) = dble(lgt(n,m)) * shmit(n,m) / f2 60 continue 50 continue c c -------------------------------------------------- c*FON Traitements effectues si ce n'est pas le 1er appel c -------------------------------------------------- c endif c c ------------------------ c*FON Calcul du tableau g(n,m) c ------------------------ c if (year .ne. tmold) then c tmold = year t = year - tzero c do 70 n = 1, kmax do 80 m = 1, kmax g(n,m) = gg(n,m) + t * ggt(n,m) 80 continue 70 continue c endif c c ---------------------------------------------------- c*FON Calcul des cosinus directeurs de la colatitude et de c*FON la longitude geocentriques c ---------------------------------------------------- c cphi = cos(phi) sphi = sin(phi) c if (thet .ne. 0.0d0 .and. thet .ne. acos(-1.0d0)) then c ct = cos(thet) st = sin(thet) c else c if (thet .ne. acos(-1.0d0)) then ct = 1.0d0 st = 1.d-15 elseif (thet .ne. 0.0d0) then ct = -1.0d0 st = 1.d-15 endif c endif c sp(2) = sphi cp(2) = cphi c c ---------------------------- c*FON Calcul des tableaux cp et sp c ---------------------------- c do 90 m = 3, kmax sp(m) = sp(2) * cp(m-1) + cp(2) * sp(m-1) cp(m) = cp(2) * cp(m-1) - sp(2) * sp(m-1) 90 continue c c --------------------------------------------------- c*FON Premiers termes des composantes du champ magnetique c --------------------------------------------------- c rkm = rre * rayt rmag = rkm / rgmt c aor = rgmt / rkm aor3 = aor * aor * aor c p(2,1) = ct dp(2,1) = -st d2p(2,1) = -ct p(2,2) = st dp(2,2) = ct d2p(2,2) = -st c br = -2.d0 * aor3 * > (g(2,1) * ct + (g(2,2) * cphi + g(1,2) * sphi) * st) c bt = aor3 * > (-g(2,1) * st + (g(2,2) * cphi + g(1,2) * sphi) * ct) c bp = aor3 * (g(1,2) * cphi - g(2,2) * sphi) c c ------------------------------------------ c*FON Calcul des 9 derivees partielles de depart c ------------------------------------------ c dbrdr = 6.d0 * aor3 * > (g(2,1) * ct + (g(2,2) * cphi + g(1,2) * sphi) * st) / rre c dbrdt = -2.d0 * aor3 * > (-g(2,1) * st + (g(2,2) * cphi + g(1,2) * sphi) * ct) c dbrdp = -2.d0 * aor3 * > (g(1,2) * cphi - g(2,2) * sphi) * st c dbtdr = -3.d0 * aor3 * > (-g(2,1) * st + (g(2,2) * cphi + g(1,2) * sphi) * ct) / rre c dbtdt = aor3 * > (-g(2,1) * ct - (g(2,2) * cphi + g(1,2) * sphi) * st) c dbtdp = aor3 * > (-g(2,2) * sphi + g(1,2) * cphi) * ct c dbpdr = -3.d0 * aor3 * > (g(1,2) * cphi - g(2,2) * sphi) / rre c dbpdt = 0.d0 c dbpdp = aor3 * (-g(2,2) * cphi - g(1,2) * sphi) c aorn = aor3 c c ----------------------------------------------------- c*FON Evaluation finale des composantes du champ magnetique c ----------------------------------------------------- c do 100 n = 3, kmax c aorn = aor * aorn daorn = -dble(n + 1) * aorn / rre c do 110 m = 1, n c if (m .eq. n) then c p(n,n) = st * p(n-1,n-1) dp(n,n) = st * dp(n-1,n-1) + ct * p(n-1,n-1) d2p(n,n) = ct * dp(n-1,n-1) + st * d2p(n-1,n-1) > - st * p(n-1,n-1) + ct * dp(n-1,n-1) c else c p(n,m) = ct * p(n-1,m) - const(n,m) * p(n-2,m) dp(n,m) = ct * dp(n-1,m) - st * p(n-1,m) - const(n,m) > * dp(n-2,m) d2p(n,m) = - st * dp(n-1,m) + ct * d2p(n-1,m) > - ct * p(n-1,m) - st * dp(n-1,m) > - const(n,m) * d2p(n-2,m) c endif c if (m .eq. 1) then c temp1 = g(n,m) dtemp1 = 0.d0 d2temp1 = 0.d0 c else c temp1 = g(n,m) * cp(m) + g(m-1,n) * sp(m) dtemp1 = fm(m) * (-g(n,m) * sp(m) + g(m-1,n) * cp(m)) d2temp1 = -fm(m) * fm(m) * temp1 c endif c br = br - temp1 * p(n,m) * aorn * fn(n) bt = bt + temp1 * dp(n,m) * aorn bp = bp + dtemp1 * p(n,m) * aorn / st c dbrdr = dbrdr - temp1 * p(n,m) * daorn * fn(n) dbrdt = dbrdt - temp1 * dp(n,m) * aorn * fn(n) dbrdp = dbrdp - dtemp1 * p(n,m) * aorn * fn(n) c dbtdr = dbtdr + temp1 * dp(n,m) * daorn dbtdt = dbtdt + temp1 * d2p(n,m) * aorn dbtdp = dbtdp + dtemp1 * dp(n,m) * aorn c dbpdr = dbpdr + dtemp1 * p(n,m) * daorn / st dbpdt = dbpdt + dtemp1 * aorn * > (dp(n,m) - ct * p(n,m) / st) / st dbpdp = dbpdp + d2temp1 * p(n,m) * aorn / st c 110 continue c 100 continue c c mise en memoire des gradients de champ c gradb(1,1) = dbrdr gradb(1,2) = dbrdt gradb(1,3) = dbrdp gradb(2,1) = dbtdr gradb(2,2) = dbtdt gradb(2,3) = dbtdp gradb(3,1) = dbpdr gradb(3,2) = dbpdt gradb(3,3) = dbpdp c bb = sqrt(bp * bp + bt * bt + br * br) c c **************** c Fin de programme c **************** c return end