subroutine chp45_70 (year,rre,thetr,phir,br,bt,bp,bb,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. CNES - JC KOSIK - janvier 1991 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 thetr (I) : colatitude geocentrique (radians) c*PAR phir (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 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, thetr, phir 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(120,6) c*LOC lg : coefficients entre 1945 et 1970 c integer kmax c*LOC kmax : nombre de pas de calcul (= 11) c integer jj,k,l,i,imin,ntot,imax,nm,inc c*LOC jj,k,l,i,imin,ntot,imax,nm,inc : indices de boucles et de tableaux c double precision sh(11,11) double precision gg(11,11) double precision hh(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 g(11,11),h(11,11) double precision tzero(6) double precision tmold,sph,st,temp,par double precision cph,ct,aor,ar,rrmag double precision delt,f1 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(i, 1), i=1, 120)/ > -30594,-2285, 5810,-1244, 2990,-1702, 1578, 477, 1282,-1834, > -499, 1255, 186, 913, -11, 944, 776, 144, 544, -276, > -421, -55, 304, -178, -253, 346, -12, 194, 95, -20, > -67, -142, -119, -82, 82, 59, 57, 6, 6, 100, > -246, 16, -25, -9, 21, -16, -104, -39, 70, -40, > -45, 0, -18, 0, 2, -29, 6, -10, 28, 15, > -17, 29, -22, 13, 7, 12, -8, -21, -5, -12, > 9, -7, 7, 2, -10, 18, 7, 3, 2, -11, > 5, -21, -27, 1, 17, -11, 29, 3, -9, 16, > 4, -3, 9, -4, 6, -3, 1, -4, 8, -3, > 11, 5, 1, 1, 2, -20, -5, -1, -1, -6, > 8, 6, -1, -4, -3, -2, 5, 0, -2, -2/ c data (lg(i, 2), i=1, 120)/ > -30554,-2250, 5815,-1341, 2998,-1810, 1576, 381, 1297,-1889, > -476, 1274, 206, 896, -46, 954, 792, 136, 528, -278, > -408, -37, 303, -210, -240, 349, 3, 211, 103, -20, > -87, -147, -122, -76, 80, 54, 57, -1, 4, 99, > -247, 33, -16, -12, 12, -12, -105, -30, 65, -55, > -35, 2, -17, 1, 0, -40, 10, -7, 36, 5, > -18, 19, -16, 22, 15, 5, -4, -22, -1, 0, > 11, -21, 15, -8, -13, 17, 5, -4, -1, -17, > 3, -7, -24, -1, 19, -25, 12, 10, 2, 5, > 2, -5, 8, -2, 8, 3, -11, 8, -7, -8, > 4, 13, -1, -2, 13, -10, -4, 2, 4, -3, > 12, 6, 3, -3, 2, 6, 10, 11, 3, 8/ c data (lg(i, 3), i=1, 120)/ > -30500,-2215, 5820,-1440, 3003,-1898, 1581, 291, 1302,-1944, > -462, 1288, 216, 882, -83, 958, 796, 133, 510, -274, > -397, -23, 290, -230, -229, 360, 15, 230, 110, -23, > -98, -152, -121, -69, 78, 47, 57, -9, 3, 96, > -247, 48, -8, -16, 7, -12, -107, -24, 65, -56, > -50, 2, -24, 10, -4, -32, 8, -11, 28, 9, > -20, 18, -18, 11, 9, 10, -6, -15, -14, 5, > 6, -23, 10, 3, -7, 23, 6, -4, 9, -13, > 4, 9, -11, -4, 12, -5, 7, 2, 6, 4, > -2, 1, 10, 2, 7, 2, -6, 5, 5, -3, > -5, -4, -1, 0, 2, -8, -3, -2, 7, -4, > 4, 1, -2, -3, 6, 7, -2, -1, 0, -3/ c data (lg(i, 4), i=1, 120)/ > -30421,-2169, 5791,-1555, 3002,-1967, 1590, 206, 1302,-1992, > -414, 1289, 224, 878, -130, 957, 800, 135, 504, -278, > -394, 3, 269, -255, -222, 362, 16, 242, 125, -26, > -117, -156, -114, -63, 81, 46, 58, -10, 1, 99, > -237, 60, -1, -20, -2, -11, -113, -17, 67, -56, > -55, 5, -28, 15, -6, -32, 7, -7, 23, 17, > -18, 8, -17, 15, 6, 11, -4, -14, -11, 7, > 2, -18, 10, 4, -5, 23, 10, 1, 8, -20, > 4, 6, -18, 0, 12, -9, 2, 1, 0, 4, > -3, -1, 9, -2, 8, 3, 0, -1, 5, 1, > -3, 4, 4, 1, 0, 0, -1, 2, 4, -5, > 6, 1, 1, -1, -1, 6, 2, 0, 0, -7/ c data (lg(i, 5), i=1, 120)/ > -30334,-2119, 5776,-1662, 2997,-2016, 1594, 114, 1297,-2038, > -404, 1292, 240, 856, -165, 957, 804, 148, 479, -269, > -390, 13, 252, -269, -219, 358, 19, 254, 128, -31, > -126, -157, -97, -62, 81, 45, 61, -11, 8, 100, > -228, 68, 4, -32, 1, -8, -111, -7, 75, -57, > -61, 4, -27, 13, -2, -26, 6, -6, 26, 13, > -23, 1, -12, 13, 5, 7, -4, -12, -14, 9, > 0, -16, 8, 4, -1, 24, 11, -3, 4, -17, > 8, 10, -22, 2, 15, -13, 7, 10, -4, -1, > -5, -1, 10, 5, 10, 1, -4, -2, 1, -2, > -3, 2, 2, 1, -5, 2, -2, 6, 4, -4, > 4, 0, 0, -2, 2, 3, 2, 0, 0, -6/ c data (lg(i,6),i=1,120)/ > -30220,-2068,5737,-1781,3000,-2047,1611, 25,1287,-2091, > -366, 1278, 251, 838,-196, 952, 800, 167, 461, -266, > -395, 26, 234, -279,-216, 359, 26, 262, 139, -42, > -139, -160, -91, -56, 83, 43, 64, -12, 15, 100, > -212, 72, 2, -37, 3, -6,-112, 1, 72, -57, > -70, 1, -27, 14, -4, -22, 8, -2, 23, 13, > -23, -2, -11, 14, 6, 7, -2, -15, -13, 6, > -3, -17, 5, 6, 0, 21, 11, -6, 3, -16, > 8, 10, -21, 2, 16, -12, 6, 10, -4, -1, > -5, 0, 10, 3, 11, 1, -2, -1, 1, -3, > -3, 1, 2, 1, -5, 3, -1, 4, 6, -4, > 4, 0, 1, -1, 0, 3, 3, 1, -1, -4/ c data kmax /11/ c data f1 /1.d0/ c data sp(1) /0.d0/ data p(1,1) /1.d0/ data cp(1) /1.d0/ data dp(1,1) /0.d0/ c data tzero /1945.d0, 1950.d0, 1955.d0, 1960.d0, > 1965.d0, 1970.d0/ c data tmold /0.d0/ c data sh(1,1) /0.d0/ c c ****************** c Debut de programme c ****************** c c --------------- c*FON Initialisations c --------------- c ier = 0 c if (sh(1,1) .ne. -1.d0) then sh(1,1) = -1.d0 c c --------------------------------------------------- c*FON Calcul des coefficients de shmidt et des constantes c --------------------------------------------------- c do 10 k = 1, kmax fn(k) = k do 11 l = 1, k fm(l) = l-1 const(k,l) = dble((k-2) * *2 - (l - 1) * *2) > / dble((2 * k - 3) * (2 * k - 5)) 11 continue 10 continue c do 15 k = 2, kmax sh(k,1) = dble(2 * k - 3) * sh(k-1,1) / dble(k-1) jj = 2 do 16 l = 2, k sh(k,l) = sh(k,l-1) * dsqrt(dble((k - l + 1) > * jj) / (k + l - 2)) jj = 1 16 continue 15 continue c endif c c ----------------------------------------------------------- c*FON Calcul des deux models de champ entourant la date de calcul c*FON et calcul de l'intervalle de temps pour l'interpolation c ----------------------------------------------------------- c if (year .ne. tmold) then tmold = year imin = idint((year - 1945.d0) / 5.d0) + 1 imax = imin + 1 delt = (year - tzero(imin)) / 5.d0 c ntot = 0 do 30 k = 2, kmax gg(k,1) = dble(lg(ntot+1,imin)) > + dble(lg(ntot+1,imax) > - lg(ntot+1,imin)) * delt l = 1 c nm = 2 * k-1 do 20 l = 2, k inc = 2 * l - 2 gg(k,l) = dble(lg(ntot+inc,imin)) > + dble(lg(ntot+inc,imax) > - lg(ntot+inc,imin)) * delt hh(k,l) = dble(lg(ntot+inc+1,imin)) > + dble(lg(ntot+inc+1,imax) > - lg(ntot+inc+1,imin)) * delt 20 continue ntot = ntot + nm 30 continue do 31 k = 2, kmax g(k,1) = gg(k,1) * sh(k,1) / f1 do 32 l = 2, k g(k,l) = gg(k,l) * sh(k,l) / f1 h(k,l) = hh(k,l) * sh(k,l) / f1 32 continue 31 continue endif c cph = dcos(phir) sph = dsin(phir) c if (thetr .ne. 0.d0) then ct = dcos(thetr) st = dsin(thetr) else ct = 1.0d0 st = 1.d-15 endif sp(2) = sph cp(2) = cph c do 90 k = 3, kmax sp(k) = sp(2) * cp(k-1) + cp(2) * sp(k-1) cp(k) = cp(2) * cp(k-1) - sp(2) * sp(k-1) 90 continue 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 c br = -(2.d0 * ar) * (g(2,1) * p(2,1) + p(2,2) > * (g(2,2) * cp(2) + h(2,2) * sp(2))) bt = ar * (g(2,1) * dp(2,1) + dp(2,2) > * (g(2,2) * cp(2) + h(2,2) * sp(2))) bp = ar * (h(2,2) * cp(2) - g(2,2) * sp(2)) * p(2,2) c k = 2 l = 12 c do 110 k = 3, kmax ar = aor * ar do 100 l = 1, k if (k .eq. l) then p(k,k) = st * p(k-1,k-1) dp(k,k) = st * dp(k-1,k-1) + ct * p(k-1,k-1) else p(k,l) = ct * p(k-1,l) - const(k,l) * p(k-2,l) dp(k,l) = ct * dp(k-1,l) - st * p(k-1,l) > - const(k,l) * dp(k-2,l) endif par = p(k,l) * ar if (l .eq. 1) then temp = g(k,l) else temp = g(k,l) * cp(l) + h(k,l) * sp(l) bp = bp-(g(k,l) * sp(l) - h(k,l) > * cp(l)) * fm(l) * par endif br = br - temp * fn(k) * par bt = bt + temp * dp(k,l) * ar 100 continue 110 continue c bp = bp / st c c --------------------------------------- c*FON Expression du champ magnetique en gauss c --------------------------------------- c bb = sqrt(br * br + bt * bt + bp * bp) c c **************** c Fin de programme c **************** c return end