Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/comconst.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/comconst.h	(revision 1188)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/comconst.h	(revision 1190)
@@ -8,5 +8,6 @@
      & dtvr,daysec,                                                     &
      & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
-     &                   ,dissip_factz,dissip_deltaz,dissip_zref
+     &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
+     &                   ,iflag_top_bound,tau_top_bound
 
 
@@ -16,4 +17,7 @@
       REAL cotot,unsim,g,omeg
       REAL dissip_factz,dissip_deltaz,dissip_zref
+      INTEGER iflag_top_bound
+      REAL tau_top_bound
+
 
 !-----------------------------------------------------------------------
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/conf_gcm.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/conf_gcm.F	(revision 1188)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/conf_gcm.F	(revision 1190)
@@ -280,4 +280,9 @@
        CALL getin('dissip_deltaz',dissip_deltaz )
        CALL getin('dissip_zref',dissip_zref )
+
+       iflag_top_bound=1
+       tau_top_bound=1.e-5
+       CALL getin('iflag_top_bound',iflag_top_bound)
+       CALL getin('tau_top_bound',tau_top_bound)
 
 !Config  Key  = coefdis
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/leapfrog.F	(revision 1188)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/leapfrog.F	(revision 1190)
@@ -386,5 +386,5 @@
 
          IF (ok_strato) THEN
-           CALL top_bound( vcov,ucov,teta, dufi,dvfi,dtetafi)
+           CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
          ENDIF
        
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/top_bound.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/top_bound.F	(revision 1188)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/top_bound.F	(revision 1190)
@@ -1,3 +1,3 @@
-      SUBROUTINE top_bound( vcov,ucov,teta, du,dv,dh )
+      SUBROUTINE top_bound( vcov,ucov,teta,masse, du,dv,dh )
       IMPLICIT NONE
 c
@@ -5,5 +5,6 @@
 #include "paramet.h"
 #include "comconst.h"
-CC#include "comgeom2.h"
+#include "comvert.h"
+#include "comgeom2.h"
 
 
@@ -27,5 +28,5 @@
 c   -------------
 
-#include "comgeom.h"
+! #include "comgeom.h"
 #include "comdissipn.h"
 
@@ -34,4 +35,5 @@
 
       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)
 
@@ -39,11 +41,14 @@
 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 :: rdamp(llm) = 
-     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) 
+      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
@@ -52,10 +57,38 @@
 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
-            vzon(j,l)=vzon(j,l)+vcov(i,j,l)/float(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
@@ -72,9 +105,22 @@
         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
-            uzon(j,l)=uzon(j,l)+ucov(i,j,l)/float(iim)
-            tzon(j,l)=tzon(j,l)+teta(i,j,l)/float(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
@@ -85,5 +131,6 @@
         do i=1,iip1
           do j=2,jjm
-            du(i,j,l)=du(i,j,l)-rdamp(l)*(ucov(i,j,l)-uzon(j,l))
+            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
