subroutine bwshff (xgsm,ygsm,zgsm,ibwsh,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.01 - 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 position du satellite par rapport a l'onde
c*ROL         de choc Fairfield.
c*
c*PAR xgsm  (I) : coordonnee solaire magnetospherique x du satellite
c*PAR           : (rayons terrestres)
c*PAR ygsm  (I) : coordonnee solaire magnetospherique y du satellite
c*PAR           : (rayons terrestres)
c*PAR zgsm  (I) : coordonnee solaire magnetospherique z du satellite
c*PAR           : (rayons terrestres)
c*
c*PAR ibwsh (O) : indicateur de positionnement du satellite par
c*PAR           : rapport a l'onde de Fairfield
c*
c*PAR ier   (O) : code de retour
c*
c*NOT ibwsh     : 0 = en deca de l'onde de choc
c*NOT ibwsh     : 1 = dans l'onde de choc
c*NOT ibwsh     : 2 = au dela de l'onde de choc
c*
c*NOT On declare le satellite dans le vent solaire a priori. Si
c*NOT xgsm est inferieur ou egal a 15 rayons geomagnetiques on
c*NOT teste s'il appartient a l'onde de choc et en deca de l'onde.
c*
c*NOT ier       : sans objet
c*
c*INF utilise   : aberrm
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.01 - 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 xgsm, ygsm, zgsm
      integer ibwsh
      integer ier
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      integer ier1
c*LOC ier1 : code retour des modules appeles
c
      double precision ga,gb,gc,gd,ge
c*LOC ga,gb,gc,gd,ge : parametre decrivant l'expression quadratique
c*LOC                : de la frontiere dans un repere annulant le terme xy
c
      double precision xgsa,ygsa,zgsa
c*LOC xgsa,ygsa,zgsa : coordonnees solaires magnetospheriques x, y et z 
c*LOC                : corrigees
c
      double precision delta
c*LOC delta : delta = 1 rayon geocentrique terrestre
c
      double precision xin,xout,yin,yout,zin,zout
c*LOC xin,xout,yin,yout,zin,zout : variables pour calcul des champs
c
      double precision bshin,bshout
c*LOC bshin,bshout : champ a la position satellite - delta, et + delta
c
      double precision x2,y2,z2
c*LOC x2,y2,z2 : coordonnees cartesiennes intermediaires
c
      double precision ro,ro2
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
      data ga,gb,gc,gd,ge
     > /-0.038311d0, 1.000211d0, -0.62933d0, 45.657607d0, -652.1d0/
c
      data delta /1.d0/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier  = 0
      ier1 = 0
c
      ibwsh = 2
c
c     -----------------------------------------------------------
c*FON Calculs effectues uniquement si le satellite est en deca ou
c*FON Dans l'onde de choc de Fairfield
c     -----------------------------------------------------------
c
      if (xgsm .le. 15.0d0) then
c
c        --------------------------------------------------------------
c*FON    Prise en compte de l'aberration due au deplacement de la terre
c*FON    l'angle d'aberration est fixe a 4. degres
c        --------------------------------------------------------------
c
         call aberrm(xgsm,ygsm,zgsm,xgsa,ygsa,zgsa,ier1)
c
c        -----------------------------------------------
c*FON    Calcul du champ a la position satellite - delta
c        -----------------------------------------------
c
         xin = xgsa - delta
         yin = ygsa - delta
         zin = zgsa - delta
c
         x2  = xin * xin
         y2  = yin * yin
         z2  = zin * zin
c
         ro2 = y2 + z2
         ro  = sqrt(ro2)
c
         bshin = ga * x2 + gb * ro2 + gc * ro + gd * xin + ge
c
c        -----------------------------------------------
c*FON    Calcul du champ a la position satellite + delta
c        -----------------------------------------------
c
         xout = xgsa + delta
         yout = ygsa + delta
         zout = zgsa + delta
c
         x2  = xout * xout
         y2  = yout * yout
         z2  = zout * zout
c   
         ro2 = y2 + z2
         ro  = sqrt(ro2)
c
         bshout = ga * x2 + gb * ro2 + gc * ro + gd * xout + ge
c
c        ---------------------------------------------------
c*FON    Positionnement du satellite par rapport a l'onde de
c*FON    choc de Fairfield
c        ---------------------------------------------------
c
         if (bshout .lt. 0.d0) then
            ibwsh = 0
         else if (bshin .le. 0.d0 .and. bshout .ge. 0.d0) then
            ibwsh = 1
         endif
c
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end