!
! $Id: $
!
SUBROUTINE dissip_loc( vcov,ucov,teta,p, dv,du,dh )
  !
  USE comdissipn_mod_h
  USE comdissnew_mod_h
  USE parallel_lmdz
  USE write_field_loc
  USE dissip_mod, ONLY: dissip_allocate
  USE comconst_mod, ONLY: dtdiss
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
  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


  ! ..  Avec nouveaux operateurs star :  gradiv2 , divgrad2, nxgraro2  ...
                              ! (  10/01/98  )

  !=======================================================================
  !
  !   Auteur:  P. Le Van
  !   -------
  !
  !   Objet:
  !   ------
  !
  !   Dissipation horizontale
  !
  !=======================================================================
  !-----------------------------------------------------------------------
  !   Declarations:
  !   -------------



  include "comgeom.h"

  !   Arguments:
  !   ----------

  REAL,INTENT(IN) :: vcov(ijb_v:ije_v,llm) ! covariant meridional wind
  REAL,INTENT(IN) :: ucov(ijb_u:ije_u,llm) ! covariant zonal wind
  REAL,INTENT(IN) :: teta(ijb_u:ije_u,llm) ! potential temperature
  REAL,INTENT(IN) :: p(ijb_u:ije_u,llmp1) ! interlayer pressure
  ! ! tendencies (.../s) on covariant winds and potential temperature
  REAL,INTENT(OUT) :: dv(ijb_v:ije_v,llm)
  REAL,INTENT(OUT) :: du(ijb_u:ije_u,llm)
  REAL,INTENT(OUT) :: dh(ijb_u:ije_u,llm)

  !   Local:
  !   ------

  REAL :: gdx(ijb_u:ije_u,llm),gdy(ijb_v:ije_v,llm)
  REAL :: grx(ijb_u:ije_u,llm),gry(ijb_v:ije_v,llm)
  REAL :: te1dt(llm),te2dt(llm),te3dt(llm)
  REAL :: deltapres(ijb_u:ije_u,llm)

  INTEGER :: l,ij

  REAL :: SSUM
  integer :: ijb,ije

  LOGICAl,SAVE :: first=.TRUE.
!$OMP THREADPRIVATE(first)

  IF (first) THEN
    CALL dissip_allocate
    first=.FALSE.
  ENDIF
  !-----------------------------------------------------------------------
  !   initialisations:
  !   ----------------

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
     te1dt(l) = tetaudiv(l) * dtdiss
     te2dt(l) = tetaurot(l) * dtdiss
     te3dt(l) = tetah(l)    * dtdiss
  ENDDO
!$OMP END DO NOWAIT
   ! CALL initial0( ijp1llm, du )
   ! CALL initial0( ijmllm , dv )
   ! CALL initial0( ijp1llm, dh )

  ijb=ij_begin
  ije=ij_end

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    du(ijb:ije,l)=0
    dh(ijb:ije,l)=0
  ENDDO
!$OMP END DO NOWAIT

  if (pole_sud) ije=ij_end-iip1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    dv(ijb:ije,l)=0
  ENDDO
!$OMP END DO NOWAIT

  !-----------------------------------------------------------------------
  !   Calcul de la dissipation:
  !   -------------------------

  !   Calcul de la partie   grad  ( div ) :
  !   -------------------------------------



  IF(lstardis) THEN
   ! IF (.FALSE.) THEN
     CALL gradiv2_loc( llm,ucov,vcov,nitergdiv,gdx,gdy )
  ELSE
      ! CALL gradiv_p ( llm,ucov,vcov,nitergdiv,gdx,gdy )
  ENDIF

IF (CPPKEY_DEBUGIO) THEN
  call WriteField_u('gdx',gdx)
  call WriteField_v('gdy',gdy)
END IF

  ijb=ij_begin
  ije=ij_end
  if (pole_sud) ije=ij_end-iip1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
     if (pole_nord) then
       DO ij = 1, iip1
          gdx(     ij ,l) = 0.
       ENDDO
     endif

     if (pole_sud) then
       DO ij = 1, iip1
          gdx(ij+ip1jm,l) = 0.
       ENDDO
     endif

     if (pole_nord) ijb=ij_begin+iip1
     DO ij = ijb,ije
        du(ij,l) = du(ij,l) - te1dt(l) *gdx(ij,l)
     ENDDO

     if (pole_nord) ijb=ij_begin
     DO ij = ijb,ije
        dv(ij,l) = dv(ij,l) - te1dt(l) *gdy(ij,l)
     ENDDO

   ENDDO
!$OMP END DO NOWAIT
  !   calcul de la partie   n X grad ( rot ):
  !   ---------------------------------------

  IF(lstardis) THEN
   ! IF (.FALSE.) THEN
     CALL nxgraro2_loc( llm,ucov, vcov, nitergrot,grx,gry )
  ELSE
      ! CALL nxgrarot_p( llm,ucov, vcov, nitergrot,grx,gry )
  ENDIF

IF (CPPKEY_DEBUGIO) THEN
  call WriteField_u('grx',grx)
  call WriteField_v('gry',gry)
END IF


  ijb=ij_begin
  ije=ij_end
  if (pole_sud) ije=ij_end-iip1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm

     if (pole_nord) then
       DO ij = 1, iip1
          grx(ij,l) = 0.
       ENDDO
     endif

     if (pole_nord) ijb=ij_begin+iip1
     DO ij = ijb,ije
        du(ij,l) = du(ij,l) - te2dt(l) * grx(ij,l)
     ENDDO

     if (pole_nord) ijb=ij_begin
     DO ij =  ijb, ije
        dv(ij,l) = dv(ij,l) - te2dt(l) * gry(ij,l)
     ENDDO

  ENDDO
!$OMP END DO NOWAIT

  !   calcul de la partie   div ( grad ):
  !   -----------------------------------


  IF(lstardis) THEN
   ! IF (.FALSE.) THEN

  ijb=ij_begin
  ije=ij_end

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
   DO l = 1, llm
      DO ij = ijb, ije
        deltapres(ij,l) = AMAX1( 0.,  p(ij,l) - p(ij,l+1) )
      ENDDO
   ENDDO
!$OMP END DO NOWAIT
     CALL divgrad2_loc( llm,teta, deltapres  ,niterh, gdx )
  ELSE
      ! CALL divgrad_p ( llm,teta, niterh, gdx        )
  ENDIF

IF (CPPKEY_DEBUGIO) THEN
  call WriteField_u('gdx',gdx)
END IF


  ijb=ij_begin
  ije=ij_end

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1,llm
     DO ij = ijb,ije
        dh( ij,l ) = dh( ij,l ) - te3dt(l) * gdx( ij,l )
     ENDDO
  ENDDO
!$OMP END DO NOWAIT

  RETURN
END SUBROUTINE dissip_loc