subroutine mcilwe (req,tgmleq,phikv,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.06.05 - V 2.0 c*VER 03.01.06 - V 2.1 c* c*AUT spec. CNES - JC KOSIK - janvier 1991 c*AUT port. CISI c* c*ROL Theme : Calculs de geophysique c*ROL Calcul du potentiel electrique en kilovolts en c*ROL fonction du rayon equatorial et de l'heure locale. c* c*PAR req (I) : distance radiale geocentrique dans le plan equatorial c* : (rayons terrestres) c*PAR tgmleq (I) : temps geomagnetique local (heures fractionnaires) c* c*PAR phikv (O) : potentiel electrique (kilovolts) c* c*PAR ier (O) : code de retour c* c*NOT ier : 0 = OK c*NOT ier : 1 = tgmleq < 0.0d0 c*NOT ier : 2 = req < 4.0d0 (rayons terrestres) c* c*NOT common : util 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 req, tgmleq double precision phikv 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 c --------------------------------- c*FON Declaration des variables locales c --------------------------------- c integer i,j c*LOC i,j : indices de boucles c double precision tga(6,20),tphi(20),tgc(20),tgb(6),td(6),tpa(6) double precision tpb(20) double precision beqnt,phikv0,som,expos1,expos2,expos,phi,cphi double precision r3,den,hrad 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 c* Constantes pour le calcul du potentiel c data ((tga(i,j),i=1,6),j=1,20)/ > 2.8d0, 5.4d0, 0.6d0, 2.9d0, -1.2d0, 0.6d0, > 6.0d0, -1.7d0, 1.7d0, -1.1d0, 0.9d0, -0.2d0, > -6.5d0, 3.2d0, -1.2d0, 1.9d0, -1.1d0, 0.4d0, > 5.7d0, -2.5d0, 1.1d0, -0.9d0, 0.5d0, -0.2d0, > -1.4d0, 1.6d0, -2.2d0, 1.3d0, -1.2d0, -0.1d0, > 4.5d0, -3.3d0, 1.0d0, -2.1d0, 0.5d0, -0.7d0, > -5.3d0, 0.7d0, -3.2d0, 0.6d0, -1.5d0, -0.1d0, > 3.6d0, -3.8d0, 0.0d0, -1.5d0, -0.2d0, -0.7d0, > -3.1d0, -2.0d0, -2.5d0, -0.6d0, -1.3d0, -0.4d0, > 1.7d0, -1.5d0, 0.2d0, -0.9d0, -0.3d0, -0.2d0, > 1.0d0, -1.5d0, -1.1d0, -0.4d0, -0.6d0, -0.2d0, > 1.7d0, -0.5d0, -0.5d0, -0.6d0, 0.3d0, -0.2d0, > 2.9d0, -2.1d0, 1.7d0, -1.2d0, 0.4d0, -0.2d0, > -0.1d0, 3.6d0, -4.9d0, 3.5d0, -2.2d0, 0.7d0, > 2.7d0, -3.4d0, 5.5d0, -0.6d0, 0.6d0, -0.5d0, > 0.7d0, 5.7d0, 2.1d0, 0.6d0, 0.1d0, -0.2d0, > 6.0d0, 2.3d0, 2.3d0, 0.3d0, 0.1d0, -0.1d0, > 3.5d0, 0.3d0, 2.6d0, -0.8d0, 0.7d0, -0.4d0, > 9.5d0, 5.1d0, 1.8d0, 1.7d0, -0.7d0, 0.3d0, > 3.0d0, -1.9d0, 2.1d0, -1.3d0, 0.9d0, -0.4d0/ c data tphi/ > 4.d0, 6.d0, 8.d0, 10.d0, 12.d0, 14.d0, 16.d0, 18.d0, > 20.d0, 21.d0, 22.d0, 22.5d0, 23.d0, 23.5d0, 0.d0, 0.5d0, > 1.d0, 1.5d0, 2.d0, 3.d0/ c data tgc /9*2.d0, 2*1.d0, 7*0.5d0, 2*1.d0/ c data tgb /0.d0, 40.d0, 100.d0, 180.d0, 280.d0, 400.d0/ c data td /30.d0, 50.d0, 70.d0, 90.d0, 110.d0, 130.d0/ c c ****************** c Debut de programme c ****************** c ier = 0 c c ----------------------------- c*FON Calcul du champ en nanoteslas c ----------------------------- c phi = tgmleq * dpi / 24.d0 cphi = cos(phi) phikv = 999.0d0 c if (req .ge. 4.0d0) then c r3 = req * req * req den = 1.d0 + 1728.d0 / r3 beqnt = 6.d0 - 24.d0 * cphi + 18.d0 * cphi * cphi / > den + 31000.d0 / r3 c if (beqnt .ge. 0.0d0) then c c --------------------- c*FON Calcul du terme isole c --------------------- c phikv0 = 10.d0 - 92.d0 * (beqnt / 31000.d0)**(1.d0 / 3.d0) c c ------------------------------ c*FON Calcul de sommation sur i et j c ------------------------------ c hrad = pi / 12.d0 som = 0.d0 c do 10 i = 1, 6 tpa(i) = log(2.d0) / (td(i) * td(i)) do 10 j = 1, 20 tpb(j) = log(2.d0) / (1.d0 - cos(tgc(j) * hrad)) expos1 = -tpa(i) * (beqnt - tgb(i))**2 expos2 = -tpb(j) * > (1.d0 - cos((tgmleq - tphi(j)) * hrad)) expos = expos1 + expos2 som = som + tga(i,j) * exp(expos) 10 continue c c -------------------- c*FON Potentiel electrique c -------------------- c phikv = phikv0 + som ier = 0 c elseif (beqnt .lt. 0.0d0) then ier = 1 endif else ier = 2 endif c c **************** c Fin de programme c **************** c return end