subroutine dgrf05 (year,rre,thet,phi,br,bt,bp,bb,ier) c* c*********************************************************************** c* c* "Copyright [c] CNES 98 - tous droits reserves" 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 05.10.10 - V 4.0 c*VER 20.02.17 - V 5.0 c* c*AUT spec. CNES - JC KOSIK - fevrier 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 2005. c* c*PAR year (I) : annee fractionnaire (>= 2005.) 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 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*INF utilise : sans objet 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 - adoption des coefficients provisoires c*HST du champ 2005 c*HST version 4.0 - 05.10.10 - Passage aux coefficiens de Schmidt c*HST définitifs et aux variations seculaires c*HST interpolees c*HST Modif de fn1 (10->100) car précision c*HST des coeff 2005 sur 2 décimales c*HST Modif de fn2 (100->1000) car précision c*HST des coeff 2005 sur 3 décimales c*HST version 5.0 - 20.02.17 - Passage aux variations séculaires c* definitives 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 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 igrf 2005 c integer kmax c*LOC kmax : nombre de pas de calcul (= 14) c integer n,m,jj c*LOC n,m,jj : indices de boucles et de tableaux c double precision gg(14,14),ggt(14,14) double precision shmit(14,14) double precision 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 tmold,tzero,sph,st,t,temp,par,f2 double precision cph,ct,f1,aor,ar,rrmag 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 c Revision septembre 2010 - coefficients lg c Revision fevrier 2017 - coefficients lgt c data lg/100, > -2955463,-233724,133630,92055,-22700,7360,7988,2480,558,-217, > 295,-215,-16, > 507799,-166905,304769,-230583,79796,35441,6956,-7446,762,976, > -612,-160,-29,-88, > -259450,-51543,165776,124639,21065,20895,7674,-165,-1173,358, > 142,-188,21,30, > -19886,26972,-52472,67251,-37986,-13654,-15134,3873,-688,-694, > -235,144,89,28, > 28207,-22523,14515,-30536,10000,-16805,-1458,1230,-1811,501, > -15,-31,-38,-43, > 4272,18025,-12345,-1957,10385,-1355,1458,937,1017,-1076,306, > 29,96,118, > -2033,5475,6363,-6353,24,5094,-8636,542,936,-125,29,-79,-30, > -37, > -6114,-2257,682,2535,1093,-2632,-464,194,-1125,876,206,53,46, > 75, > 1120,-2088,983,-1971,1622,761,-1276,-6,-487,-666,377,180,-35, > -26, > -2011,1269,1267,-672,-816,810,292,-773,601,-922,-21,16,-36,35, > 219,10,446,476,-658,-101,-347,-86,-231,-793,-209,96,8,-5, > 26,144,-77,-227,90,-58,-269,-108,-158,-190,-139,399,-49,41, > -55,23,238,-263,61,40,1,2,28,-87,-34,88,-8,-10, > -76,33,172,-54,-107,-4,63,21,53,38,-22,-57,-82,-18 > / c data lgt/1000, > 11612,-11764,710,-1578,-774,-164,112,-78,-16,46,20,6,14, > -26746,16526,-4270,-4142,2202,576,-174,-108,118,-62,-24,24,16, > -2, > -22808,-12060,2082,-2858,-8814,-1738,-164,-580,-554,-26,-106, > -30,18,2, > 7692,-3594,-2462,-7756,4606,-902,1988,1302,258,334,256,42,30, > 28, > 882,2840,3862,-872,-2120,976,-1650,340,-246,-376,-2,-40,-50,-4, > 372,1752,1078,3912,-562,1104,-296,218,288,-324,-122,50, -2,-20, > -114,-2114,-418,-546,556,892,1654,-756,298,98,-124,0,38,12, > 668,274,-56,-78,-780,-258,272,596,-560,-66,14,-32,12,6, > -72,170,400,460,98,-130,404,340,266,-352,-136,-2,-8,16, > -86,-236,16,-84,148,-26,-156,330,200,-172,-164,-8,-2,6, > 108,-40,50,-64,-128,10,-96,-226,68,-76,-142,-42,26,14, > -26,46,22,102,-10,38,36,-38,-106,-8,-94,-48,-56,2, > -64, 8,-50,28,-24,38,-2,22,-2,2,22,-2,24,-32, > -22,-6,-12,-10,-14,-6,-18,-22,-8,12,-6,8,6,-16 > / c data f1 /100.d0/, f2 /1000.d0/ 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 /2005.d0/ c c ****************** c Debut de programme c ****************** c ier = 0 c c --------------------------------------------------------------- c*FON Initialisation a 0 des composantes du champ 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 ---------------------------------- 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 cph = cos(phi) sph = 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) = sph cp(2) = cph 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 rrmag = rre * rayt / rgmt aor = 1.d0 / rrmag ar = aor * aor * aor c p(2,1) = ct dp(2,1) = -st p(2,2) = st dp(2,2) = ct c br = - (ar + ar) * (g(2,1) * p(2,1) + p(2,2) * (g(2,2) * cp(2) + > g(1,2) * sp(2))) bt = ar * (g(2,1) * dp(2,1) + dp(2,2) * (g(2,2) * cp(2) + > g(1,2) * sp(2))) bp = ar * (g(1,2) * cp(2) - g(2,2) * sp(2)) * p(2,2) c c ----------------------------------------------------- c*FON Evaluation finale des composantes du champ magnetique c ----------------------------------------------------- c do 100 n = 3, kmax ar = aor * ar do 110 m = 1, n c if (m .eq. n) then p(n,n) = st * p(n-1,n-1) dp(n,n) = st * dp(n-1,n-1) + ct * p(n-1,n-1) else if (const(n,m) .eq. 0) then p(n,m) = ct * p(n-1,m) dp(n,m) = ct * dp(n-1,m) - st * p(n-1,m) else 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) endif endif c par = p(n,m) * ar c if (m .eq. 1) then temp = g(n,m) else temp = g(n,m) * cp(m) + g(m-1,n) * sp(m) bp = bp - (g(n,m) * sp(m) - g(m-1,n) * cp(m)) * > fm(m) * par endif c br = br - temp * fn(n) * par bt = bt + temp * dp(n,m) * ar c 110 continue c 100 continue c bp = bp / st c c --------------------------------------- c*FON Expression du champ magnetique en gauss c --------------------------------------- c br = br / 100000.d0 bt = bt / 100000.d0 bp = bp / 100000.d0 c bb = sqrt(bp * bp + bt * bt + br * br) c c **************** c Fin de programme c **************** c return end