Index: trunk/LMDZ.GENERIC/libf/dyn3d/comdissip.h
===================================================================
--- trunk/LMDZ.GENERIC/libf/dyn3d/comdissip.h	(revision 993)
+++ trunk/LMDZ.GENERIC/libf/dyn3d/comdissip.h	(revision 1006)
@@ -1,7 +1,7 @@
-c-----------------------------------------------------------------------
-c INCLUDE dissip.h
+!-----------------------------------------------------------------------
+! INCLUDE dissip.h
 
-      COMMON/comdissip/
-     $    niterdis,coefdis,tetavel,tetatemp,gamdissip
+      COMMON/comdissip/                                                 &
+     &    niterdis,coefdis,tetavel,tetatemp,gamdissip
 
 
@@ -10,3 +10,3 @@
       REAL tetavel,tetatemp,coefdis,gamdissip
 
-c-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
Index: trunk/LMDZ.GENERIC/libf/dyn3d/defrun_new.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/dyn3d/defrun_new.F	(revision 993)
+++ trunk/LMDZ.GENERIC/libf/dyn3d/defrun_new.F	(revision 1006)
@@ -37,4 +37,5 @@
 ! to use  'getin'
       USE ioipsl_getincom
+      use sponge_mod,only: callsponge,nsponge,mode_sponge,tetasponge
       IMPLICIT NONE
 
@@ -45,5 +46,4 @@
 #include "serre.h"
 #include "comdissnew.h"
-#include "sponge.h"
 c
 c   arguments:
@@ -451,5 +451,5 @@
 
         WRITE(tapeout,*) ""
-        WRITE(tapeout,*) "Avec sponge layer"
+        WRITE(tapeout,*) "Use a sponge layer?"
         callsponge=.true. ! default value
         call getin("callsponge",callsponge)
@@ -457,6 +457,16 @@
 
         WRITE(tapeout,*) ""
-        WRITE(tapeout,*) "Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), ",
-     & "mode2(u=umoy,v=vmoy)"
+        WRITE(tapeout,*) "Sponge: number of layers over which",
+     &                    " sponge extends"
+        nsponge=3 ! default value
+        call getin("nsponge",nsponge)
+        WRITE(tapeout,*)" nsponge = ",nsponge
+
+        WRITE(tapeout,*)""
+        WRITE(tapeout,*)"Sponge mode: (forcing is towards ..."
+        WRITE(tapeout,*)"  over upper nsponge layers)"
+        WRITE(tapeout,*)"  0: (h=hmean,u=v=0)"
+        WRITE(tapeout,*)"  1: (h=hmean,u=umean,v=0)"
+        WRITE(tapeout,*)"  2: (h=hmean,u=umean,v=vmean)"
         mode_sponge=2 ! default value
         call getin("mode_sponge",mode_sponge)
@@ -464,11 +474,7 @@
 
         WRITE(tapeout,*) ""
-        WRITE(tapeout,*) "Sponge:  hauteur de sponge (km)"
-        hsponge=90.0 ! default value
-        call getin("hsponge",hsponge)
-        WRITE(tapeout,*)" hsponge = ",hsponge
-
-        WRITE(tapeout,*) ""
-        WRITE(tapeout,*) "Sponge:  tetasponge (secondes)"
+        WRITE(tapeout,*) "Sponge: characteristic time scale tetasponge"
+        WRITE(tapeout,*) "(seconds) at topmost layer (time scale then "
+        WRITE(tapeout,*) " doubles with decreasing layer index)."
         tetasponge=50000.0
         call getin("tetasponge",tetasponge)
Index: trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F	(revision 993)
+++ trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F	(revision 1006)
@@ -1,4 +1,5 @@
       PROGRAM gcm
 
+      use sponge_mod,only: callsponge,mode_sponge,sponge
       IMPLICIT NONE
 
@@ -46,5 +47,4 @@
 #include "serre.h"
 #include "tracstoke.h"
-#include "sponge.h"
 #include"advtrac.h"
 
@@ -444,9 +444,8 @@
 c        Sponge layer
 c        ~~~~~~~~~~~~
-         DO ij=1, ip1jmp1
-            pext(ij)=ps(ij)*aire(ij)
-         ENDDO
+
          IF (callsponge) THEN
-            CALL sponge(ucov,vcov,teta,pext,dtdiss,mode_sponge)
+           pext(1:ip1jmp1)=ps(1:ip1jmp1)*aire(1:ip1jmp1)
+           CALL sponge(ucov,vcov,teta,pext,dtdiss,mode_sponge)
          ENDIF
 
