subroutine calphi (xgse,ygse,zgse,rgsa,thetd,phid,
     >                   ival,phi1,phi2,ier)
c*
c***********************************************************************
c*
c*      "Copyright [c] CNES 98 - tous droits reserves"
c*      **********************************************
c*
c*PRO MAGLIB
c*
c*VER 01.10.23 - V 2.0
c*VER 03.01.06 - V 2.1
c*
c*AUT spec. CNES - JC KOSIK - octobre 2001
c*AUT port. CISI
c*
c*ROL Theme : Astronomie et calculs d'orbite
c*ROL       Recherche des deux fuseaux qui encadrent le point donne.
c*ROL       Ces fuseaux sont donnes par les angles phi1 et phi2
c*ROL       comptes a partir de l'axe ygsa.
c*
c*PAR xgse  (I) : composante solaire ecliptque en x (degres)
c*PAR ygse  (I) : composante solaire ecliptque en y (degres)
c*PAR zgse  (I) : composante solaire ecliptque en z (degres)
c*
c*PAR rgsa  (I) : composante solaire corrigee 
c
c*PAR thetd (I) : colatitude geocentrique du dipole (degres)
c*PAR phidp (I) : longitude geocentrique du dipole (degres)

c*PAR ival  (I) : numero du fuseau
c*PAR phi1  (O) : angle phi1 du premier meridien (degres)
c*PAR phi2  (O) : angle phi2 du second meridien (degres)
c*
c*PAR ier   (O) : code retour
c*
c*NOT ier      : 0 = OK
c*NOT ier      : 1 = thetd superieur a 165 degres
c*
c*NOT common   : util
c*
c*INF utilise  : aberrm
c*
c*HST version 2.0 - 01.10.23 - Enrichissement de la maglib au CDPP
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
      integer ival, ier
      double precision xgse,ygse,zgse,rgsa,thetd,phid,phi1,phi2
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
      double precision thetmax
c*LOC thetmax : valeur maximale de l'angle theta
c
      double precision phiinf, phisup
c*LOC phiinf, phisup : bornes inferieure et superieure de l'angle phi
c
      double precision xgsa,ygsa,zgsa
c*LOC xgsa,ygsa,zgsa : composantes solaires ecliptiques corrigees
c
      integer i
c*LOC i :indice de boucles
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
      data thetmax/165.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier = 0
c
c     ------------------------------------------------------
c*FON Calcul des coordonnees dans le systeme du vent solaire
c*FON i.e. avec aberration
c     ------------------------------------------------------
c
      call aberrm (xgse,ygse,zgse,xgsa,ygsa,zgsa,ier)
c
c     ---------------------------------------------------------------
c*FON Calcul des coordonnees thet phi du systeme Olson
c*FON thet est compte en degres depuis l'axe xgsa, phi est
c*FON compte positivementdepuis l'axe ygsa de 0 a 90 vers
c*FON l'axe zgsa positif et negativement vers l'axe zgsa
c*FON negatif. Pour ygsa negatif on fait ygsa = -ygsa
c*FON (symetrie par rapport a l'axe zgsa)
c     ---------------------------------------------------------------
c
c     changement de signe de ygsa si ygsa est negatif
c
      if (ygsa .lt. 0.d0) then
         ygsa = - ygsa
      endif
c
      rgsa = dsqrt(xgsa*xgsa + ygsa*ygsa + zgsa*zgsa)
c
      thetd = dacos(xgsa / rgsa) * deg
c
      phid  = datan(zgsa / ygsa) * deg
c
      if (thetd .le. thetmax) then
c
         ier = 1
c
         do 20  i = 1, 12

            phiinf = -90.d0  + dble(i-1) * 15.d0
            phisup = -90.d0  + dble(i)   * 15.d0

            if (phid .ge. phiinf .and. phid .le. phisup) then
c
               ival  = i
               ier   = 0
               phi1  = phiinf
               phi2  = phisup
c
            endif
c
20       continue
      endif
c
c     ****************
c     Fin de programme
c     ****************
c

      return
      end