Index: trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1000)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1010)
@@ -370,15 +370,21 @@
        CALL getin('dissip_pupstart',dissip_pupstart )
 
-! Parametres controlant la couche eponge
-! Actifs uniquement avec ok_strato=y
-!   mode = 0 : pas de sponge
-!   mode = 1 : u et v -> 0
-!   mode = 2 : u et v -> moyenne zonale
-!   mode = 3 : u, v et h -> moyenne zonale
+! top_bound sponge: only active if ok_strato=.true. and iflag_top_bound!=0
+!                   iflag_top_bound=0 for no sponge
+!                   iflag_top_bound=1 for sponge over 4 topmost layers
+!                   iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
        iflag_top_bound=1
+       CALL getin('iflag_top_bound',iflag_top_bound)
+
+! mode_top_bound : fields towards which sponge relaxation will be done:
+!                  mode_top_bound=0: no relaxation
+!                  mode_top_bound=1: u and v relax towards 0
+!                  mode_top_bound=2: u and v relax towards their zonal mean
+!                  mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
        mode_top_bound=3
+       CALL getin('mode_top_bound',mode_top_bound)
+
+! top_bound sponge : inverse of charactericstic relaxation time scale for sponge
        tau_top_bound=1.e-5
-       CALL getin('iflag_top_bound',iflag_top_bound)
-       CALL getin('mode_top_bound',mode_top_bound)
        CALL getin('tau_top_bound',tau_top_bound)
 
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1000)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1010)
@@ -113,8 +113,9 @@
 
 c   tendances top_bound (sponge layer)
-      REAL,SAVE :: dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
-      REAL,SAVE :: dtetatop(ip1jmp1,llm)
-      REAL,SAVE :: dptop(ip1jmp1)
-      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqtop
+c      REAL,SAVE :: dvtop(ip1jm,llm)
+      REAL,SAVE :: dutop(ip1jmp1,llm)
+c      REAL,SAVE :: dtetatop(ip1jmp1,llm)
+c      REAL,SAVE :: dptop(ip1jmp1)
+c      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqtop
 
 c   TITAN : tendances due au forces de marees */s
@@ -247,5 +248,5 @@
             ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
             ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
-            ALLOCATE(dqtop(ip1jmp1,llm,nqtot))
+c            ALLOCATE(dqtop(ip1jmp1,llm,nqtot))
          END IF
 c$OMP END MASTER      
@@ -262,8 +263,8 @@
         dtetadis(:,:)=0.
         dutop(:,:)   =0.
-        dvtop(:,:)   =0.
-        dtetatop(:,:)=0.
-        dqtop(:,:,:) =0.
-        dptop(:)     =0.
+c        dvtop(:,:)   =0.
+c        dtetatop(:,:)=0.
+c        dqtop(:,:,:) =0.
+c        dptop(:)     =0.
         dufi(:,:)   =0.
         dvfi(:,:)   =0.
@@ -1001,11 +1002,14 @@
 c      -------------------
          IF (ok_strato) THEN
-           CALL top_bound_p( vcov,ucov,teta,phi,masse,
-     $                       dutop,dvtop,dtetatop)
-           CALL addfi_p( dtphys, leapf, forward   ,
-     $                  ucov, vcov, teta , q   ,ps ,
-     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
-
-         ENDIF
+           CALL top_bound_p(vcov,ucov,teta,masse,dtphys,
+     $                       dutop)
+        ijb=ij_begin
+        ije=ij_end
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)        
+        DO l=1,llm
+          dutop(ijb:ije,l)=dutop(ijb:ije,l)/dtphys ! convert to a tendency in (m/s)/s
+        ENDDO
+c$OMP END DO NOWAIT        
+         ENDIF ! of IF (ok_strato)
        
 c$OMP BARRIER
@@ -1124,31 +1128,13 @@
         ! Sponge layer (if any)
         IF (ok_strato) THEN
