Index: trunk/LMDZ.COMMON/libf/dyn3d/bilan_dyn.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/bilan_dyn.F	(revision 1012)
+++ trunk/LMDZ.COMMON/libf/dyn3d/bilan_dyn.F	(revision 1017)
@@ -20,4 +20,5 @@
 
       USE control_mod, ONLY: planet_type 
+      USE cpdet_mod, only: tpot2t
 
       IMPLICIT NONE
Index: trunk/LMDZ.COMMON/libf/dyn3d/calfis.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/calfis.F	(revision 1012)
+++ trunk/LMDZ.COMMON/libf/dyn3d/calfis.F	(revision 1017)
@@ -31,5 +31,5 @@
       USE infotrac
       USE control_mod
- 
+      USE cpdet_mod, only: t2tpot,tpot2t
 
       IMPLICIT NONE
Index: trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 1012)
+++ trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 1017)
@@ -15,4 +15,5 @@
       USE infotrac, ONLY : type_trac
       use assert_m, only: assert
+      use sponge_mod, only: callsponge,mode_sponge,nsponge,tetasponge
 
       IMPLICIT NONE
@@ -361,4 +362,33 @@
        tau_top_bound=1.e-5
        CALL getin('tau_top_bound',tau_top_bound)
+
+! the other possible sponge layer (sponge_mod)
+       callsponge=.false. ! default value; don't use the sponge
+       call getin("callsponge",callsponge)
+       ! check that user is not trying to use both sponge models
+       if ((iflag_top_bound.ge.1).and.callsponge) then
+         write(lunout,*)'Bad choice of options:'
+         write(lunout,*)' iflag_top_bound=',iflag_top_bound
+         write(lunout,*)' and callsponge=.true.'
+         write(lunout,*)'But both sponge models should not be',
+     &                  ' used simultaneously!'
+         stop
+       endif
+       
+! nsponge: number of atmospheric layers over which the sponge extends
+       nsponge=3 ! default value
+       call getin("nsponge",nsponge)
+
+! mode_sponge: (quenching is towards ... over the upper nsponge layers)
+!      0: (h=hmean,u=v=0)
+!      1: (h=hmean,u=umean,v=0)
+!      2: (h=hmean,u=umean,v=vmean)"
+       mode_sponge=2 ! default value
+       call getin("mode_sponge",mode_sponge)
+
+! tetasponge: characteristic time scale (seconds) at topmost layer
+!            (time scale then doubles with decreasing layer index)."
+       tetasponge=50000.0
+       call getin("tetasponge",tetasponge)
 
 ! FOR TITAN: tidal forces
