C=======================================================================
       subroutine GEOVIS(lat,lon,nbp,azi,zen,crdu)
C
C=============================================[ documentation       ]===
C---------------------------------------------[ identification      ]---
C     auteur            : CMS,   mise a jour le 20/11/96
C---------------------------------------------[ role                ]---
c      calcule angles azimutaux et zenithaux du satellite pour n points
C---------------------------------------------[ domaine de validite ]---
C      le common pog doit avoir ete initialise par la subroutine defgeo
C---------------------------------------------[ modifications       ]---
C=============================================[ declarations        ]===
      implicit none
C---------------------------------------------[ parameters          ]---
#include "p_earth.h"
	double precision pi
	parameter ( pi = 3.14159265358979311600d0)
C---------------------------------------------[ arguments           ]---
      integer nbp		! nombre de points
      double precision lat(nbp)	! latitude en radians [-pi/2, pi/2] >0 au nord
      double precision lon(nbp)	! longitude en radians [-pi, pi]    >0 a l'est
      double precision azi(nbp)	! azimut en rad. [-pi, pi] 0 au nord, >0 a l'est
      double precision zen(nbp)	! angle zenithal en radians [-pi/2, pi/2]
      integer crdu       	! compte-rendu
c                         	  crdu=0  OK
c                          	  crdu=-1 common pog pas initialise
C---------------------------------------------[ commons             ]---
#include "c_pog.h"     
C---------------------------------------------[ variables locales   ]---
       double precision slats,clats,slatp,clatp,clon,slon
       double precision rayterre,geocentp
       double precision x,y,z,dsp
       integer i
c---------------------------------------------[ equivalences        ]---
c---------------------------------------------[ fonctions appelees  ]---
       double precision datan3
c---------------------------------------------[ data                ]---
c=============================================[ instructions        ]===
c
c      test initialisation common geosat
c
       if (rsat.eq.0) then
         crdu=-1
       else
c
c      calcul angles
c
          slats= dsin(rflat)
          clats= dcos(rflat)

          do 100 i=1,nbp

            clatp=dcos(lat(i))
            slatp=dsin(lat(i))
            clon = dcos(satlon-lon(i))
            slon = dsin(satlon-lon(i))

            geocentp=datan(rp2_re2*dtan(lat(i)))
            rayterre=rpo/dsqrt(1.d0-exct2*dcos(geocentp)**2)

c      coordonnees du vecteur point terrestre-satellite dans le repere
c      local du point    x sud  y est  z verticale
            x = rsat*(slatp * clats * clon - clatp * slats )
     >          - rayterre * dsin(lat(i)-geocentp)
            y = rsat * clats * slon
            z = rsat*(clatp * clats * clon + slatp * slats )
     >          - rayterre * dcos(lat(i)-geocentp)

c      calcul distance satellite point-vise, azimut et angle zenithal
            dsp = dsqrt( x*x + y*y + z*z )
            azi(i) = pi - datan3( y , x )
            if(azi(i).gt.pi) azi(i) = azi(i) -2*pi
            zen(i) = dacos( z / dsp )
100       continue

          crdu=0
      endif

      return
c=============================================[ formats             ]===
      end