subroutine dgrf45_70 (year,rre,thet,phi,br,bt,bp,b,ier) c* c*********************************************************************** c* c*PRO MAGLIB c* c*VER 99.03.31 - V 1.0 c*VER 01.06.05 - V 2.0 c*VER 03.01.06 - V 2.1 c*VER 01.01.07 - V 3.0 c* c*AUT spec. N A S A c*AUT spec. CNES - JC KOSIK - janvier 1996 c*AUT port. CISI c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul du champ magnetique pour la periode c*ROL comprise entre 1945 et 1970. c* c*PAR year (I) : annee fractionnaire (entre 1945. et 1970.) c* c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) c*PAR thet (I) : colatitude geocentrique (radians), la colatitude c*PAR : est calculee du pole nord au pole sud c*PAR phi (I) : longitude geocentrique (radians), calculee c*PAR : positive vers l'est de 0 a 2*PI 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 b (O) : module du champ magnetique (gauss) c* c*PAR ier (O) : code de retour c* c*NOT ier : sans objet c c*NOT Un test est effectue pour les colatitudes inferieurs a 1.d-15 c*NOT ou superieures a PI - 1.d-15 et le cosinus est remplace par 1.d0 c*NOT et le sinus par 1.d-15 ou bien le cosinus est remplace par -1.d0 c*NOT et le sinus par 1.d-15. c* c*NOT common : util, util2 c* c*INF utilise : sans objet c* c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP c*HST version 2.0 - 01.06.05 - 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 - remplace dgrfold de la version 2.1 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 double precision rre, thet, phi double precision br, bt, bp, b 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(11,11,6) c*LOC lg : coefficients entre 1945 et 1970 c integer kmax c*LOC kmax : nombre de pas de calcul (= 11) c integer n,m,jj,km c*LOC n,m,jj,km : indices de boucles et de tableaux c double precision cof1,cof2,delt double precision shmit(11,11) double precision g(11,11) double precision const(11,11),fn(11),fm(11) double precision p(11,11),dp(11,11),sp(11),cp(11) double precision tzero(6) double precision tmold,sph,st,temp,par double precision cph,ct,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* coefficients pour 1945,1950,1955,1960,1965,1970 c data lg/ > 1,-30594,-1244, 1282, 944, -253, 59, 70, 13, 5, -3, > 5810, -2285, 2990,-1834, 776, 346, 57, -40, 7, -21, 11, > -1702, 477, 1578, 1255, 544, 194, 6, 0, -8, 1, 1, > -499, 186, -11, 913, -421, -20, -246, 0, -5, -11, 2, > 144, -276, -55, -178, 304, -142, -25, -29, 9, 3, -5, > -12, 95, -67, -119, 82, -82, 21, -10, 7, 16, -1, > 6, 100, 16, -9, -16, -39, -104, 15, -10, -3, 8, > -45, -18, 2, 6, 28, -17, -22, 29, 7, -4, -1, > 12, -21, -12, -7, 2, 18, 3, -11, 2, -3, -3, > -27, 17, 29, -9, 4, 9, 6, 1, 8, -4, 5, > 5, 1, -20, -1, -6, 6, -4, -2, 0, -2, -2, > > 1,-30554,-1341, 1297, 954, -240, 54, 65, 22, 3, -8, > 5815, -2250, 2998,-1889, 792, 349, 57, -55, 15, -7, 4, > -1810, 381, 1576, 1274, 528, 211, 4, 2, -4, -1, -1, > -476, 206, -46, 896, -408, -20, -247, 1, -1, -25, 13, > 136, -278, -37, -210, 303, -147, -16, -40, 11, 10, -4, > 3, 103, -87, -122, 80, -76, 12, -7, 15, 5, 4, > -1, 99, 33, -12, -12, -30, -105, 5, -13, -5, 12, > -35, -17, 0, 10, 36, -18, -16, 19, 5, -2, 3, > 5, -22, 0, -21, -8, 17, -4, -17, -1, 3, 2, > -24, 19, 12, 2, 2, 8, 8, -11, -7, 8, 10, > 13, -2, -10, 2, -3, 6, -3, 6, 11, 8, 3, > > 1,-30500,-1440, 1302, 958, -229, 47, 65, 11, 4, -3, > 5820, -2215, 3003,-1944, 796, 360, 57, -56, 9, 9, -5, > -1898, 291, 1581, 1288, 510, 230, 3, 2, -6, -4, -1, > -462, 216, -83, 882, -397, -23, -247, 10, -14, -5, 2, > 133, -274, -23, -230, 290, -152, -8, -32, 6, 2, -3, > 15, 110, -98, -121, 78, -69, 7, -11, 10, 4, 7, > -9, 96, 48, -16, -12, -24, -107, 9, -7, 1, 4, > -50, -24, -4, 8, 28, -20, -18, 18, 6, 2, -2, > 10, -15, 5, -23, 3, 23, -4, -13, 9, 2, 6, > -11, 12, 7, 6, -2, 10, 7, -6, 5, 5, -2, > -4, 0, -8, -2, -4, 1, -3, 7, -1, -3, 0, > > 1,-30421,-1555, 1302, 957, -222, 46, 67, 15, 4, 1, > 5791, -2169, 3002,-1992, 800, 362, 58, -56, 6, 6, -3, > -1967, 206, 1590, 1289, 504, 242, 1, 5, -4, 0, 4, > -414, 224, -130, 878, -394, -26, -237, 15, -11, -9, 0, > 135, -278, 3, -255, 269, -156, -1, -32, 2, 1, -1, > 16, 125, -117, -114, 81, -63, -2, -7, 10, 4, 4, > -10, 99, 60, -20, -11, -17, -113, 17, -5, -1, 6, > -55, -28, -6, 7, 23, -18, -17, 8, 10, -2, 1, > 11, -14, 7, -18, 4, 23, 1, -20, 8, 3, -1, > -18, 12, 2, 0, -3, 9, 8, 0, 5, -1, 2, > 4, 1, 0, 2, -5, 1, -1, 6, 0, -7, 0, > > 1,-30334,-1662, 1297, 957, -219, 45, 75, 13, 8, -2, > 5776, -2119, 2997,-2038, 804, 358, 61, -57, 5, 10, -3, > -2016, 114, 1594, 1292, 479, 254, 8, 4, -4, 2, 2, > -404, 240, -165, 856, -390, -31, -228, 13, -14, -13, -5, > 148, -269, 13, -269, 252, -157, 4, -26, 0, 10, -2, > 19, 128, -126, -97, 81, -62, 1, -6, 8, -1, 4, > -11, 100, 68, -32, -8, -7, -111, 13, -1, -1, 4, > -61, -27, -2, 6, 26, -23, -12, 1, 11, 5, 0, > 7, -12, 9, -16, 4, 24, -3, -17, 4, 1, 2, > -22, 15, 7, -4, -5, 10, 10, -4, 1, -2, 2, > 2, 1, 2, 6, -4, 0, -2, 3, 0, -6, 0, > > 1,-30220,-1781, 1287, 952, -216, 43, 72, 14, 8, -3, > 5737, -2068, 3000,-2091, 800, 359, 64, -57, 6, 10, -3, > -2047, 25, 1611, 1278, 461, 262, 15, 1, -2, 2, 2, > -366, 251, -196, 838, -395, -42, -212, 14, -13, -12, -5, > 167, -266, 26, -279, 234, -160, 2, -22, -3, 10, -1, > 26, 139, -139, -91, 83, -56, 3, -2, 5, -1, 6, > -12, 100, 72, -37, -6, 1, -112, 13, 0, 0, 4, > -70, -27, -4, 8, 23, -23, -11, -2, 11, 3, 1, > 7, -15, 6, -17, 6, 21, -6, -16, 3, 1, 0, > -21, 16, 6, -4, -5, 10, 11, -2, 1, -1, 3, > 1, 1, 3, 4, -4, 0, -1, 3, 1, -4, -1/ c data kmax /11/ 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 /1945.d0, 1950.d0, 1955.d0, 1960.d0, > 1965.d0, 1970.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 b = 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, 11 fn(n) = n do 20 m = 1, 11 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, 11 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 Traitements effectues si c'est le 1er appel c ------------------------------------------- c endif c c -------------------------------------------------- c*FON Calcul du g(n,m), interpolation entre deux jeux de c*FON coefficients avec un delai de 5 ans entre eux c -------------------------------------------------- c if (year .ne. tmold) then c tmold = year c km = idint((year - 1945.d0) / 5.d0) + 1 c delt = (year - tzero(km)) / 5.d0 c do 70 n = 1, kmax do 80 m = 1, kmax cof1 = dble(lg(n,m,km)) cof2 = dble(lg(n,m,km+1)) g(n,m) = cof1 + (cof2 - cof1) * delt 80 continue 70 continue c do 77 n = 1, kmax do 78 m = 1, kmax g(n,m) = g(n,m) * shmit(n,m) 78 continue 77 continue c endif c c ---------------------------------------------------- c*FON Calcul des cosinus directeurs de la colatitude et de c*FON la longitude. On test si la colatitude est comprise c*FON entre 0 et PI 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 Calcul des trois premiers termes c -------------------------------- c rrmag = rre * rayt / rgmt aor = 1.d0 / rrmag ar = aor * aor * aor p(2,1) = ct dp(2,1) = -st p(2,2) = st dp(2,2) = ct 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 Calcul final 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 b = sqrt(bp * bp + bt * bt + br * br) c c **************** c Fin de programme c **************** c return end