Index: trunk/LMDZ.COMMON/libf/dyn3d/cpdet.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/cpdet.F	(revision 1012)
+++ 	(revision )
@@ -1,134 +1,0 @@
-! ADAPTATION GCM POUR CP(T)
-c======================================================================
-c S. Lebonnois, 10/2010
-c
-c Cp doit être calculé par cpdet(t) pour être valable partout
-c
-c La fonction d'Exner reste pk = RCPD*(play/pref)**RKAPPA
-c (RCPD=cpp, RKAPPA=kappa)
-c
-c On passe de T a teta (temperature potentielle) par t2tpot(t,teta,pk)
-c On passe de teta a T par tpot2t(teta,t,pk)
-c
-c======================================================================
-
-      SUBROUTINE ini_cpdet
-      
-      USE control_mod, ONLY: planet_type 
-      IMPLICIT none
-c======================================================================
-c Initialisation de nu_venus et t0_venus
-c======================================================================
-
-! for cpp, nu_venus and t0_venus:
-#include "comconst.h"
-
-      if (planet_type.eq."venus") then
-          nu_venus=0.35
-          t0_venus=460.
-      else
-          nu_venus=0.
-          t0_venus=0.
-      endif
-
-      return
-      end
-
-c======================================================================
-c======================================================================
-
-      FUNCTION cpdet(t)
-
-      USE control_mod, ONLY: planet_type 
-      IMPLICIT none
-
-! for cpp, nu_venus and t0_venus:
-#include "comconst.h"
-
-      real cpdet,t
-
-      if (planet_type.eq."venus") then
-          cpdet = cpp*(t/t0_venus)**nu_venus
-      else
-          cpdet = cpp
-      endif
-
-      return
-      end
-      
-c======================================================================
-c======================================================================
-
-      SUBROUTINE t2tpot(npoints, yt, yteta, ypk)
-c======================================================================
-c Arguments:
-c
-c yt   --------input-R- Temperature
-c yteta-------output-R- Temperature potentielle
-c ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
-c
-c======================================================================
-
-      USE control_mod, ONLY: planet_type 
-      IMPLICIT NONE
-      
-! for cpp, nu_venus and t0_venus:
-#include "comconst.h"
-
-      integer npoints
-      REAL    yt(npoints), yteta(npoints), ypk(npoints)
-      
-      if (planet_type.eq."venus") then
-          yteta = yt**nu_venus                                          &
-     &            - nu_venus * t0_venus**nu_venus * log(ypk/cpp)
-          yteta = yteta**(1./nu_venus)
-      else
-          yteta = yt * cpp/ypk
-      endif
-
-      return
-      end
-
-c======================================================================
-c======================================================================
-
-      SUBROUTINE tpot2t(npoints,yteta, yt, ypk)
-c======================================================================
-c Arguments:
-c
-c yteta--------input-R- Temperature potentielle
-c yt   -------output-R- Temperature
-c ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
-c
-c======================================================================
-
-      USE control_mod, ONLY: planet_type 
-      IMPLICIT NONE
-
-! for cpp, nu_venus and t0_venus:
-#include "comconst.h"
-
-      integer npoints
-      REAL    yt(npoints), yteta(npoints), ypk(npoints)
-      
-      if (planet_type.eq."venus") then
-         yt = yteta**nu_venus                                           &
-     &       + nu_venus * t0_venus**nu_venus * log(ypk/cpp)
-         yt = yt**(1./nu_venus)
-      else
-          yt = yteta * ypk/cpp
-      endif
-  
-      return
-      end
-
-c======================================================================
-c======================================================================
-c
-c ATTENTION
-c
-c Si un jour on a besoin, il faudra coder les routines 
-c    dt2dtpot / dtpto2dt 
-c
-c======================================================================
-c======================================================================
Index: trunk/LMDZ.COMMON/libf/dyn3d/cpdet_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/cpdet_mod.F90	(revision 1017)
+++ trunk/LMDZ.COMMON/libf/dyn3d/cpdet_mod.F90	(revision 1017)
@@ -0,0 +1,144 @@
+      module cpdet_mod
+
+      implicit none
+
+! ADAPTATION OF GCM TO CP(T)
+!======================================================================
+! S. Lebonnois, 10/2010
+!
+! Cp must be computed using cpdet(t) to be valid
+!
+! The Exner function is still pk = RCPD*(play/pref)**RKAPPA
+! (RCPD=cpp, RKAPPA=kappa)
+!
+! One goes from T to teta (potential temperature) using t2tpot(t,teta,pk)
+! One goes from teta to T using tpot2t(teta,t,pk)
+!
+!======================================================================
+
+      contains
+
+      SUBROUTINE ini_cpdet
+      
+      USE control_mod, ONLY: planet_type 
+      IMPLICIT none
+!======================================================================
+! Initialization of nu_venus and t0_venus
+!======================================================================
+
+! for cpp, nu_venus and t0_venus:
+#include "comconst.h"
+
+      if (planet_type.eq."venus") then
+          nu_venus=0.35
+          t0_venus=460.
+      else
+          nu_venus=0.
+          t0_venus=0.
+      endif
+
+      return
+      end subroutine ini_cpdet
+
+!======================================================================
+!======================================================================
+
+      FUNCTION cpdet(t)
+
+      USE control_mod, ONLY: planet_type 
+      IMPLICIT none
+
+! for cpp, nu_venus and t0_venus:
+#include "comconst.h"
+
+      real,intent(in) :: t
+      real cpdet
+
+      if (planet_type.eq."venus") then
+          cpdet = cpp*(t/t0_venus)**nu_venus
+      else
+          cpdet = cpp
+      endif
+
+      return
+      end function cpdet
+      
+!======================================================================
+!======================================================================
+
+      SUBROUTINE t2tpot(npoints, yt, yteta, ypk)
+!======================================================================
+! Arguments:
+!
+! yt   --------input-R- Temperature
+! yteta-------output-R- Temperature potentielle
+! ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
+!
+!======================================================================
+
+      USE control_mod, ONLY: planet_type 
+      IMPLICIT NONE
+      
+! for cpp, nu_venus and t0_venus:
+#include "comconst.h"
+
+      integer,intent(in) :: npoints
+      REAL,intent(in) :: yt(npoints), ypk(npoints)
+      REAL,intent(out) :: yteta(npoints)
+      
+      if (planet_type.eq."venus") then
+          yteta = yt**nu_venus                                          &
+     &            - nu_venus * t0_venus**nu_venus * log(ypk/cpp)
+          yteta = yteta**(1./nu_venus)
+      else
+          yteta = yt * cpp/ypk
+      endif
+
+      return
+      end subroutine t2tpot
+
+!======================================================================
+!======================================================================
+
+      SUBROUTINE tpot2t(npoints,yteta, yt, ypk)
+!======================================================================
+! Arguments:
+!
+! yteta--------input-R- Temperature potentielle
+! yt   -------output-R- Temperature
+! ypk  --------input-R- Fonction d'Exner: RCPD*(pplay/pref)**RKAPPA
+!
+!======================================================================
+
+      USE control_mod, ONLY: planet_type 
+      IMPLICIT NONE
+
+! for cpp, nu_venus and t0_venus:
+#include "comconst.h"
+
+      integer,intent(in) :: npoints
+      REAL,intent(in) :: yteta(npoints), ypk(npoints)
+      REAL,intent(out) :: yt(npoints)
+      
+      if (planet_type.eq."venus") then
+         yt = yteta**nu_venus                                           &
+     &       + nu_venus * t0_venus**nu_venus * log(ypk/cpp)
+         yt = yt**(1./nu_venus)
+      else
+          yt = yteta * ypk/cpp
+      endif
+  
+      return
+      end subroutine tpot2t
+
+!======================================================================
+!======================================================================
+!
+! ATTENTION
+!
+! Si un jour on a besoin, il faudra coder les routines 
+!    dt2dtpot / dtpto2dt 
+!
+!======================================================================
+!======================================================================
+      end module cpdet_mod
Index: trunk/LMDZ.COMMON/libf/dyn3d/gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 1012)
+++ trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 1017)
@@ -16,4 +16,5 @@
       USE infotrac
       USE control_mod
