source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dissip_loc.f90 @ 5153

! $Id: $

SUBROUTINE dissip_loc(vcov, ucov, teta, p, dv, du, dh)
  !
  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 lmdz_comdissipn, ONLY: tetaudiv, tetaurot, tetah, cdivu, crot, cdivh
  USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
          tetagrot, tetatemp, coefdis, vert_prof_dissip
  USE lmdz_comgeom

  IMPLICIT NONE


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

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

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


END SUBROUTINE dissip_loc