source: LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_dissip.f90 @ 5218

MODULE lmdz_dissip
  IMPLICIT NONE; PRIVATE
  PUBLIC dissip

CONTAINS

  SUBROUTINE dissip(vcov, ucov, teta, p, dv, du, dh)
    USE comconst_mod, ONLY: dtdiss
    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

    USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
    USE lmdz_paramet
    IMPLICIT NONE


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

    !=======================================================================

    !   Auteur:  P. Le Van
    !   -------

    !   Objet:
    !   ------

    !   Dissipation horizontale

    !=======================================================================
    !-----------------------------------------------------------------------
    !   Declarations:
    !   -------------




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

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

    !   Local:
    !   ------

    REAL :: gdx(ip1jmp1, llm), gdy(ip1jm, llm)
    REAL :: grx(ip1jmp1, llm), gry(ip1jm, llm)
    REAL :: te1dt(llm), te2dt(llm), te3dt(llm)
    REAL :: deltapres(ip1jmp1, llm)

    INTEGER :: l, ij

    !-----------------------------------------------------------------------
    !   initialisations:
    !   ----------------

    DO l = 1, llm
      te1dt(l) = tetaudiv(l) * dtdiss
      te2dt(l) = tetaurot(l) * dtdiss
      te3dt(l) = tetah(l) * dtdiss
    ENDDO
    du = 0.
    dv = 0.
    dh = 0.

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

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

    IF(lstardis) THEN
      CALL gradiv2(llm, ucov, vcov, nitergdiv, gdx, gdy)
    ELSE
      CALL gradiv (llm, ucov, vcov, nitergdiv, gdx, gdy)
    ENDIF

    DO l = 1, llm

      DO ij = 1, iip1
        gdx(ij, l) = 0.
        gdx(ij + ip1jm, l) = 0.
      ENDDO

      DO ij = iip2, ip1jm
        du(ij, l) = du(ij, l) - te1dt(l) * gdx(ij, l)
      ENDDO
      DO ij = 1, ip1jm
        dv(ij, l) = dv(ij, l) - te1dt(l) * gdy(ij, l)
      ENDDO

    ENDDO

    !   calcul de la partie   n X grad ( rot ):
    !   ---------------------------------------

    IF(lstardis) THEN
      CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry)
    ELSE
      CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry)
    ENDIF

    DO l = 1, llm
      DO ij = 1, iip1
        grx(ij, l) = 0.
      ENDDO

      DO ij = iip2, ip1jm
        du(ij, l) = du(ij, l) - te2dt(l) * grx(ij, l)
      ENDDO
      DO ij = 1, ip1jm
        dv(ij, l) = dv(ij, l) - te2dt(l) * gry(ij, l)
      ENDDO
    ENDDO

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

    IF(lstardis) THEN

      DO l = 1, llm
        DO ij = 1, ip1jmp1
          deltapres(ij, l) = AMAX1(0., p(ij, l) - p(ij, l + 1))
        ENDDO
      ENDDO

      CALL divgrad2(llm, teta, deltapres, niterh, gdx)
    ELSE
      CALL divgrad (llm, teta, niterh, gdx)
    ENDIF

    DO l = 1, llm
      DO ij = 1, ip1jmp1
        dh(ij, l) = dh(ij, l) - te3dt(l) * gdx(ij, l)
      ENDDO
    ENDDO

  END SUBROUTINE dissip


END MODULE lmdz_dissip