subroutine corgm (magin,year,magout,indgm,indval,tilt,rb, > rggsm,rgsmg,rgdip,tetdip,phidip, > rre,theto,phio,rmax,np, > tr,tthet,tphi,tdr,tdtet,tdphi,ieq,dreq, > dteteq,dphieq,corthet,corlat,corlon,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*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 13.10.10 - V 4.0 c*VER 2017.02.23 - V 5.0 c* c*AUT spec. CNES - JC KOSIK - janvier 1991 c*AUT port. CISI c*AUT adapt. AKKA c* c*ROL Theme : Calculs de geophysique c*ROL Calcul des coordonnees geomagnetiques corrigees, c*ROL latitude et longitude, tenant compte de champ magnetique c*ROL interne (de type DGRF 2000 ou DGRF 2005 ou DGRF 2010 ou IGRF 2015) c*ROL et externe (de type Tsyganenko 1897 et Tsyganenko 1989 et Kosik 97). c* c*PAR magin (I) : type de champ magnetique interne c* c*PAR year (I) : annee fractionnaire c* c*PAR magout (I) : type de champ magnetique externe c* c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae c* c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ c* c*PAR tilt (I) : angle de tilt (radians) c* c*PAR rb (I) : distance subsolaire de la magnetopause c*PAR : (rayons terrestres) rb = 10. * re c* c*PAR rggsm (I) : matrice (3,3) de passage du repere geographique au c*PAR : repere solaire magnetospherique c*PAR rgsmg (I) : matrice (3,3) de passage du repere solaire c*PAR : magnetospherique au repere geographique c* c*PAR rgdip (I) : matrice (3,3) de passage du repere geographique c*PAR : au repere dipolaire c* c*PAR tetdip (I) : colatitude geocentrique du dipole (radians) c*PAR phidip (I) : longitude geocentrique du dipole (radians) c* c*PAR dir (I) : direction du trace des lignes de champ c* c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) c* c*PAR theto (I) : colatitude geocentrique (radians) c*PAR phio (I) : longitude geocentrique (radians) c* c*PAR rmax (I) : distance geocentrique maximum autorisee dans le trace c*PAR : (rayons terrestres) c* c*PAR np (O) : nombre de points calcules c* c*PAR tr (O) : tableau des distances geocentriques des point de la c*PAR : ligne de champ c* c*PAR tthet (O) : tableau des colatitudes des points calcules de la c*PAR : ligne de champ c*PAR tphi (O) : tableau des longitudes des points calcules de la c*PAR : ligne de champ c* c*PAR tdr (I) : tableau des distances dipolaires (r. terrestres) c* c*PAR tdtet (I) : tableau des colatitudes dipolaires (radians) c*PAR tdphi (I) : tableau des longitudes dipolaires (radians) c* c*PAR ieq (O) : indicateur de passage a l'equateur c* c*PAR dreq (O) : distance radiale geocentrique du point equatorial c*PAR : de la ligne de champ (rayons terrestres) c*PAR dteteq (O) : colatitude geocentrique du point equatorial c*PAR : de la ligne de champ (radians) c*PAR dphieq (O) : longitude geocentrique du point equatorial c*PAR : de la ligne de champ (radians) c* c*PAR corthet (O) : colatitude geomagnetique corrigee (radians) c*PAR corlat (O) : latitude geomagnetique corrigee (radians) c*PAR corlon (O) : longitude geomagnetique corrigee (radians) c* c*PAR ier (O) : code de retour c* c*NOT magin : 1 ---> champ dipolaire c*NOT magin : 2 ---> champ DGRF 2000 c*NOT magin : 3 ---> champ DGRF 2005 c*NOT magin : 4 ---> champ DGRF 2010 c*NOT magin : 5 ---> champ IGRF 2015 c* c*NOT year : year doit etre >= 2000. c* c*NOT magout : 0 = pas de champ externe c*NOT magout : 1 = Tsyganenko 87 c*NOT magout : 2 = Tsyganenko 89 c*NOT magout : 3 = Kosik 97 c* c*NOT indgm : 1 ---> indice geomagnetique Kp c*NOT indval : 1 ---> Kp = 0 , 0+ c*NOT indval : 2 ---> Kp = 1- , 1 , 1+ c*NOT indval : 3 ---> Kp = 2- , 2 , 2+ c*NOT indval : 4 ---> Kp = 3- , 3 , 3+ c*NOT indval : 5 ---> Kp = 4- , 4 , 4+ c*NOT indval : 6 ---> Kp > 5- c* c*NOT indgm : 2 ---> indice geomagnetique Ae c*NOT indval : 1 ---> Ae = 0 - 50 c*NOT indval : 2 ---> Ae = 50 - 100 c*NOT indval : 3 ---> Ae = 100 - 150 c*NOT indval : 4 ---> Ae = 150 - 250 c*NOT indval : 5 ---> Ae = 250 - 400 c*NOT indval : 6 ---> Ae >= 400 c* c*NOT dir : +1 = trace vers les altitudes plus elevees c*NOT : (depuis la surface vers l hemisphere oppose) c* c*NOT rfin : point d'arret de la ligne de champ choisi c*NOT : pour la surface rfin=1 c* c*NOT rmax : point d arret des lignes de champ ouvertes c* c*NOT ieq : 1 = passage a l'equateur c*NOT : 0 = non passage a l'equateur c* c*NOT ier : 0 = OK c*NOT : 1 = la ligne de champ ne croise pas le c*NOT : plan equatorial c* c*NOT common : util, iter c* c*INF utilise : magtot, iterat, mpause, equdip, angdip, corr c* c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP c*HST version 2.0 - 01.05.30 - Ajout indice Kosik 97 c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 c*HST version 3.0 - 01.01.07 - adoption des coefficients definitifs c*HST du champ 2000 c*HST version 4.0 - 13.10.10 - Adoption des coefficients definitifs c*HST du champ 2005 et des coefficients provisoires c*HST 2010 c*HST version 5.0 - 2017.02.23 - Adoption des coefficients definitifs c*HST du champ 2010 et des coefficients provisoires c*HST 2015 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 integer magin double precision year integer magout, indgm, indval double precision tilt, rb double precision rggsm(3,3), rgsmg(3,3), rgdip(3,3) double precision tetdip, phidip, dir double precision rre, theto, phio double precision rmax integer np double precision tr(500), tthet(500), tphi(500) double precision tdr(500), tdtet(500), tdphi(500) integer ieq double precision dreq, dteteq, dphieq double precision corlat, corthet, corlon integer ier c c ---------------------------------- c*FON Declaration des variables communes c ---------------------------------- c double precision y(3),yold(3) double precision br,bt,bp,b,st,sgn,ds c*COM y : composantes spheriques c*COM yold : composantes spheriques precedantes c*COM br,bt,bp : composantes radiale, tangentielle et azimuthale du champ c*COM b : module du champ c*COM st : sinus theta c*COM sgn : sens du calcul (+1, -1) c*COM ds : intervalle d'integration le long de la ligne de champ c common/iter/y,yold,br,bt,bp,b,st,sgn,ds c 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 c --------------------------------- c*FON Declaration des variables locales c --------------------------------- c integer l c*LOC l : pas d'integration c integer npeq c*LOC npeq : indice courant du point de la ligne de champ c integer ipop c*LOC ipop : indicateur de positionnement par rapport a la magnetosphere c* integer ier1,ier2,ier3,ier4,ier5,ier6 c*LOC ier1,ier2,ier3,ier4,ier5,ier6 : codes retour des modules appeles c double precision brm,btm,bpm,bm c*LOC brm,btm,bpm,bm : coordonnees spheriques double precision cglond,singml c*LOC cglond,singml : variables intermediaires pour la determination du ds c double precision altgc c*LOC altgc : altitude en km du 1er point de la ligne de champ c double precision epsr c*LOC epsr : epsilon pour validite de l'altitude c double precision r,thet,phi,rp,thetp,phip c*LOC r,thet,phi,rp,thetp,phip : coordonnees spheriques c double precision tetdp,phidp,theti c*LOC tetdp,phidp,theti : coordonnees spheriques c double precision rdp,prod 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 Debut de programme c ****************** c c ------------------------- c*FON Initialisations de depart c ------------------------- c ier = 0 ier1 = 0 ier2 = 0 ier3 = 0 ier4 = 0 ier5 = 0 ier6 = 0 corlat = 999.d0 * rad corlon = 999.d0 * rad c dir = +1.d0 sgn = dir ieq = 0 altgc = (rre - 1.d0) epsr = 0.01d0 c c ------------------- c*FON Determination du ds c ------------------- c cglond = cos(phio - phidip) singml = cos(tetdip) * cos(theto) + sin(tetdip) > * sin(theto) * cglond ds = .06d0 * rre / (1.d0 - singml**2) - 0.058d0 c if (ds .gt. 1.d0) then ds = 0.3d0 endif c if (dir .lt. 0.d0) then if (altgc .gt. epsr .and. altgc .lt. 2.0d0) then ds = altgc / 20.d0 else if (altgc .le. epsr)then ier = 1 endif endif c if (ier .eq. 1) goto 20 c r = rre thet = theto phi = phio l = 0 ieq = 0 c c -------------------------------------- c*FON Boucle de traitement tant que r > rfin c*FON et que = 0 c -------------------------------------- c 10 continue c c --------------------- c*FON Calcul du champ total c --------------------- c st = sin(thet) c call magtot(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,r, > thet,phi,brm,btm,bpm,bm,ier1) c c ---------------------------------------- c*FON Integration le long de la ligne de champ c ---------------------------------------- c if (l .eq. 0) then yold(1) = 0.d0 yold(2) = 0.d0 yold(3) = 0.d0 sgn = sign(1.d0, brm * dir) endif l = l + 1 c br = brm bt = btm bp = bpm b = bm y(1) = r y(2) = thet y(3) = phi call iterat(l,ier2) r = y(1) thet = y(2) phi = y(3) rp = yold(1) thetp = yold(2) phip = yold(3) c if (l .gt. 500 .or. rp .gt. rmax) ier = 1 c if (ier .eq. 1) go to 20 c c ----------------------------------------------- c*FON Calcul et test d'appartenance a la magnetopause c*FON de type Shabansky c ----------------------------------------------- c call mpause(rb,rggsm,rp,thetp,phip,ipop,ier3) c if (ipop .eq. 2) then ier = 1 go to 20 endif c np = l tr(np) = rp c if (thet .gt. pi) then thet = dpi - thet phi = pi + phi endif c if (thetp .gt. pi) then thetp = dpi - thetp phip = pi + phip endif c c ----------------------------------------------------- c*FON Calcul de la colatitude et de la longitude dipolaires c ----------------------------------------------------- c call angdip(rgdip,rp,thetp,phip,rdp,tetdp,phidp,ier4) c tdr(np) = rdp tthet(np) = thetp tdtet(np) = tetdp tphi(np) = phip tdphi(np) = phidp prod = (tdtet(np-1) - pid) * (tdtet(np) - pid) c if (np .gt. 8 .and. prod .lt. 0.d0) then ieq = 1 npeq = np c c -------------------------------------------------- c*FON Calcul de l'intersection de la ligne de champ avec c*FON l'equateur geomagnetique dipolaire incline c -------------------------------------------------- c call equdip(npeq,tdtet,tdphi,tdr,dteteq,dphieq,dreq,ier5) theti = tdtet(1) c c -------------------------------------------------------- c*FON calcul de la colatitude geomagnetique et de la longitude c*FON geomagnetique corrigees c -------------------------------------------------------- c call corr(altgc,theti,dreq,dphieq,corlat,corlon,ier6) c corthet = pid - corlat c endif c if (ieq .eq. 0. and .ier .eq. 0) go to 10 c c ----------------------------------------------------- c*FON Calcul de la colatitude et de la longitude dipolaires c ----------------------------------------------------- c call angdip(rgdip,rp,thetp,phip,rdp,tetdp,phidp,ier4) c tdtet(np) = tetdp tdphi(np) = phidp tdr(np) = rdp 20 continue c c **************** c Fin de programme c **************** c 999 continue return end