! $Header$

SUBROUTINE divgrad2(klevel, h, deltapres, lh, divgra)
  !
  ! P. Le Van
  !
  !   ***************************************************************
  !
  ! .....   calcul de  (div( grad ))   de (  pext * h ) .....
  !   ****************************************************************
  !   h ,klevel,lh et pext  sont des arguments  d'entree pour le s-prg
  !     divgra     est  un argument  de sortie pour le s-prg
  !
  USE lmdz_ssum_scopy, ONLY: scopy
  USE lmdz_comdissipn, ONLY: tetaudiv, tetaurot, tetah, cdivu, crot, cdivh
  USE lmdz_comgeom2

  IMPLICIT NONE
  !
  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"

  !    .......    variables en arguments   .......
  !
  INTEGER :: klevel
  REAL :: h(ip1jmp1, klevel), deltapres(ip1jmp1, klevel)
  REAL :: divgra(ip1jmp1, klevel)
  !
  !    .......    variables  locales    ..........
  !
  REAL :: signe, nudivgrs, sqrtps(ip1jmp1, llm)
  INTEGER :: l, ij, iter, lh
  !    ...................................................................

  !
  signe = (-1.)**lh
  nudivgrs = signe * cdivh

  CALL SCOPY (ip1jmp1 * klevel, h, 1, divgra, 1)

  !
  CALL laplacien(klevel, divgra, divgra)

  DO l = 1, klevel
    DO ij = 1, ip1jmp1
      sqrtps(ij, l) = SQRT(deltapres(ij, l))
    ENDDO
  ENDDO
  !
  DO l = 1, klevel
    DO ij = 1, ip1jmp1
      divgra(ij, l) = divgra(ij, l) * sqrtps(ij, l)
    ENDDO
  ENDDO

  !    ........    Iteration de l'operateur  laplacien_gam    ........
  !
  DO  iter = 1, lh - 2
    CALL laplacien_gam (klevel, cuvscvgam2, cvuscugam2, unsair_gam2, &
            unsapolnga2, unsapolsga2, divgra, divgra)
  ENDDO
  !
  !    ...............................................................

  DO l = 1, klevel
    DO ij = 1, ip1jmp1
      divgra(ij, l) = divgra(ij, l) * sqrtps(ij, l)
    ENDDO
  ENDDO
  !
  CALL laplacien (klevel, divgra, divgra)
  !
  DO l = 1, klevel
    DO ij = 1, ip1jmp1
      divgra(ij, l) = nudivgrs * divgra(ij, l) / deltapres(ij, l)
    ENDDO
  ENDDO

  RETURN
END SUBROUTINE divgrad2