subroutine posmag (magout,isatex,year,rre,thetr,phir,alfag,alfas,
     >                   deltas,tilt,rgsm,rggsm,rgsmg,rgdip,rgse,
     >                   tetdip,phidip,xgsm,ygsm,zgsm,xgse,ygse,zgse,
     >                   tgl,flg,xlamb,tglc,hsl,clatgm,clongm,
     >                   iposmg,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*
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI
c*
c*ROL Theme : Frontieres et regions
c*ROL         Calcul de la position d'un satellite dans la magnetosphere
c*ROL         par rapport aux differentres regions.
c*ROL         Calcul aussi des parametres essentiels tgl, L et mlt,
c*ROL         et des coordonnees dans divers systemes d axes.
c*ROL         L'indicateur de positionnement iposmg est renseigne.
c*
c*PAR magout (I) : type de champ magnetique externe
c*
c*PAR isatex (I) : type de satellite
c*
c*PAR year   (I) : annee fractionnaire >= 2000
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 alfag  (I) : ascension droite de Greenwich (radians)
c*PAR alfas  (I) : ascension droite du soleil (radians)
c*PAR deltas (I) : declinaison du soleil (radians)
c*
c*PAR tilt   (I) : angle de tilt (radians)
c*
c*PAR rgsm   (I) : matrice de passage du repere geographique au
c*PAR            : repere solaire magnetique
c*
c*PAR rggsm  (I) : matrice de passage du repere geographique au
c*PAR            : repere solaire magnetospherique
c*PAR rgsmg  (I) : matrice de passage du repere magnetospherique
c*PAR            : au repere geographique
c*
c*PAR rgdip  (I) : matrice (3,3) de passage du repere geographique
c*PAR            : au repere dipolaire
c*
c*PAR rgse   (I) : matrice de passage du repere geographique
c*PAR            : au repere solaire ecliptique
c*
c*PAR tetdip (I) : colatitude geocentrique du dipole (radians)
c*PAR phidip (I) : longitude geocentrique du dipole (radians)
c*
c*PAR xgsm   (O) : coordonnee solaire magnetospherique en x (rayons terrestres)
c*PAR ygsm   (O) : coordonnee solaire magnetospherique en y (rayons terrestres)
c*PAR zgsm   (O) : coordonnee solaire magnetospherique en z (rayons terrestres)
c*
c*PAR xgse   (O) : coordonnee solaire ecliptique en x (rayons terrestres)
c*PAR ygse   (O) : coordonnee solaire ecliptique en y (rayons terrestres)
c*PAR zgse   (O) : coordonnee solaire ecliptique en z (rayons terrestres)
c*
c*PAR tgl    (O) : temps geomagnetique local du satellite
c*               : (heures fractionnaires)
c*PAR flg    (O) : parametre L de Galperin
c*PAR xlamb  (O) : latitude invariante (radians)
c*
c*PAR tglc   (O) : temps geomagnetique local du point conjugue
c*               : (heures fractionnaires)
c*
c*PAR hsl    (O) : hauteur de l'ombre (kilometres)
c*
c*PAR clatgm (O) : latitude geomagnetique  (radians)
c*PAR clongm (O) : longitude geomagnetique (radians)
c*PAR iposmg (O) : tableau (dim=15) des indicateurs de positions
c*
c*PAR ier    (O) : code de retour
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 isatex     : 0 = satellite proche      (apogee < 8. r. terre.)
c*NOT isatex     : 1 = satellite excentrique (apogee > 8. r. terre.)
c*
c*NOT iposmg     : 1 = le satellite appartient a la region
c*NOT iposmg     : 0 = le satellite n'appartient pas a la region
c*
c*NOT iposmg(1)  : calotte polaire
c*NOT iposmg(2)  : ovale auroral
c*NOT iposmg(3)  : cusp
c*NOT iposmg(4)  : zone auroral diffuse
c*NOT iposmg(5)  : ceinture de Van Allen pour xbelt = 6. 
c*NOT iposmg(6)  : plasmasphere
c*NOT iposmg(7)  : couche neutre
c*NOT iposmg(8)  : couche de plasma
c*NOT iposmg(9)  : magnetosphere
c*NOT iposmg(10) : magnetopause
c*NOT iposmg(11) : magnetogaine
c*NOT iposmg(12) : onde de choc
c*NOT iposmg(13) : vent solaire
c*NOT iposmg(14) : ombre de la terre (satellite auroral)
c*NOT iposmg(15) : ceinture de Van Allen pour xbelt = 3.5  
c*
c*NOT ier        : sans objet
c*
c*NOT common     : util
c*
c*INF utilise    : spcar, geogsm, tgml, geose, bwshff, mpsib, msheath
c*INF utilise    : dconjr, invar, invlat, ggeom, oval, calpol, cusp
c*INF utilise    : gussen, rbelt, chapel, cahsl, posnsh, pospsh 
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.05.30 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*HST version 3.0 - 01.01.07 - utilisation du champ IGRF 2005
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 magout, isatex
      double precision year
      double precision rre, thetr, phir
      double precision alfag, alfas, deltas, tilt
      double precision rgsm(3,3), rggsm(3,3), rgsmg(3,3)
      double precision rgdip(3,3), rgse(3,3)
      double precision tetdip, phidip
      double precision tr(500), tthet(500), tphi(500)
      double precision xgsm, ygsm, zgsm
      double precision xgse, ygse, zgse
      double precision tgl, flg, xlamb, tglc, hsl
      double precision clatgm, clongm
      integer iposmg(15)
      double precision znsh
      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 imagn,ibwsh,imp,isheath,isolw,ioval,ipol
      integer icusp,iguss,iallen,ichapp,ihsl,insh,ipsh
c*LOC Indicateurs d'appartenance aux zones
c
      integer i
c*LOC i : indice de boucle
c
      integer magin
c*LOC magin : type de champ magnetique interne (= 3 : IGRF 2005)
c
      integer indgm
c*LOC indgm : type d'indice geomagnetique (=1 : Kp)
c
      integer indval
c*LOC indval : indice geomagnetique : niveau d'amplitude du champ (1 a 6)
c
      integer np
c*LOC np : nombre de points calcules
c
      integer ier1,ier2,ier3,ier4,ier5,ier6,ier7,ier9
      integer ier10,ier11,ier12,ier13,ier14,ier15,ier16
      integer ier17,ier18,ier19
c*LOC ier1 a ier19 : codes retour des modules appeles
c
      double precision rmagc,thetc,phic
c*LOC rmagc,thetc,phic : coordonnees spheriques 
c*
      double precision xmag,ymag,zmag
c*LOC xmag,ymag,zmag : composantes cartesiennes en x, y et z
c
      double precision xbelt
c*LOC xbelt : valeur de la limite externe de la ceinture
c
      double precision rfin,rmax
c*LOC rfin,rmax : distances geocentriques d'arret et maximum autorisee
c
      double precision dir
c*LOC dir : direction du trace des lignes de champ
c
      double precision rb
c*LOC rb : distance subsolaire de la magnetopause
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier   = 0
      ier1  = 0
      ier2  = 0
      ier3  = 0
      ier4  = 0
      ier5  = 0
      ier6  = 0
      ier7  = 0
      ier9  = 0
      ier10 = 0
      ier11 = 0
      ier12 = 0
      ier13 = 0
      ier14 = 0
      ier15 = 0
      ier16 = 0
      ier17 = 0
      ier18 = 0
      ier19 = 0
c
c     --------------------------------------------------
c*FON Mise a 0 du tableau d'indicateur de positionnement
c*FON Initialisations de valeurs par defaut
c     --------------------------------------------------
c
      do 10 i = 1, 15
         iposmg(i) = 0
10    continue
c
c     dans tous les cas mais surtout pour satellites auroraux :
      iposmg(9) = 1
c
      tgl    = 999.0d0
      flg    = 999.0d0
      xlamb  = 999.0d0 * rad
      tglc   = 999.0d0
      clatgm = 999.0d0 * rad
      clongm = 999.0d0 * rad
c
c     -------------------------------------------------
c*FON Calcul des coordonnees solaires magnetospheriques
c     -------------------------------------------------
c
      call spcar(rre,thetr,phir,xmag,ymag,zmag,ier1)
c
      call geogsm(rggsm,xmag,ymag,zmag,xgsm,ygsm,zgsm,ier2)
c
c     ---------------------------------------------------
c*FON Calcul du temps geomagnetique local du satellite tgl
c     ---------------------------------------------------
c
      call tgml(rgsm,thetr,phir,tgl,ier3)
c
c     ------------------------------------------
c*FON Calcul des coordonnees solaire-ecliptiques
c     ------------------------------------------
c
      call geose(rgse,xmag,ymag,zmag,xgse,ygse,zgse,ier4)
c
c     --------------------------------------------------------
c*FON Si nosat est superieur a 2, verification si le satellite
c*FON Appartient a la magnetosphere de sibeck
c     --------------------------------------------------------
c
      imagn = 1
c
      if (isatex .eq. 1) then
c
c     ---------------------------------------------------------------
c*FON Calcul de la position par rapport a l onde de choc de Fairfield
c*FON et la magnetopause de Sibeck
c     ---------------------------------------------------------------
c
         call bwshff(xgsm,ygsm,zgsm,ibwsh,ier5)
c
         if (ibwsh .eq. 1) then
            iposmg(12) = ibwsh
         endif
c
         call mpsib(xgsm,ygsm,zgsm,imp,ier6)
c
         if (imp .eq. 1) then
            iposmg(10) = imp
         endif
c
         call msheath (xgsm,ygsm,zgsm,imagn,isheath,isolw,ier7)
c
         iposmg(13) = isolw
         iposmg(11) = isheath
         iposmg(9)  = imagn
c
      endif
c
      if (imagn .eq. 1) then
c
c     ------------------------------------------------------
c*FON Calcul du point conjugue situe dans le meme hemisphere
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 DGRF 2015
c*NOT magin      : 6 = champ IGRF 2020

c     ------------------------------------------------------
c
      if (year .ge. 2000 .and. year .lt. 2005) then
c
         magin  =  2 
c
      else if (year .ge. 2005 .and. year .lt. 2010) then
c
         magin  =  3
c
      else if (year .ge. 2010 .and. year .lt. 2015) then
c
         magin  =  4
c
      else if (year .ge. 2015 .and. year .lt. 2020) then
c
         magin  =  5
c
      else if (year .ge. 2020) then
c
         magin  =  6
c
      else
         ier = 1
      endif
c
c
         indgm  =   1
         indval =   3
         rb     =  11.d0
         rfin   =   1.d0
         rmax   =  40.d0
         dir    =  -1.d0
c
         call dconjr(magin,year,magout,indgm,indval,tilt,rb,rggsm,
     >               rgsmg,dir,rre,thetr,phir,rfin,rmax,np,tr,
     >               tthet,tphi,ier1)
c
         if (ier1 .eq. 0) then
c
            rmagc  = tr(np)
            thetc  = tthet(np)
            phic   = tphi(np)
c   
            call invar(rmagc,thetc,phic,flg,ier2)
c
            call invlat(flg,xlamb,ier3)
c
c     ---------------------------------------
c*FON Calcul d appartenance a l ovale auroral
c     ---------------------------------------
c
c     ----------------------------------------------------------
c*FON Calcul des coordonnees geomagnetiques et mlt (appele tglc)
c*FON du pied de la ligne de champ passant par le satellite
c     ----------------------------------------------------------
c
            call ggeom(rgdip,thetc,phic,clatgm,clongm,ier9)
            call tgml(rgsm,thetc,phic,tglc,ier10)

            call oval(tglc,clatgm,ioval,ier11)
            call calpol(tglc,clatgm,ipol,ier12)
c
c     -------------------------------------------
c*FON Calcul de la position par rapport a la CUSP
c     -------------------------------------------
c
            call cusp(clatgm,tglc,icusp,ier13)
c
c     ------------------------------------------------------
c*FON Calcul de la frontiere sud de la zone aurorale diffuse
c     ------------------------------------------------------
c
            call gussen(tglc,xlamb,iguss,ier14)
c
            iposmg(1) = ipol
            iposmg(2) = ioval
            iposmg(3) = icusp
            iposmg(4) = iguss
c
         endif
c
c     -------------------------------------------------------
c*FON Si le satellite est a moins de 7 rayons terrestres.....
c     -------------------------------------------------------
c
         if (rre .lt. 7.d0) then
c
c     ------------------------------------------------
c*FON Calcul d appartenance a la ceinture de Van Allen
c     ------------------------------------------------
c
c     Ceinture a 3.5 r. terr.
c     -----------------------
c
            xbelt = 3.5d0
            call rbelt(xbelt,tetdip,phidip,rre,thetr,
     +                 phir,iallen,ier15)
            iposmg(15) = iallen
c
c     Ceinture a 6.0 r. terr.
c     -----------------------
c
            xbelt = 6.d0
            call rbelt(xbelt,tetdip,phidip,rre,thetr,
     +                 phir,iallen,ier15)
            iposmg(5) = iallen
c
c     ---------------------------------------------------
c*FON Calcul d appartenance a la plasmasphere de Chappell
c     ---------------------------------------------------
c
            call chapel(tgl,flg,ichapp,ier16)
            iposmg(6) = ichapp
c
c     --------------------------------------------------------
c*FON Calcul de la distance du satellite a l ombre de la terre
c     --------------------------------------------------------
c
            call cahsl(alfag,alfas,deltas,rre,thetr,phir,
     +                 hsl,ihsl,ier17)
            iposmg(14) = ihsl
c
         endif
c
      endif
c
c     -----------------------------------------------------------------
c*FON Calcul de la distance a la couche neutre et a la couche de plasma 
c*FON si satellite excentrique et xgsm inferieur a -9 rt
c     -----------------------------------------------------------------
c
      if (isatex .eq. 1 .and. xgsm .le. -9.d0) then
c
c     ----------------------------------------------------
c*FON Calcul de la position par rapport a la couche neutre
c     ----------------------------------------------------
c
         call posnsh(tilt,xgsm,ygsm,zgsm,insh,znsh,ier18)
         iposmg(7) = insh
c
c     -------------------------------------------------------
c*FON Calcul de la position par rapport a la couche de plasma
c     -------------------------------------------------------
c
         call pospsh(tilt,xgsm,ygsm,zgsm,ipsh,ier19)
         iposmg(8) = ipsh
c
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end