subroutine chp95_20 (year,rre,thetr,phir,br,bt,bp,bb,ier) c* c*********************************************************************** c* c*PRO MAGLIB c* c*VER 20.07.13 - V 1.0 c* c*AUT spec. CNES - JC KOSIK - janvier 1991 c*AUT port. AKKA c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul du champ magnetique pour la periode c*ROL comprise entre 1995 et 2020. c* c*PAR year (I) : annee fractionnaire (entre 1995. et 2020.) 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 - 20.07.13 - Modele 2015 definitif et ajout c*HST du modele provisoire 2020 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 cos(-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(195,6) c*LOC lg : coefficients entre 1995 et 2020 c integer kmax c*LOC kmax : nombre de pas de calcul (= 14) c integer i,k,l,jj,imin,imax,ntot,nm,inc c*LOC i,k,l,jj,imin,imax,ntot,nm,inc : indices de boucles et de tableaux c double precision sh(14,14),gg(14,14),hh(14,14) double precision g(14,14),h(14,14),fn(14),fm(14) double precision cp(14),sp(14),p(14,14),dp(14,14),const(14,14) double precision tzero(6) double precision tmold, delt, cph, sph, ct, st, rrmag, ar, aor double precision temp, par, fmin, fmax 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 octobre 2006, nouveaux lg c Champ 1995 révisé octobre 2006 c data (lg(i,1),i=1,195)/ > -296920, -17840, 53060, -22000, 30700, -23660, 16810, -4130, > 13350,-22670, > -2620, 12490, 3020, 7590, -4270, 9400, 7800, 2620, 2900, -2360, > -4180, 970, 1220, -3060, -2140, 3520, 460, 2350, 1650, -1180, > -1430, -1660, -550, -170, 1070, 680, 670, -170, 680, 720, > -1700, 670, -10, -580, 190, 10, -930, 360, 770, -720, > -690, 10, -250, 280, 40, 50, 240, 40, 170, 80, > -240, -20, -60, 250, 60, 110, -60, -210, -90, 80, > -140, -230, 90, 150, 60, 110, -50, -160, -70, -40, > 40, 90, -200, 30, 150, -100, 120, 80, -60, -80, > -80, -10, 80, 100, 50, -20, -80, -80, 30, -30, > -60, 10, 20, 0, -40, 40, -10, 50, 40, -50, > 20, -10, 20, -20, 50, 10, 10, -20, 0, -70,75*0/ c c Champ 2000 révisé octobre 2006 c data (lg(i,2),i=1,195)/ > -296194,-17282,51861,-22677,30684,-24816,16709,-4580,13396, > -22880, > -2276,12521,2934,7145,-4911,9323,7868,2726,2500,-2319, > -4030,1198,1113,-3038,-2188,3514,438,2223,1719,-1304, > -1331,-1686,-393,-129,1063,723,682,-174,742,637, > -1609,651,-59,-612,169,7,-904,438,790,-740, > -646,0,-242,333,62,91,240,69,148,73, > -254,-12,-58,244,66,119,-92,-215,-79,85, > -166,-215,91,155,70,89,-79,-149,-70,-21, > 50,94,-197,30,134,-84,125,63,-62,-89, > -84,-15,84,93,38,-43,-82,-82,48,-26, > -60,17,17,0,-31,40,-5,49,37,-59, > 10,-12,20,-29,42,2,3,-22,-11,-74, > 27,-17,1,-19,13,15,-9,-1,-26,1, > 9,-7,-7,7,-28,17,-9,1,-12,12, > -19,40,-9,-22,-3,-4,2,3,9,25, > -2,-26,9,7,-5,3,3,0,-3,0, > -4,3,-1,-9,-2,-4,-4,8,-2,-9, > -9,3,2,1,18,-4,-4,13,-10,-4, > -1,7,7,-4,3,3,6,-1,3,4, > -2,0,-5,1,-9/ c c Champ 2005 coeff. definitifs révisés 2010.10.05 c data (lg(i,3),i=1,195)/ > -2955463,-166905,507799,-233724,304769,-259450,165776,-51543, > 133630,-230583,-19886,124639,26972,67251, > -52472,92055,79796,28207,21065,-22523,-37986,14515,10000, > -30536,-22700,35441,4272,20895, > 18025,-13654,-12345,-16805,-1957,-1355,10385,7360,6956,-2033, > 7674,5475,-15134,6363, > -1458,-6353,1458,24,-8636,5094,7988,-7446,-6114,-165,-2257, > 3873,682,1230, > 2535,937,1093,542,-2632,194,-464,2480,762,1120,-1173,-2088, > -688,983, > -1811,-1971,1017,1622,936,761,-1125,-1276,-487,-6,558,976, > -2011,358, > 1269,-694,1267,501,-672,-1076,-816,-125,810,876,292,-666,-773, > -922, > 601,-217,-612,219,142,10,-235,446,-15,476,306,-658,29,-101, > 206,-347,377,-86,-21,-231,-209,-793,295,-160,26,-188,144,144, > -77,-31,-227,29,90,-79,-58,53,-269,180,-108,16,-158,96, > -190,399,-139,-215,-29,-55,21,23,89,238,-38,-263,96,61, > -30,40,46,1,-35,2,-36,28,8,-87,-49,-34,-8,88, > -16,-88,-76,30,33,28,172,-43,-54,118,-107,-37,-4,75, > 63,-26,21,35,53,-5,38,41,-22,-10,-57,-18,-82 >/ c c Champ 2010 coeff. definitifs révisés 2017.02.23 c data (lg(i,4),i=1,195)/ > -2949657,-158642,494426,-239606,302634,-270854,166817,-57573, > 133985,-232654,-16040,123210,25175, > 63373,-53703,91266,80897,28648,16658,-21103,-35683,16446,8940, > -30972,-23087,35729,4458, > 20026,18901,-14105,-11806,-16317,-1,-803,10104,7278,6869, > -2090,7592,4418,-14140, > 6154,-2283,-6626,1310,302,-7809,5540,8044,-7500,-5780,-455, > -2120,4524,654, > 1400,2496,1046,703,164,-2761,492,-328,2441,821,1084,-1450, > -2003,-559, > 1183,-1934,-1741,1161,1671,1085,696,-1405,-1074,-354,164,550, > 945,-2054, > 345,1151,-527,1275,313,-714,-1238,-742,-76,797,843,214,-842, > -608, > -1008,701,-194,-624,273,89,-10,-107,471,-16,444,245,-722,-33, > -96,213,-395,309,-199,-103,-197,-280,-831,305,-148,13,-203,167, > 165,-66,-51,-176,54,85,-79,-39,37,-251,179,-127,12,-211, > 75,-194,375,-186,-212,-21,-87,30,27,104,213,-63,-249,95, > 49,-11,59,52,0,-39,13,-37,27,21,-86,-77,-23,4, > 87,-9,-89,-87,31,30,42,166,-45,-59,108,-114,-31,-7, > 78,54,-18,10,38,49,2,44,42,-25,-26,-53,-26,-79 > / c c Champ 2015 coeff. definitifs révisés 2020.07.13 c data (lg(i,5),i=1,195)/ > -2944146,-150177,479599,-244588,301220,-284541,167635,-64217, > 135033,-235226,-11529,122585,24504, > 58169,-53870,90742,81368,28354,12049,-18843,-33485,18095,7038, > -32923,-23291,36014,4698, > 19235,19698,-14094,-11914,-15740,1598,430,10012,6955,6757,-2061, > 7279,3330,-12985, > 5874,-2893,-6664,1314,735,-7085,6241,8129,-7599,-5427,-679, > -1953,5182,559, > 1507,2445,932,327,-288,-2750,661,-232,2398,889,1004,-1678,-1826, > -316, > 1318,-2056,-1460,1333,1616,1176,569,-1598,-910,-202,226,533,883, > -2177, > 302,1076,-322,1174,67,-674,-1320,-688,-10,779,868,104,-906,-389, > -1054,844,-201,-626,328,17,-40,55,455,-55,440,170,-792,-67, > -61,213,-416,233,-285,-180,-112,-359,-872,300,-140,0,-230,211, > 208,-60,-79,-105,58,76,-70,-20,14,-212,170,-144,-22,-257, > 44,-201,349,-234,-209,-16,-108,46,37,123,175,-89,-219,85, > 27,10,72,54,-9,-37,29,-43,23,22,-89,-94,-16,-3, > 72,-2,-92,-88,42,49,63,156,-42,-50,96,-124,-19,-10, > 81,42,-13,-4,38,48,8,48,46,-30,-35,-43,-36,-71 > / c c Champ 2020 coeff. provisoires révisés 2020.07.13 c data (lg(i,6),i=1,195)/ > -2940480,-145090,465250,-249960,298200,-299160,167700,-73460, > 136320,-238120,-8210,123620,24190, > 52570,-54340,90300,80950,28190,8630,-15840,-30940,19970,4800, > -34970,-23430,36320,4770, > 18780,20830,-14070,-12120,-15120,3230,1350,9890,6600,6550, > -1910,7290,2510,-12150, > 5280,-3620,-6450,1350,890,-6470,6810,8060,-7670,-5150,-820, > -1690,5650,220, > 1580,2350,640,-220,-720,-2720,980,-180,2370,970,840,-1760,-1530, > -50, > 1280,-2110,-1170,1530,1490,1370,360,-1650,-690,-30,280,500,840, > -2340, > 290,1100,-150,980,-110,-510,-1320,-630,110,780,880,40,-930,-140, > -1190,960,-190,-620,340,-10,-20,170,360,-90,480,70,-860,-90, > -10,190,-430,140,-340,-240,-10,-380,-880,300,-140,00,-250,250, > 230,-60,-90,-40,30,60,-70,-20,-10,-170,140,-160,-60,-300, > 20,-200,310,-260,-200,-10,-120,50,50,130,140,-120,-180,70, > 10,30,80,50,-20,-30,60,-50,20,10,-90,-110,00,-30, > 50,10,-90,-90,50,60,70,140,-30,-40,80,-130,00,-10, > 80,30,00,-10,40,50,10,50,50,-40,-50,-40,-40,-60 > / c data kmax /14/ c 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 /1995.d0,2000.d0,2005.d0,2010.d0,2015.d0,2020.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 20 k = 1, kmax fn(k) = k do 10 l = 1, k fm(l) = l - 1 const(k,l) = dble((k-2)**2 - (l - 1)**2) > / dble((2 * k - 3) * (2 * k - 5)) 10 continue 20 continue c do 40 k = 2, kmax sh(k,1) = dble(2 * k - 3) * sh(k-1,1) / dble(k-1) jj = 2 do 30 l = 2, k sh(k,l) = sh(k,l-1) * dsqrt(dble((k - l + 1) * jj) > / (k + l - 2)) jj = 1 30 continue 40 continue c endif c c ------------------------------------------------------------ c*FON Calcul des deux modeles 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 - 1995.d0) / 5.d0) + 1 imax = imin + 1 c *********************************************** c Facteurs de normalisation des modèle min et max c *********************************************** if (tzero(imin) .ge. 2005.d0) then fmin = 100.d0 else fmin = 10.d0 endif if (tzero(imax) .ge. 2005.d0) then fmax = 100.d0 else fmax = 10.d0 endif delt = (year - tzero(imin)) / 5.d0 c ntot = 0 do 60 k = 2, kmax gg(k,1) = dble(lg(ntot+1,imin))/fmin > + ((dble(lg(ntot+1,imax))/fmax) > - (dble(lg(ntot+1,imin))/fmin)) > * delt l = 1 c nm = 2 * k - 1 do 50 l = 2, k inc = 2 * l - 2 gg(k,l) = dble(lg(ntot+inc,imin))/fmin > + ((dble(lg(ntot+inc,imax))/fmax) > - (dble(lg(ntot+inc,imin))/fmin)) > * delt hh(k,l) = dble(lg(ntot+inc+1,imin))/fmin > + ((dble(lg(ntot+inc+1,imax))/fmax) > - (dble(lg(ntot+inc+1,imin))/fmin)) > * delt c 50 continue ntot = ntot + nm 60 continue do 80 k = 2, kmax g(k,1) = gg(k,1) * sh(k,1) do 70 l = 2, k g(k,l) = gg(k,l) * sh(k,l) h(k,l) = hh(k,l) * sh(k,l) 70 continue 80 continue c 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 c sp(k) = sp(2) * cp(k-1) + cp(2) * sp(k-1) cp(k) = cp(2) * cp(k-1) - sp(2) * sp(k-1) c 90 continue 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 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 c ar = aor * ar c do 100 l = 1, k c 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 c par = p(k,l) * ar c 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 c 110 continue c bp = bp / st c bb = dsqrt(br * br + bt * bt + bp * bp) c c **************** c Fin de programme c **************** c return end