+      use cpdet_mod, only: ini_cpdet
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Index: trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 1012)
+++ trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 1017)
@@ -16,4 +16,6 @@
       USE write_field
       USE control_mod
+      use cpdet_mod, only: cpdet,tpot2t,t2tpot
+      use sponge_mod, only: callsponge,mode_sponge,sponge
       IMPLICIT NONE
 
@@ -197,10 +199,11 @@
       ! for CP(T)
       real :: dtec
-      real,external :: cpdet
       real :: ztetaec(ip1jmp1,llm)
 
 c dummy: sinon cette routine n'est jamais compilee...
       if(1.eq.0) then
+#ifdef CPP_PHYS
         CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
+#endif
       endif
 
@@ -562,4 +565,8 @@
       IF(apdiss) THEN
 
+        ! sponge layer
+        if (callsponge) then
+          CALL sponge(ucov,vcov,teta,ps,dtdiss,mode_sponge)
+        endif
 
 c   calcul de l'energie cinetique avant dissipation
Index: trunk/LMDZ.COMMON/libf/dyn3d/sponge_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/sponge_mod.F90	(revision 1017)
+++ trunk/LMDZ.COMMON/libf/dyn3d/sponge_mod.F90	(revision 1017)
@@ -0,0 +1,162 @@
+module sponge_mod
+
+implicit none
+
+! sponge parameters (set/read via conf_gcm.F)
+logical,save :: callsponge  ! do we use a sponge on upper layers
+integer,save :: mode_sponge ! sponge mode
+integer,save :: nsponge ! number of sponge layers 
+real,save :: tetasponge  ! sponge time scale (s) at topmost layer
+
+
+contains
+
+      subroutine sponge(ucov,vcov,h,ps,dt,mode)
+
+! Sponge routine: Quench ucov, vcov and potential temperature near the
+!                 top of the model
+! Depending on 'mode' relaxation of variables is towards:
+! mode = 0 : h -> h_mean , ucov -> 0 , vcov -> 0
+! mode = 1 : h -> h_mean , ucov -> ucov_mean , vcov -> 0
+! mode >= 2 : h -> h_mean , ucov -> ucov_mean , vcov -> vcov_mean
+! Number of layer over which sponge is applied is 'nsponge' (read from def file)
+! Time scale for quenching at top level is given by 'tetasponge' (read from
+! def file) and doubles as level indexes decrease.
+! Quenching is modeled as: A(t)=Am+A0exp(-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*dt))
+
+
+      implicit none
+#include "dimensions.h"
+#include "paramet.h"
+#include "comdissip.h"
+#include "comvert.h"
+#include "comgeom2.h"
+#include "iniprint.h"
+
+! Arguments:
+!------------
+      real,intent(inout) :: ucov(iip1,jjp1,llm) ! covariant zonal wind
+      real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
+      real,intent(inout) :: h(iip1,jjp1,llm) ! potential temperature
+!      real,intent(in) :: pext(iip1,jjp1) ! extensive pressure
+      real,intent(in) :: ps(iip1,jjp1) ! surface pressure (Pa)
+      real,intent(in) :: dt   ! time step
+      integer,intent(in) :: mode  ! sponge mode
+
+!   Local:
+!   ------
+
+      real,save :: sig_s(llm) !sigma au milieu des couches
+      REAL vm,um,hm,ptot(jjp1)
+      real,save :: cst(llm)
+      real :: pext(iip1,jjp1) ! extensive pressure
+
+      INTEGER l,i,j
+      integer,save :: l0 ! layer down to which sponge is applied
+
+      real ssum
+
+      real zkm
+      logical,save :: firstcall=.true.
+
+
+
+      if (firstcall) then
+
+       ! build approximative sigma levels at midlayer
+        do l=1,llm
+          sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
+        enddo
+
+        l0=llm-nsponge+1
+ 
+        write(lunout,*)
+        write(lunout,*)'sponge mode',mode
+        write(lunout,*)'nsponge tetasponge ',nsponge,tetasponge
+        write(lunout,*)'Coeffs for the sponge layer'
+        write(lunout,*)'Z (km)         tau             cst'
+        do l=llm,l0,-1
+          ! double time scale with every level, starting from the top
+          cst(l)=1.-exp(-dt/(tetasponge*2**(llm-l)))
+        enddo
+
+        do l=l0,llm
+           zkm=-scaleheight*log(sig_s(l))
+           print*,zkm,tetasponge*2**(llm-l),cst(l)
+        enddo
+        PRINT*
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+
+!-----------------------------------------------------------------------
+!   compute sponge relaxation:
+!   -------------------------
+
+      pext(1:iip1,1:jjp1)=ps(1:iip1,1:jjp1)*aire(1:iip1,1:jjp1)
+
+      do j=1,jjp1
+         ptot(j)=sum(pext(1:iim,j))
+      enddo
+ 
+!   potential temperature
+      do l=l0,llm
+         do j=1,jjp1
+            ! compute zonal average
+            hm=0.
+            do i=1,iim
+               hm=hm+h(i,j,l)*pext(i,j)
+            enddo
+            hm=hm/ptot(j)
+            ! update h()
+            do i=1,iim
+               h(i,j,l)=h(i,j,l)-cst(l)*(h(i,j,l)-hm)
+            enddo
+            h(iip1,j,l)=h(1,j,l)
+         enddo
+      enddo
+
+!   zonal wind
+      do l=l0,llm
+         do j=2,jjm
+            um=0.
+            if(mode.ge.1) then ! compute zonal average
+               do i=1,iim
+                  um=um+0.5*ucov(i,j,l)*(pext(i,j)+pext(i+1,j))/cu(i,j)
+               enddo
+               um=um/ptot(j)
+            endif
+            ! update ucov()
+            do i=1,iim
+               ucov(i,j,l)=ucov(i,j,l)-cst(l)*(ucov(i,j,l)-um*cu(i,j))
+            enddo
+            ucov(iip1,j,l)=ucov(1,j,l)
+         enddo
+      enddo
+
+!  meridional wind
+      do l=l0,llm
+         do j=1,jjm
+            vm=0.
+            if(mode.ge.2) then ! compute zonal average
+               do i=1,iim
+                  vm=vm+vcov(i,j,l)*(pext(i,j)+pext(i,j+1))/cv(i,j)
+               enddo
+               vm=vm/(ptot(j)+ptot(j+1))
+            endif
+            ! update vcov
+            do i=1,iim
+               vcov(i,j,l)=vcov(i,j,l)-cst(l)*(vcov(i,j,l)-vm*cv(i,j))
+            enddo
+            vcov(iip1,j,l)=vcov(1,j,l)
+         enddo
+      enddo
+
+      end subroutine sponge
+      
+end module sponge_mod
+
Index: trunk/LMDZ.COMMON/libf/dyn3d/vlspltqs.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/vlspltqs.F	(revision 1012)
+++ trunk/LMDZ.COMMON/libf/dyn3d/vlspltqs.F	(revision 1017)
@@ -22,4 +22,5 @@
 c     pk exner au milieu des couches necessaire pour calculer Qsat
 c   --------------------------------------------------------------------
+      use cpdet_mod, only: tpot2t
       IMPLICIT NONE
 c