-          ! set dutop,dvtop,... to zero
+          CALL top_bound_p(vcov,ucov,teta,masse,dtvr,
+     $                     dutop)
           ijb=ij_begin
           ije=ij_end
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-          do l=1,llm
-            dutop(ijb:ije,l)=0
-            dtetatop(ijb:ije,l)=0
-            dqtop(ijb:ije,l,1:nqtot)=0
-          enddo
-!$OMP END DO
-!$OMP SINGLE
-          dptop(ijb:ije)=0
-!$OMP END SINGLE
-          ijb=ij_begin
-          ije=ij_end
-          if (pole_sud) ije=ije-iip1
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-          do l=1,llm
-            dvtop(ijb:ije,l)=0
-          enddo
-!$OMP END DO
-
-          CALL top_bound_p(vcov,ucov,teta,phi,masse,
-     $                     dutop,dvtop,dtetatop)
-          CALL addfi_p( dtvr, leapf, forward   ,
-     $                  ucov, vcov, teta , q   ,ps ,
-     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)        
+          DO l=1,llm
+            dutop(ijb:ije,l)=dutop(ijb:ije,l)/dtvr ! convert to a tendency in (m/s)/s
+          ENDDO
+c$OMP END DO NOWAIT        
 !$OMP BARRIER
         ENDIF ! of IF (ok_strato) 
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/top_bound_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/top_bound_p.F	(revision 1000)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/top_bound_p.F	(revision 1010)
@@ -1,3 +1,6 @@
-      SUBROUTINE top_bound_p( vcov,ucov,teta,phi,masse, du,dv,dh )
+!
+! $Id: top_bound_p.F 1793 2013-07-18 07:13:18Z emillour $
+!
+      SUBROUTINE top_bound_p(vcov,ucov,teta,masse,dt,ducov)
       USE parallel
       IMPLICIT NONE
@@ -25,17 +28,40 @@
 c
 c=======================================================================
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
+
+! top_bound sponge layer model:
+! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*dt)
+! where Am is the zonal average of the field (or zero), and lambda the inverse
+! of the characteristic quenching/relaxation time scale
+! Thus, assuming Am to be time-independent, field at time t+dt is given by:
+! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*t))
+! Moreover lambda can be a function of model level (see below), and relaxation
+! can be toward the average zonal field or just zero (see below).
+
+! NB: top_bound sponge is only called from leapfrog if ok_strato=.true.
+
+! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst.h)
+!    iflag_top_bound=0 for no sponge
+!    iflag_top_bound=1 for sponge over 4 topmost layers
+!    iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
+!    mode_top_bound=0: no relaxation
+!    mode_top_bound=1: u and v relax towards 0
+!    mode_top_bound=2: u and v relax towards their zonal mean
+!    mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
+!    tau_top_bound : inverse of charactericstic relaxation time scale at
+!                       the topmost layer (Hz)
+
 
 #include "comdissipn.h"
+#include "iniprint.h"
 
 c   Arguments:
 c   ----------
 
-      REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm)
-      REAL phi(iip1,jjp1,llm)                  ! geopotentiel
-      REAL masse(iip1,jjp1,llm)
-      REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm)
+      real,intent(inout) :: ucov(iip1,jjp1,llm) ! covariant zonal wind
+      real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
+      real,intent(inout) :: teta(iip1,jjp1,llm) ! potential temperature
+      real,intent(in) :: masse(iip1,jjp1,llm) ! mass of atmosphere 
+      real,intent(in) :: dt ! time step (s) of sponge model
+      real,intent(out) :: ducov(iip1,jjp1,llm) ! increment on ucov due to sponge
 
 c   Local:
@@ -44,64 +70,58 @@
       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/) 
+      REAL,SAVE :: rdamp(llm) ! quenching coefficient
+      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
       LOGICAL,SAVE :: first=.true.
-      
-      REAL zkm
       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.
+! sponge quenching over the topmost 4 atmospheric layers
+             lambda(:)=0.
+             lambda(llm)=tau_top_bound
+             lambda(llm-1)=tau_top_bound/2.
+             lambda(llm-2)=tau_top_bound/4.
+             lambda(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
+! sponge quenching over topmost layers down to pressures which are
+! higher than 100 times the topmost layer pressure
+             lambda(:)=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)'
+
+! quenching coefficient rdamp(:)
+!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
+         rdamp(:)=1.-exp(-lambda(:)*dt)
+
+         write(lunout,*)'TOP_BOUND mode',mode_top_bound
+         write(lunout,*)'Sponge layer coefficients'
+         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
          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)
