subroutine pconjr (magin,year,magout,indgm,indval,tilt,rb,rggsm, > rgsmg,tetdip,phidip,dir,rre,theto,phio,rfin, > rmax,np,tr,tthet,tphi,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.24- 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 du point conjugue d'un point donne en tenant c*ROL compte d'une combinaison de champs internes et externes. c* c*PAR magin (I) : type de champ interne c* c*PAR year (I) : annee fractionnaire de calcul c* c*PAR magout (I) : type de champ externe c* c*PAR indgm (I) : type d'indice geomagnetique c*PAR indval (I) : indice geomagnetique pour Kp ou Ae c* c*PAR tilt (I) : angle de tilt (r adians) c* c*PAR rb (I) : distance subsolaire de la magnetopause c*PAR : (rayons terrestres) rb = 11. * 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 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 force c* c*PAR rre (I) : distance radiale geocentrique du point de depart c*PAR : (rayons terrestres) c* c*PAR theto (I) : colatitude geocentrique du point de depart (radians) c*PAR phio (I) : longitude geocentrique du point de depart (radians) c* c*PAR rfin (I) : distance geocentrique d'arret (rayons terrestres) 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 force c*PAR tthet (O) : tableau des colatitudes des points calcules de la c*PAR : ligne de force c*PAR tphi (O) : tableau des longitudes des points calcules de la c*PAR : ligne de force c* c*PAR ier (O) : code de retour c* c*NOT magin : 1 = dipole 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 ou vers l hemisphere oppose) c*NOT : -1 = trace vers les altitudes plus basses c*NOT : (vers la surface du meme hemisphere) c* c*NOT rfin : point d arret de la ligne de force choisi c*NOT : pour la surface rfin = 1 c*NOT rmax : point d arret des lignes de force 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 ier : 1 = l > 500 ou r > rmax ou (rre - 1) <= 0.01 c*NOT ier : 2 = sortie de magnetopause c* c*NOT common : iter c*NOT common : util c* c*INF utilise : magtot, iterat, mpause 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 - utilisation du champ 2005 c*HST version 4.0 - 13.10.10 - introduction du champ igrf10 et passage c*HST de igrf05 à dgrf05 c*HST version 5.0 - 2017.02.24 - introduction du champ igrf15 et passage c*HST de igrf10 à dgrf10 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 double precision rb double precision rggsm(3,3), rgsmg(3,3) double precision tetdip, phidip, dir double precision rre, theto, phio double precision rfin double precision rmax integer np double precision tr(500) double precision tthet(500) double precision tphi(500) 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 ipop c*LOC ipop : indicateur de positionnement par rapport a la magnetosphere c* integer l c*LOC l : pas d'integration c integer ier1,ier2,ier3 c*LOC ier1,ier2,ier : codes retour des modules appeles c double precision r,thet,phi,rp,thetp,phip c*LOC r,thet,phi,rp,thetp,phip : coordonnees spheriques c double precision phipp,pphi,pphid,pthet,rpp,thetpp c*LOC phipp,pphi,pphid,pthet,rpp,thetpp : coordonnees spheriques c double precision cglond,singml c*LOC cglond,singml : variables intermediaires pour la determination du ds c double precision a1,a2,a3,a4,a5,a6,a7 c*LOC a1,a2,a3,a4,a5,a6,a7 : termes des equations c double precision altgc,epsr,epsang,ds0 c*LOC Variables de travail intermediaires c c ----------------------------------------- c*FON Declaration de la fonction interne pinter c ----------------------------------------- c double precision pinter c SAVE c c ---------------------------------------- c*FON Definition de la fonction interne pinter c ---------------------------------------- c pinter(a1,a2,a3,a4,a5,a6,a7) = > ((a2 - a3) * (a7 - a2) * (a7 - a3) * a4 > - (a1 - a3) * (a7 - a1) * (a7 - a3) * a5 > + (a1 - a2) * (a7 - a2) * (a7 - a1) * a6) > / ((a1 - a2) * (a2 - a3) * (a1 - a3)) c c --------------------------------- c*FON Affectation identificateur rcs_id c --------------------------------- c data rcs_id /" >$id$"/ c c -------------------------- c*FON Affectation des constantes c -------------------------- c data epsr/0.001d0/ data epsang/0.00174532d0/ c c ****************** c Debut de programme c ****************** c c --------------- c*FON Initialisations c --------------- c ier = 0 ier1 = 0 ier2 = 0 ier3 = 0 ipop = 0 c sgn = dir c altgc = (rre - 1.0d0) if (theto .lt. epsang) then theto = epsang endif c c ------------------- c*FON Determination du ds c ------------------- c cglond = cos(phio - phidip) singml = cos(tetdip) * cos(theto) + sin(tetdip) > * sin(theto) * cglond ds0 = .06d0 * rre / (1.d0 - singml**2) - 0.058d0 if (ds0 .gt. 0.3d0) ds0 = 0.3d0 if (dir .lt.0. 0d0) then if (altgc .gt. epsr .and. altgc .lt. 2.0d0) then ds0 = altgc / 20.0d0 elseif (altgc .le. epsr) then ier = 0 np = 1 tr(1) = rre tthet(1) = theto tphi(1) = phio goto 30 endif endif ds = ds0 c r = rre phi = phio thet = theto l = 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) call magtot(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,r, > thet,phi,br,bt,bp,b,ier1) c c ---------------------------------- c*FON Memorisation du champ a l'equateur c ---------------------------------- c if (l .eq. 0) then sgn = sign(1.d0,br * dir) endif c l = l + 1 c thetpp = thetp phipp = phip rpp = rp c c ---------------------------------------- c*FON Integration le long de la ligne de champ c ---------------------------------------- c if (l .gt. 30 .and. tr(l) .lt. 2.d0) then ds = ds0 / 2.d0 endif c y(1) = r y(2) = thet y(3) = phi yold(1) = rp yold(2) = thetp yold(3) = phip 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 .ge. 500 .or. rp .gt. rmax) then ier = 1 ipop = 1 endif 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 = 2 go to 20 endif c tr(l) = 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 tthet(l) = thetp tphi(l) = phip c if (r .gt. rfin .and. ipop .eq. 0) go to 10 c c ------------- c*FON Fin de boucle c ------------- c pthet = pinter(rpp,rp,r,thetpp,thetp,thet,rfin) pphi = pinter(rpp,rp,r,phipp,phip,phi,rfin) pphi = pphi + dpi pphi = mod(pphi,dpi) np = l + 1 tr(np) = rfin tthet(np) = pthet tphi(np) = pphi c 20 continue c c* ----------------------------------------------------- c*FON Initialisation des parametes en sortie a des valeurs c*FON arbitraires en cas de non calcul de la ligne de force c* ----------------------------------------------------- c if (ipop .eq. 1 .or. ipop .eq. 2) then c np = l pthet = 999.0d0 * rad pphid = 999.0d0 * rad c tr(np) = 999.0d0 tthet(np) = pthet tphi(np) = pphid c endif c 30 continue c c **************** c Fin de programme c **************** c return end