!
! $Header$
!
SUBROUTINE divergst(klevel,x,y,div)
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
  USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm

  IMPLICIT NONE
  !
  ! P. Le Van
  !
  !  ******************************************************************
  !  ... calcule la divergence a tous les niveaux d'1 vecteur de compos. x et y...
  !       x et y  etant des composantes contravariantes   ...
  !  ****************************************************************
  !  x  et  y  sont des arguments  d'entree pour le s-prog
  !    div      est  un argument  de sortie pour le s-prog
  !
  !
  !   -------------------------------------------------------------------
  !
  INCLUDE "comgeom.h"

  INTEGER :: klevel
  REAL :: x( ip1jmp1,klevel ),y( ip1jm,klevel ),div( ip1jmp1,klevel )
  INTEGER :: ij,l,i
  REAL :: aiy1( iip1 ) , aiy2( iip1 )
  REAL :: sumypn,sumyps

  REAL :: SSUM
  !
  !
  DO l = 1,klevel
  !
  DO ij = iip2, ip1jm - 1
  div( ij + 1, l ) = x(ij+1,l) - x(ij,l)+ y(ij-iim,l)-y(ij+1,l)
  END DO
  !
  ! ....  correction pour  div( 1,j,l)  ......
  ! ....   div(1,j,l)= div(iip1,j,l) ....
  !
  !DIR$ IVDEP
  DO ij = iip2,ip1jm,iip1
  div( ij,l ) = div( ij + iim,l )
  END DO
  !
  ! ....  calcul  aux poles  .....
  !
  !
  DO i  = 1,iim
  aiy1(i)= y(i,l)
  aiy2(i)= y(i+ip1jmi1,l)
  END DO
  sumypn = SSUM ( iim,aiy1,1 )
  sumyps = SSUM ( iim,aiy2,1 )
  DO i = 1,iip1
  div(     i    , l ) = - sumypn/iim
  div( i + ip1jm, l ) =   sumyps/iim
  END DO
  !
  END DO
  RETURN
END SUBROUTINE divergst