! $Header$

SUBROUTINE divergst(klevel, x, y, div)
  USE lmdz_ssum_scopy, ONLY: ssum
  USE lmdz_comgeom

  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 "dimensions.h"
  INCLUDE "paramet.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
  !
  !
  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