+             write(lunout,'(6(1pe12.4,1x))')
+     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
+     &           1./lambda(l),lambda(l)
            endif
          enddo
+         first=.false.
 c$OMP END MASTER
 c$OMP BARRIER
-      endif
+      endif ! of if (first)
 
 
       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
-
+      ! compute zonal average of vcov (or set it to zero)
       if (mode_top_bound.ge.2) then
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+       jjb=jj_begin
+       jje=jj_end
+       IF (pole_sud) jje=jj_end-1
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
        do l=1,llm
         do j=jjb,jje
@@ -109,7 +129,6 @@
           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
+! NB: we can work using vcov zonal mean rather than v since the
+! cv coefficient (which relates the two) only varies with latitudes 
             vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
             zm=zm+masseby(i,j,l)
@@ -118,40 +137,20 @@
         enddo
        enddo
-!$OMP END DO NOWAIT   
+c$OMP END DO NOWAIT   
       else
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-       do l=1,llm
-        do j=jjb,jje
-          vzon(j,l)=0.
-        enddo
-       enddo
-!$OMP END DO NOWAIT
-      endif
-
-C   AMORTISSEMENTS LINEAIRES:
-
-      if (mode_top_bound.ge.1) then
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-       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
-!$OMP END DO NOWAIT
-      endif
-
-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)      
+       do l=1,llm
+         vzon(:,l)=0.
+       enddo
+c$OMP END DO NOWAIT
+      endif ! of if (mode_top_bound.ge.2)
+
+      ! compute zonal average of u (or set it to zero)
       if (mode_top_bound.ge.2) then
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+       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)      
        do l=1,llm
         do j=jjb,jje
@@ -165,17 +164,20 @@
         enddo
        enddo
-!$OMP END DO NOWAIT
+c$OMP END DO NOWAIT
       else
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-       do l=1,llm
-        do j=jjb,jje
-          uzon(j,l)=0.
-        enddo
-       enddo
-!$OMP END DO NOWAIT
-      endif
-
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+       do l=1,llm
+         uzon(:,l)=0.
+       enddo
+c$OMP END DO NOWAIT
+      endif ! of if (mode_top_bound.ge.2)
+
+      ! compute zonal average of potential temperature, if necessary
       if (mode_top_bound.ge.3) then
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+       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)   
        do l=1,llm
         do j=jjb,jje
@@ -189,34 +191,61 @@
         enddo
        enddo
-!$OMP END DO NOWAIT
-      endif
-
-C   AMORTISSEMENTS LINEAIRES:
+c$OMP END DO NOWAIT
+      endif ! of if (mode_top_bound.ge.3)
 
       if (mode_top_bound.ge.1) then
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
+       ! Apply sponge quenching on vcov:
+       jjb=jj_begin
+       jje=jj_end
+       IF (pole_sud) jje=jj_end-1
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
        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
-!$OMP END DO NOWAIT
-      endif
-      
-      if (mode_top_bound.ge.3) then
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+            vcov(i,j,l)=vcov(i,j,l)
+     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
+          enddo
+        enddo
+       enddo
+c$OMP END DO NOWAIT
+
+       ! Apply sponge quenching on ucov:
+       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)      
        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
-!$OMP END DO NOWAIT
-      endif      
-
-!$OMP BARRIER
-      RETURN
+            ducov(i,j,l)=-rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
+            ucov(i,j,l)=ucov(i,j,l)
+     &                  +ducov(i,j,l)
+          enddo
+       enddo
+       enddo
+c$OMP END DO NOWAIT
+      endif ! of if (mode_top_bound.ge.1)
+
+      if (mode_top_bound.ge.3) then   
+       ! Apply sponge quenching on teta:
+       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)      
+       do l=1,llm
+        do j=jjb,jje
+          do i=1,iip1
+            teta(i,j,l)=teta(i,j,l)
+     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
+          enddo
+       enddo
+       enddo
+c$OMP END DO NOWAIT
+      endif ! of if (mode_top_bond.ge.3)
+
       END
