source: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/top_bound.F @ 1233

      SUBROUTINE top_bound( vcov,ucov,teta,masse, du,dv,dh )
      IMPLICIT NONE
c
#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comvert.h"
#include "comgeom2.h"


c ..  DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,
C     F. LOTT DEC. 2006
c                                 (  10/12/06  )

c=======================================================================
c
c   Auteur:  F. LOTT  
c   -------
c
c   Objet:
c   ------
c
c   Dissipation linéaire (ex top_bound de la physique)
c
c=======================================================================
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

! #include "comgeom.h"
#include "comdissipn.h"

c   Arguments:
c   ----------

      REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm)
      REAL masse(iip1,jjp1,llm)
      REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm)

c   Local:
c   ------

      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
      
      INTEGER NDAMP
      PARAMETER (NDAMP=4)
      integer i
      REAL,SAVE :: rdamp(llm)
!     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) 

      LOGICAL,SAVE :: first=.true.

      INTEGER j,l


C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
      
      if (iflag_top_bound.eq.0) return

      if (first) then
         if (iflag_top_bound.eq.1) then
! couche eponge dans les 4 dernieres couches du modele
             rdamp(:)=0.
             rdamp(llm)=tau_top_bound
             rdamp(llm-1)=tau_top_bound/2.
             rdamp(llm-2)=tau_top_bound/4.
             rdamp(llm-3)=tau_top_bound/8.
         else if (iflag_top_bound.eq.2) then
! couce eponge dans toutes les couches de pression plus faible que
! 100 fois la pression de la derniere couche
             rdamp(:)=tau_top_bound
     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
         endif
         first=.false.
         print*,'TOP_BOUND rdamp=',rdamp
      endif

      CALL massbar(masse,massebx,masseby)

      do l=1,llm
        do j=1,jjm
          vzon(j,l)=0.
          zm=0.
          do i=1,iim
! Rm: on peut travailler directement avec la moyenne zonale de vcov
! plutot qu'avec celle de v car le coefficient cv qui relie les deux
! ne varie qu'en latitude
            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
            zm=zm+masseby(i,j,l)
          enddo
          vzon(j,l)=vzon(j,l)/zm
        enddo
      enddo

      do l=1,llm
        do i=1,iip1
          do j=1,jjm
            dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
          enddo
        enddo
      enddo

      do l=1,llm
        do j=2,jjm
          uzon(j,l)=0.
          zm=0.
          do i=1,iim
            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
            zm=zm+massebx(i,j,l)
          enddo
          uzon(j,l)=uzon(j,l)/zm
        enddo
      enddo

      do l=1,llm
        do j=2,jjm
          zm=0.
          tzon(j,l)=0.
          do i=1,iim
            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
            zm=zm+masse(i,j,l)
          enddo
          tzon(j,l)=tzon(j,l)/zm
        enddo
      enddo

C   AMORTISSEMENTS LINEAIRES:

      do l=1,llm
        do i=1,iip1
          do j=2,jjm
            du(i,j,l)=du(i,j,l)
     s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
            dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
          enddo
        enddo
      enddo
      

      RETURN
      END