subroutine gsfc65 (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 99.03.31 - V 1.0 c* c*AUT spec. GODDARD Space Flight Center c*AUT adap. CNES - JC KOSIK - janvier 1991 c*AUT port. CISI c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul du champ magnetique terrestre avec les c*ROL coefficients de 1965 correspondant au programme c*ROL original de Mac Ilwain. c* c*PAR year (I) : annee fractionnaire (>= 1965.) 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 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* 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(13,13),lgt(13,13) c*LOC lg, lgt : coefficients 1965 correspondant au programme c*LOC : original de Mac Ilwain c integer kmax c*LOC kmax : nombre de pas de calcul (= 10) c integer n,m,jj c*LOC n,m,jj : indices de boucles et de tableaux c double precision gg(13,13),ggt(13,13) double precision shmit(13,13) double precision g(13,13) double precision const(13,13),fn(13),fm(13) double precision p(13,13),dp(13,13),sp(13),cp(13) double precision aor,ar,cph,ct,f1,f2,par,sph,st,t double precision temp,tmold,tzero,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 data lg/ > 10,-304249,-15361, 13009, 9576,-2277, 498, 709, 48, 99,3*0, >57748, -21616, 30002,-19870, 8028, 3595, 607,-572, 67, 29,3*0, >-19498, 2043, 15853, 12904, 5026, 2313, 45, 56, -88, 74,3*0, > -4310, 2308, -1300, 8712,-3940, -312,-2417, 75,-138,-156,3*0, > 1520, -2684, 29, -2505, 2714,-1573, -12,-244, -33, 114,3*0, > 86, 1212, -1160, -1104, 799, -652, 5, -15, 71,111,3*0, > -119, 1028, 609, -272, -124, -116,-1091, 141, -56, 10,3*0, > -540, -244, -91, 22, 276, -211, -201, 58, 117, 4*0, > 69, -122, 58, -170, 26, 236, -25, -160, 64, 16,3*0, > -220, 156, 51, -35, -18, 96, 121, 2, -25, 15,3*0, > 39*0/ c data lgt/ > 100, 2059, -2907, 266, -86, 255, -70, 6*0, > -394, 602, 121, -1003, 194, -8, 99, 6*0, > -1369, -1578, -70, 163, -117, 153, 85, 6*0, > 649, 293, -924, -130, -54, -42, 211, 6*0, > -177, -154, 318, -548, -417, -72, 157, 6*0, > 304, 288, -186, 125, 80, 164, -9, 6*0, > -139, 12, 153, -73, -6, 45, 6, 6*0, > 78*0/ c data kmax /10/ c data p(1,1) /1.0d0/ data cp(1) /1.0d0/ data dp(1,1) /0.0d0/ data sp(1) /0.0d0/ c data shmit(1,1) /0.d0/ c data tmold /0.d0/ c data tzero /1960.d0/ c c ****************** c Debut de programme c ****************** c ier = 0 c if (shmit(1,1) .eq. -1.d0) go to 40 c c -------------------------------------------- c*FON Initialisation des tableaux au premier appel c -------------------------------------------- c shmit(1,1) = -1.d0 c do 10 n = 1, 13 fn(n) = dble(n) do 15 m = 1, 13 fm(m) = dble(m - 1) const(n,m) = dble((n - 2)**2 - (m - 1)**2) > / ((2 * n - 3) * (2 * n - 5)) 15 continue 10 continue c do 20 n = 2, 13 shmit(n,1) = dble(2 * n - 3) * shmit(n-1,1) / dble(n - 1) jj = 2 do 25 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 25 continue 20 continue c f1 = dble(lg(1,1)) f2 = dble(lgt(1,1)) c do 30 n = 1, kmax do 35 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 35 continue 30 continue c 40 continue c if (year .eq. tmold) go to 55 tmold = year t = year - tzero c do 45 n = 1, kmax do 50 m = 1, kmax g(n,m) = gg(n,m) + t * ggt(n,m) 50 continue 45 continue c c --------------------------------------------------------------- c*FON Calcul des coordonnees du dipole boreal et du moment magnetique c*FON (em) - les calculs commencent ici c --------------------------------------------------------------- c 55 continue c ct = cos(thet) st = sin(thet) cph = cos(phi) sph = sin(phi) sp(2) = sph cp(2) = cph c do 60 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) 60 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 = - (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 do 65 n = 3, kmax ar = aor * ar do 64 m = 1, n if (m .eq. n) go to 75 if (const(n,m) .eq. 0.0d0) go to 70 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) go to 80 c 70 continue c p(n,m) = ct * p(n-1,m) dp(n,m) = ct * dp(n-1,m) - st * p(n-1,m) go to 80 c 75 continue 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) c 80 continue c par = p(n,m) * ar if (m .eq. 1) go to 85 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 go to 90 c 85 continue c temp = g(n,m) c 90 continue c br = br - temp * fn(n) * par bt = bt + temp * dp(n,m) * ar c 64 continue 65 continue c c --------------------------------------- c*FON Expression du champ magnetique en gauss c --------------------------------------- c bp = bp / st / 100000.d0 br = br / 100000.d0 bt = bt / 100000.d0 bb = sqrt(br * br + bt * bt + bp * bp) c c **************** c Fin de programme c **************** c return end