Index: trunk/LMDZ.GENERIC/libf/dyn3d/sponge.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/dyn3d/sponge.F	(revision 993)
+++ 	(revision )
@@ -1,137 +1,0 @@
-      subroutine sponge(ucov,vcov,h,pext,dt,mode)
-      implicit none
-#include "dimensions.h"
-#include "paramet.h"
-#include "comdissip.h"
-#include "comvert.h"
-#include "comgeom2.h"
-#include "sponge.h"
-      real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm)
-      real h(iip1,jjp1,llm),pext(iip1,jjp1)
-      real dt
-      real sig_s(llm) !sigma au milieu des couches
-      save sig_s
-
-c   Local:
-c   ------
-
-      REAL vm,um,hm,ptot(jjp1)
-      real cst(llm)
-      integer mode
-
-      INTEGER l,i,j,l0
-
-      real ssum
-
-      real echelle,zkm
-      logical firstcall
-      save firstcall,cst,l0
-      data firstcall/.true./
-
-      if (firstcall) then
-       do l=1,llm
-          sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
-       enddo
-
-       IF (mode.eq.3) then
-         l0=llm-2
-         echelle=10.
- 
-         PRINT*
-         print*,'sponge mode',mode
-         print*,'hsponge',hsponge
-         print*,'tetasponge n intervient pas'
-         print*,'Constantes de dissipations fixees comme les anglais'
-         print*,'Coeffs pour la couche en eponge'
-         print*,'Z (km)  tau'
-         cst(llm)=dt/(0.5*88775)
-         cst(llm-1)=dt/(88775)
-         cst(l0)=dt/(2*88775)
-         do l=l0,llm
-            zkm=-echelle*log(sig_s(l))
-            print*,zkm,dt/cst(l),cst(l)
-         enddo
-         firstcall=.false.
-         PRINT*
-       ELSE
-         l0=1
-         echelle=10.
- 
-         PRINT*
-         print*,'sponge mode',mode
-         print*,'hsponge tetasponge ',hsponge,tetasponge
-         print*,'Coeffs pour la couche en eponge'
-         print*,'Z (km)  tau'
-         do l=1,llm
-            zkm=-echelle*log(sig_s(l))
-            cst(l)=.5*(1.+tanh((zkm-hsponge)/echelle))
-            cst(l)= max(tetasponge*1.e-15,cst(l))
-            print*,zkm,1./cst(l)*tetasponge,cst(l)*dt/tetasponge
-            cst(l)=cst(l)*dt/tetasponge
-         enddo
-         firstcall=.false.
-         PRINT*
-       ENDIF
-      endif
-
-c-----------------------------------------------------------------------
-c   calcul de la dissipation:
-c   -------------------------
-
-      do j=1,jjp1
-         ptot(j)=ssum(iim,pext(1,j),1)
-      enddo
- 
-c   temperature potentielle
-      do l=l0,llm
-         do j=1,jjp1
-            hm=0.
-            do i=1,iim
-               hm=hm+h(i,j,l)*pext(i,j)
-            enddo
-            hm=hm/ptot(j)
-            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
-
-c   vent zonal
-      do l=l0,llm
-         do j=2,jjm
-            um=0.
-            if(mode.ge.1) then
-               do i=1,iim
-                  um=um+0.5*ucov(i,j,l)*(pext(i,j)+pext(i+1,j))
-     s               /cu(i,j)
-               enddo
-               um=um/ptot(j)
-            endif
-            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
-
-c  vent meridien
-      do l=l0,llm
-         do j=1,jjm
-            vm=0.
-            if(mode.ge.2) then
-               do i=1,iim
-                  vm=vm+vcov(i,j,l)*(pext(i,j)+pext(i,j+1))
-     s               /cv(i,j)
-               enddo
-               vm=vm/(ptot(j)+ptot(j+1))
-            endif
-            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
-
-      RETURN
-      end
Index: trunk/LMDZ.GENERIC/libf/dyn3d/sponge.h
===================================================================
--- trunk/LMDZ.GENERIC/libf/dyn3d/sponge.h	(revision 993)
+++ 	(revision )
@@ -1,9 +1,0 @@
-c-----------------------------------------------------------------------
-c INCLUDE 'sponge.h'
-
-      COMMON/com_sponge/callsponge,mode_sponge,hsponge,tetasponge
-
-      LOGICAL   callsponge
-      INTEGER   mode_sponge
-      REAL   hsponge,tetasponge
-c-----------------------------------------------------------------------
Index: trunk/LMDZ.GENERIC/libf/dyn3d/sponge_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/dyn3d/sponge_mod.F90	(revision 1006)
+++ trunk/LMDZ.GENERIC/libf/dyn3d/sponge_mod.F90	(revision 1006)
@@ -0,0 +1,157 @@
+module sponge_mod
+
+implicit none
+
+! sponge parameters (set/read via defrun_new.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,pext,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*t))
+
+
+      implicit none
+#include "dimensions.h"
+#include "paramet.h"
+#include "comdissip.h"
+#include "comvert.h"
+#include "comgeom2.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) :: 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)
+
+      INTEGER l,i,j
+      integer,save :: l0 ! layer down to which sponge is applied
+
+      real ssum
+
+      real echelle,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
+ 
+        PRINT*
+        print*,'sponge mode',mode
+        print*,'nsponge tetasponge ',nsponge,tetasponge
+        print*,'Coeffs for the sponge layer'
+        print*,'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
+
+        echelle=10.
+        do l=l0,llm
+           zkm=-echelle*log(sig_s(l))
+           print*,zkm,tetasponge*2**(llm-l),cst(l)
+        enddo
+        PRINT*
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+
+!-----------------------------------------------------------------------
+!   compute sponge relaxation:
+!   -------------------------
+
+      do j=1,jjp1
+         ptot(j)=ssum(iim,pext(1,j),1)
+      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
