subroutine gussen (tgls,xinvla,iguss,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 : Frontieres et regions
c*ROL         Calcul de la presence d'un satellite dans la zone diffuse.
c*ROL         Calcul de la frontiere sud de la zone aurorale diffuse
c*ROL         d'apres le modele de Gussenhoven pour l'activite
c*ROL         magnetique Kp = 3.
c*
c*PAR tgls   (I) : temps geomagnetique local (heures fractionnaires)
c*
c*PAR xinvla (I) : latitude invariante (radians)
c*
c*PAR iguss  (O) : indicateur de presence dans la zone diffuse
c*
c*PAR ier    (O) : code de retour
c*
c*NOT iguss      : 1 = le satellite est dans la zone diffuse
c*NOT iguss      : 0 = le satellite est en dehors
c*
c*NOT ier        : sans objet
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 tgls
      double precision xinvla
      integer iguss
      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,k
c*LOC i,k : indices de boucles
c
      double precision txinv(24),a1s(2),a2s(2),a3s(2),tthets(2)
      double precision t,alfas,xinvli,xoval,xthets,xinvm,tginf,tgsup
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*     Donnees de latitude invariante pour Kp=3
c
      data (txinv(i),i=1,24)/
     >   60.13d0, 60.45d0, 61.39d0, 62.33d0, 63.26d0, 62.19d0,
     >   62.50d0, 63.17d0, 63.69d0, 64.43d0, 65.27d0, 66.35d0,
     >   66.88d0, 67.41d0, 67.94d0, 68.47d0, 67.76d0, 67.17d0,
     >   65.98d0, 64.91d0, 63.73d0, 63.02d0, 62.56d0, 61.59d0/
c
c*    Donnees de la frontiere sud de l'oval auroral
c
      data (a1s(i),i=1,2) /17.36d0, 23.18d0/
      data (a2s(i),i=1,2) / 3.03d0,  4.85d0/
      data (a3s(i),i=1,2) /  3.3d0,   3.3d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
c     ----------------------------------------------------------
c*FON Calcul de la latitude geom.frontiere sud de l oval auroral
c     ----------------------------------------------------------
c
      t = tgls * 15.d0
      t = mod(t,360.d0)
c
      do 10 k = 1, 2
         alfas     = t + a3s(k)
         tthets(k) = a1s(k) + a2s(k) * cos(alfas * rad)
10    continue
c
      xthets = (tthets(1) + tthets(2)) * rad / 2.d0
      xoval  = pid - xthets
c
c     ------------------------------------------
c*FON Calcul de l appartenance a la zone diffuse
c     ------------------------------------------
c
      iguss = 0
c
      xinvm = 60.13d0 * rad
      if (xinvla .ge. xinvm) then
c
         do 15 i = 1, 24
c
            tginf = dble(i - 1)
            tgsup = dble(i)
c
            if (tgls .ge. tginf .and. tgls .le. tgsup) then
c
               xinvli = txinv(i) * rad
c
               if (xinvla .gt. xinvli .and. xinvla .lt. xoval) then
                  iguss = 1
               endif
c
            endif
c
15       continue
c
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end