SUBROUTINE top_bound_p( vcov,ucov,teta,masse, du,dv,dh )
      USE parallel
      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 "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,jjb,jje


      if (iflag_top_bound == 0) return
      if (first) then
c$OMP BARRIER
c$OMP MASTER
         if (iflag_top_bound == 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 == 2) then
! couche 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 mode',mode_top_bound
         print*,'Coeffs pour la couche eponge a l equateur'
         print*,'p (Pa)  z(km)  tau (s)   dt*rdamp'
         do l=1,llm
           if (rdamp(l).ne.0.) then
            zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
          print*,presnivs(l),zkm,
     .          1./rdamp(l),
     .          dt*rdamp(l)
           endif
         enddo
c$OMP END MASTER
c$OMP BARRIER
      endif


      CALL massbar_p(masse,massebx,masseby)

c   mode = 0 : pas de sponge
c   mode = 1 : u et v -> 0
c   mode = 2 : u et v -> moyenne zonale
c   mode = 3 : u, v et h -> moyenne zonale

C POUR V

C  CALCUL DES CHAMPS EN MOYENNE ZONALE:

      jjb=jj_begin
      jje=jj_end
      IF (pole_sud) jje=jj_end-1

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
      if (mode_top_bound.ge.2) then
       do l=1,llm
        do j=jjb,jje
          zm=0.
          vzon(j,l)=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
      else
       do l=1,llm
        do j=jjb,jje
          vzon(j,l)=0.
        enddo
       enddo
      endif
c$OMP END DO NOWAIT   

C   AMORTISSEMENTS LINEAIRES:

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
      if (mode_top_bound.ge.1) then
       do l=1,llm
        do j=jjb,jje
          do i=1,iip1
            dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
          enddo
        enddo
       enddo
      endif
c$OMP END DO NOWAIT

C POUR U ET H

C  CALCUL DES CHAMPS EN MOYENNE ZONALE:

      jjb=jj_begin
      jje=jj_end
      IF (pole_nord) jjb=jj_begin+1
      IF (pole_sud)  jje=jj_end-1

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
      if (mode_top_bound.ge.2) then
       do l=1,llm
        do j=jjb,jje
          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
      else
       do l=1,llm
        do j=jjb,jje
          uzon(j,l)=0.
        enddo
       enddo
      endif
c$OMP END DO NOWAIT

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
      if (mode_top_bound.ge.3) then
       do l=1,llm
        do j=jjb,jje
          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
      endif
c$OMP END DO NOWAIT

C   AMORTISSEMENTS LINEAIRES:

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
      if (mode_top_bound.ge.1) then
       do l=1,llm
        do j=jjb,jje
          do i=1,iip1
            du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
          enddo
        enddo
       enddo
      endif
c$OMP END DO NOWAIT
      
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
      if (mode_top_bound.ge.3) then
       do l=1,llm
        do j=jjb,jje
          do i=1,iip1
            dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
          enddo
        enddo
       enddo
      endif
c$OMP END DO NOWAIT
      

      RETURN
      END