Index: LMDZ5/branches/testing/libf/dyn3dpar/top_bound_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/top_bound_p.F	(revision 1663)
+++ LMDZ5/branches/testing/libf/dyn3dpar/top_bound_p.F	(revision 1795)
@@ -1,3 +1,6 @@
-      SUBROUTINE top_bound_p( vcov,ucov,teta,masse, du,dv,dh )
+!
+! $Id$
+!
+      SUBROUTINE top_bound_p(vcov,ucov,teta,masse,dt)
       USE parallel
       IMPLICIT NONE
@@ -25,16 +28,39 @@
 c
 c=======================================================================
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
+
+! top_bound sponge layer model:
+! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*t)
+! 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 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
 
 c   Local:
@@ -43,9 +69,7 @@
       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.
       INTEGER j,l,jjb,jje
@@ -53,43 +77,57 @@
 
       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
-! couce 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
+
+! 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
+             write(lunout,'(6(1pe12.4,1x))')
+     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
+     &           1./lambda(l),lambda(l)
+           endif
+         enddo
          first=.false.
-         print*,'TOP_BOUND rdamp=',rdamp
 c$OMP END MASTER
 c$OMP BARRIER
-      endif
+      endif ! of if (first)
 
 
       CALL massbar_p(masse,massebx,masseby)
-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)      
-      do l=1,llm
+
+      ! compute zonal average of vcov (or set it to zero)
+      if (mode_top_bound.ge.2) then
+       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
           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
+! 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)
@@ -97,24 +135,22 @@
           vzon(j,l)=vzon(j,l)/zm
         enddo
-      enddo
+       enddo
 c$OMP END DO NOWAIT   
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      do l=1,llm
-        do j=jjb,jje
-          do i=1,iip1
-            dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
-          enddo
-        enddo
-      enddo
-c$OMP END DO NOWAIT
-
-      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
+      else
+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
+       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
           uzon(j,l)=0.
@@ -126,9 +162,22 @@
           uzon(j,l)=uzon(j,l)/zm
         enddo
-      enddo
-c$OMP END DO NOWAIT
-
+       enddo
+c$OMP END DO NOWAIT
+      else
+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
+       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 l=1,llm
         do j=jjb,jje
           zm=0.
@@ -140,22 +189,61 @@
           tzon(j,l)=tzon(j,l)/zm
         enddo
-      enddo
-c$OMP END DO NOWAIT
-
-C   AMORTISSEMENTS LINEAIRES:
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      do l=1,llm
+       enddo
+c$OMP END DO NOWAIT
+      endif ! of if (mode_top_bound.ge.3)
+
+      if (mode_top_bound.ge.1) then
+       ! 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)=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
-c$OMP END DO NOWAIT
-      
-
-      RETURN
+            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
+            ucov(i,j,l)=ucov(i,j,l)
+     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(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
