Index: /trunk/DOC/chantiers/commit_importants.log
===================================================================
--- /trunk/DOC/chantiers/commit_importants.log	(revision 1016)
+++ /trunk/DOC/chantiers/commit_importants.log	(revision 1017)
@@ -1233,2 +1233,17 @@
   a Forward Euler approximation.
 - Sponge tendencies are now added within top_bound.
+
+**********************
+**** commit_v1017 ****
+**********************
+Ehouarn: Improved cpdet routines in // and additional sponge mode:
+- Additionnal sponge mode (trigered with "callsponge" flag), in line with the one used
+  in the Generic and Martian GCM. This sponge is called whenever there is a dissipation
+  step.
+- Improvement of the // cpdet routines : created routines tpot2t_glo_p and t2tpot_glo_p
+  which handle fields on the whole dynamics (scaler) grid, which are more efficient than
+  calling tpot2t_p or t2tpot_p with slabs of data (generated use of intermediate copies
+  of these chunks of data at every call)
+- Turned cpdet.F into a module cpdet_mod.F90 (and correspondingly adapted all routines
+  in the Venus physics).
+
Index: /trunk/LMDZ.COMMON/libf/dyn3d/bilan_dyn.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/bilan_dyn.F	(revision 1016)
+++ /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 1016)
+++ /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 1016)
+++ /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: unk/LMDZ.COMMON/libf/dyn3d/cpdet.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/cpdet.F	(revision 1016)
+++ 	(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 1016)
+++ /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 1016)
+++ /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 1016)
+++ /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
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 1016)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 1017)
@@ -39,4 +39,5 @@
       USE infotrac
       USE control_mod
+      USE cpdet_mod, only: tpot2t_p, t2tpot_p
 
       IMPLICIT NONE
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1016)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1017)
@@ -19,4 +19,5 @@
       USE infotrac, ONLY : type_trac
       use assert_m, only: assert
+      use sponge_mod_p, only: callsponge,mode_sponge,nsponge,tetasponge
       IMPLICIT NONE
 c-----------------------------------------------------------------------
@@ -388,4 +389,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: unk/LMDZ.COMMON/libf/dyn3dpar/cpdet.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/cpdet.F	(revision 1016)
+++ 	(revision )
@@ -1,199 +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 t2tpot_p(nlon,nlev, yt, yteta, ypk)
-! Parallel version of t2tpot
-      USE parallel
-      USE control_mod, only : planet_type
-      IMPLICIT none
-
-! for cpp, nu_venus and t0_venus:
-#include "comconst.h"
-
-      integer,intent(in) :: nlon,nlev
-      real,intent(in) :: yt(nlon,nlev)
-      real,intent(out) :: yteta(nlon,nlev)
-      real,intent(in) :: ypk(nlon,nlev)
-! local variable:
-      integer :: l
-      
-      if (planet_type.eq."venus") then
-        do l=1,nlev
-          yteta(:,l)=yt(:,l)**nu_venus                                  &
-     &                     -nu_venus*t0_venus**nu_venus*                &
-     &                          log(ypk(:,l)/cpp)
-          yteta(:,l)=yteta(:,l)**(1./nu_venus)
-        enddo
-      else
-        do l=1,nlev
-          yteta(:,l)=yt(:,l)*cpp/ypk(:,l)
-        enddo
-      endif ! of if (planet_type.eq."venus")
-
-      END
-
-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======================================================================
-      SUBROUTINE tpot2t_p(nlon,nlev,yteta,yt,ypk)
-! Parallel version of tpot2t
-      USE parallel
-      USE control_mod, only : planet_type
-      IMPLICIT none
-! for cpp, nu_venus and t0_venus:
-#include "comconst.h"
-
-      integer,intent(in) :: nlon,nlev
-      real,intent(out) :: yt(nlon,nlev)
-      real,intent(in) :: yteta(nlon,nlev)
-      real,intent(in) :: ypk(nlon,nlev)
-
-! local variable:
-      integer :: l
-
-      if (planet_type.eq."venus") then
-        do l=1,nlev
-          yt(:,l)=yteta(:,l)**nu_venus                                  &
-     &                  +nu_venus*t0_venus**nu_venus*                   &
-     &                       log(ypk(:,l)/cpp)
-          yt(:,l)=yt(:,l)**(1./nu_venus)
-        enddo
-      else
-        do l=1,nlev
-          yt(:,l)=yteta(:,l)*ypk(:,l)/cpp
-        enddo
-      endif ! of if (planet_type.eq."venus")
-      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/dyn3dpar/cpdet_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/cpdet_mod.F90	(revision 1017)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/cpdet_mod.F90	(revision 1017)
@@ -0,0 +1,303 @@
+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 t2tpot_p(nlon,nlev, yt, yteta, ypk)
+! Parallel version of t2tpot, for an arbitrary number of columns
+      USE control_mod, only : planet_type
+      IMPLICIT none
+
+! for cpp, nu_venus and t0_venus:
+#include "comconst.h"
+
+      integer,intent(in) :: nlon,nlev
+      real,intent(in) :: yt(nlon,nlev)
+      real,intent(out) :: yteta(nlon,nlev)
+      real,intent(in) :: ypk(nlon,nlev)
+! local variable:
+      integer :: l
+      
+      if (planet_type.eq."venus") then
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,nlev
+          yteta(:,l)=yt(:,l)**nu_venus                                  &
+     &                     -nu_venus*t0_venus**nu_venus*                &
+     &                          log(ypk(:,l)/cpp)
+          yteta(:,l)=yteta(:,l)**(1./nu_venus)
+        enddo
+!$OMP END DO
+      else
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,nlev
+          yteta(:,l)=yt(:,l)*cpp/ypk(:,l)
+        enddo
+!$OMP END DO
+      endif ! of if (planet_type.eq."venus")
+
+      end subroutine t2tpot_p
+
+!======================================================================
+!======================================================================
+
+      SUBROUTINE t2tpot_glo_p(yt, yteta, ypk)
+! Parallel version of t2tpot, over the full dynamics (scalar) grid
+! (more efficient than multiple calls to t2tpot_p() with slices of data)
+      USE parallel, only : jj_begin,jj_end
+      USE control_mod, only : planet_type
+      IMPLICIT none
+! for iip1, jjp1 and llm
+#include "dimensions.h"
+#include "paramet.h"
+! for cpp, nu_venus and t0_venus:
+#include "comconst.h"
+
+      real,intent(in) :: yt(iip1,jjp1,llm)
+      real,intent(out) :: yteta(iip1,jjp1,llm)
+      real,intent(in) :: ypk(iip1,jjp1,llm)
+! local variable:
+      integer :: j,l
+      integer :: jjb,jje
+      
+      jjb=jj_begin
+      jje=jj_end
+
+      if (planet_type.eq."venus") then
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,llm
+          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)**nu_venus                  &
+     &                     -nu_venus*t0_venus**nu_venus*                &
+     &                          log(ypk(:,jjb:jje,l)/cpp)
+          yteta(:,jjb:jje,l)=yteta(:,jjb:jje,l)**(1./nu_venus)
+        enddo
+!$OMP END DO
+      else
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,llm
+          yteta(:,jjb:jje,l)=yt(:,jjb:jje,l)*cpp/ypk(:,jjb:jje,l)
+        enddo
+!$OMP END DO
+      endif ! of if (planet_type.eq."venus")
+
+      end subroutine t2tpot_glo_p
+
+!======================================================================
+
+      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
+
+!======================================================================
+!======================================================================
+
+      SUBROUTINE tpot2t_p(nlon,nlev,yteta,yt,ypk)
+! Parallel version of tpot2t, for an arbitrary number of columns
+      USE control_mod, only : planet_type
+      IMPLICIT none
+! for cpp, nu_venus and t0_venus:
+#include "comconst.h"
+
+      integer,intent(in) :: nlon,nlev
+      real,intent(out) :: yt(nlon,nlev)
+      real,intent(in) :: yteta(nlon,nlev)
+      real,intent(in) :: ypk(nlon,nlev)
+
+! local variable:
+      integer :: l
+
+      if (planet_type.eq."venus") then
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,nlev
+          yt(:,l)=yteta(:,l)**nu_venus                                  &
+     &                  +nu_venus*t0_venus**nu_venus*                   &
+     &                       log(ypk(:,l)/cpp)
+          yt(:,l)=yt(:,l)**(1./nu_venus)
+        enddo
+!$OMP END DO
+      else
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,nlev
+          yt(:,l)=yteta(:,l)*ypk(:,l)/cpp
+        enddo
+!$OMP END DO
+      endif ! of if (planet_type.eq."venus")
+      end subroutine tpot2t_p
+
+!======================================================================
+!======================================================================
+
+      SUBROUTINE tpot2t_glo_p(yteta,yt,ypk)
+! Parallel version of tpot2t, over the full dynamics (scalar) grid
+! (more efficient than multiple calls to tpot2t_p() with slices of data)
+      USE parallel, only : jj_begin,jj_end
+      USE control_mod, only : planet_type
+      IMPLICIT none
+! for iip1, jjp1 and llm
+#include "dimensions.h"
+#include "paramet.h"
+! for cpp, nu_venus and t0_venus:
+#include "comconst.h"
+
+      real,intent(out) :: yt(iip1,jjp1,llm)
+      real,intent(in) :: yteta(iip1,jjp1,llm)
+      real,intent(in) :: ypk(iip1,jjp1,llm)
+! local variable:
+      integer :: j,l
+      integer :: jjb,jje
+      
+      jjb=jj_begin
+      jje=jj_end
+
+      if (planet_type.eq."venus") then
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,llm
+          yt(:,jjb:jje,l)=yteta(:,jjb:jje,l)**nu_venus                  &
+     &                  +nu_venus*t0_venus**nu_venus*                   &
+     &                       log(ypk(:,jjb:jje,l)/cpp)
+          yt(:,jjb:jje,l)=yt(:,jjb:jje,l)**(1./nu_venus)
+        enddo
+!$OMP END DO
+      else
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,llm
+          yt(:,jjb:jje,l)=yteta(:,jjb:jje,l)*ypk(:,jjb:jje,l)/cpp
+        enddo
+!$OMP END DO
+      endif ! of if (planet_type.eq."venus")
+      end subroutine tpot2t_glo_p
+
+!======================================================================
+!======================================================================
+!
+! ATTENTION
+!
+! Si un jour on a besoin, il faudra coder les routines 
+!    dt2dtpot / dtpto2dt 
+!
+!======================================================================
+!======================================================================
+end module cpdet_mod
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1016)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1017)
@@ -19,4 +19,5 @@
       USE filtreg_mod
       USE control_mod
+      use cpdet_mod, only: ini_cpdet
 
 ! Ehouarn: the following are needed with (parallel) physics:
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1016)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1017)
@@ -21,4 +21,6 @@
        USE getparam
        USE control_mod
+       use cpdet_mod, only: cpdet,tpot2t_glo_p,t2tpot_glo_p
+       use sponge_mod_p, only: callsponge,mode_sponge,sponge_p
 
       IMPLICIT NONE
@@ -186,5 +188,4 @@
       ! for CP(T)  -- Aymeric
       real :: dtec
-      real,external :: cpdet
       real,save :: ztetaec(ip1jmp1,llm)  !!SAVE ???
 
@@ -211,5 +212,7 @@
 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
 
@@ -550,8 +553,7 @@
 c   --------------------------------
 ! ADAPTATION GCM POUR CP(T)
+      call tpot2t_glo_p(teta,temp,pk)
       ijb=ij_begin
       ije=ij_end
-      call tpot2t_p(ije-ijb+1,llm,teta(ijb:ije,:),temp(ijb:ije,:),
-     &                            pk(ijb:ije,:))
 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       do l=1,llm
@@ -1158,6 +1160,5 @@
 
       IF(apdiss) THEN
-cc$OMP  PARALLEL DEFAULT(SHARED) 
-cc$OMP+          PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
+
 c$OMP MASTER
         call suspend_timer(timer_caldyn)
@@ -1169,4 +1170,22 @@
         call VTb(VThallo)
 c$OMP END MASTER
+
+        ! sponge layer
+        if (callsponge) then
+          call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
+     &                          jj_Nb_dissip,1,1,Request_dissip)
+          call Register_Hallo(ps,ip1jm,1,1,1,1,1,Request_Dissip)
+          call SendRequest(Request_Dissip)
+c$OMP BARRIER
+          call WaitRequest(Request_Dissip)
+c$OMP BARRIER
+c$OMP MASTER
+          call VTe(VThallo)
+          call VTb(VThallo)
+c$OMP END MASTER
+c$OMP BARRIER
+          CALL sponge_p(ucov,vcov,teta,ps,dtdiss,mode_sponge)
+        endif
+
 
 c$OMP BARRIER
@@ -1205,8 +1224,5 @@
 c   dissipation
 ! ADAPTATION GCM POUR CP(T)
-            ijb=ij_begin
-            ije=ij_end
-        call tpot2t_p(ije-ijb+1,llm,teta(ijb:ije,:),temp(ijb:ije,:),
-     &                            pk(ijb:ije,:))
+        call tpot2t_glo_p(teta,temp,pk)
 
 !        CALL FTRACE_REGION_BEGIN("dissip")
@@ -1266,6 +1282,7 @@
             enddo
 c$OMP END DO 
-        call t2tpot_p(ije-ijb+1,llm,temp(ijb:ije,:),ztetaec(ijb:ije,:),
-     &                            pk(ijb:ije,:))
+!        call t2tpot_p(ije-ijb+1,llm,temp(ijb:ije,:),ztetaec(ijb:ije,:),
+!     &                            pk(ijb:ije,:))
+         call t2tpot_glo_p(temp,ztetaec,pk)
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)            
             do l=1,llm
@@ -1576,8 +1593,7 @@
 
 ! ADAPTATION GCM POUR CP(T)
+      call tpot2t_glo_p(teta,temp,pk)
       ijb=ij_begin
       ije=ij_end
-      call tpot2t_p(ije-ijb+1,llm,teta(ijb:ije,:),temp(ijb:ije,:),
-     &                            pk(ijb:ije,:))
 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       do l=1,llm
@@ -1801,8 +1817,7 @@
 
 ! ADAPTATION GCM POUR CP(T)
+                call tpot2t_glo_p(teta,temp,pk)
                 ijb=ij_begin
                 ije=ij_end
-           call tpot2t_p(ije-ijb+1,llm,teta(ijb:ije,:),temp(ijb:ije,:),
-     &                            pk(ijb:ije,:))
 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
                 do l=1,llm      
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/sponge_mod_p.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/sponge_mod_p.F90	(revision 1017)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/sponge_mod_p.F90	(revision 1017)
@@ -0,0 +1,188 @@
+module sponge_mod_p
+
+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_p(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))
+
+      USE Write_Field_p
+      use parallel, only: pole_sud,pole_nord,jj_begin,jj_end
+      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 :: jjb,jje
+      integer,save :: l0 ! layer down to which sponge is applied
+
+      real ssum
+
+      real zkm
+      logical,save :: firstcall=.true.
+
+
+
+      if (firstcall) then
+!$OMP BARRIER
+!$OMP MASTER
+       ! 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.
+!$OMP END MASTER
+!$OMP BARRIER
+      endif ! of if (firstcall)
+
+!-----------------------------------------------------------------------
+!   compute sponge relaxation:
+!   -------------------------
+      jjb=jj_begin
+      jje=jj_end+1 ! +1 because vcov(j) computations require pext(j+1) & ptot(j+1)
+      IF (pole_sud)  jje=jj_end-1+1
+!$OMP MASTER
+      pext(1:iip1,jjb:jje)=ps(1:iip1,jjb:jje)*aire(1:iip1,jjb:jje)
+      do j=jjb,jje
+         ptot(j)=sum(pext(1:iim,j))
+      enddo
+!$OMP END MASTER
+!$OMP BARRIER
+
+
+!   potential temperature
+      jjb=jj_begin
+      jje=jj_end
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+      do l=l0,llm
+         do j=jjb,jje
+            ! 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
+!$OMP END DO NOWAIT
+
+!   zonal wind
+      jjb=jj_begin
+      jje=jj_end
+      IF (pole_nord) jjb=jj_begin+1
+      IF (pole_sud)  jje=jj_end-1
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+      do l=l0,llm
+         do j=jjb,jje
+            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
+!$OMP END DO NOWAIT
+
+!  meridional wind
+      jjb=jj_begin
+      jje=jj_end
+      IF (pole_sud) jje=jj_end-1
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+      do l=l0,llm
+         do j=jjb,jje
+            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
+!$OMP END DO
+
+      end subroutine sponge_p
+      
+end module sponge_mod_p
+
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltqs_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltqs_p.F	(revision 1016)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltqs_p.F	(revision 1017)
@@ -25,4 +25,5 @@
       USE mod_hallo
       USE VAMPIR
+      use cpdet_mod, only: tpot2t_glo_p
       IMPLICIT NONE
 
@@ -97,6 +98,6 @@
 ! ADAPTATION GCM POUR CP(T)
 ! probablement a revoir...
-      call tpot2t_p(ip1jmp1,llm,teta,tempe,pk)
-
+!      call tpot2t_p(ip1jmp1,llm,teta,tempe,pk)
+      call tpot2t_glo_p(teta,tempe,pk)
 
       call SetTag(MyRequest1,100)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/YOMCST.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/YOMCST.h	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/YOMCST.h	(revision 1017)
@@ -20,6 +20,6 @@
       REAL RKAPPA,RETV
 ! ADAPTATION GCM POUR CP(T)
-      real cpdet
-      external cpdet
+!      real cpdet
+!      external cpdet
 ! A1.5,6 Thermodynamic liquid,solid phases
       REAL RCW,RCS
Index: /trunk/LMDZ.VENUS/libf/phyvenus/ajsec.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/ajsec.F	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/ajsec.F	(revision 1017)
@@ -7,4 +7,5 @@
 
       use dimphy
+      use cpdet_mod, only: t2tpot, tpot2t
       IMPLICIT none
 c======================================================================
Index: /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F	(revision 1017)
@@ -35,4 +35,5 @@
       USE interface_surf
       use dimphy
+      use cpdet_mod, only: t2tpot
       IMPLICIT none
 c======================================================================
@@ -490,4 +491,5 @@
       USE interface_surf
       use dimphy
+      use cpdet_mod, only: t2tpot,tpot2t,cpdet
 
       IMPLICIT none
@@ -898,4 +900,5 @@
 
       use dimphy
+      use cpdet_mod, only: cpdet,t2tpot
       IMPLICIT none
 c======================================================================
@@ -1181,4 +1184,5 @@
 
       use dimphy
+      use cpdet_mod, only: cpdet
       IMPLICIT none
 c======================================================================
Index: /trunk/LMDZ.VENUS/libf/phyvenus/comcstfi.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/comcstfi.h	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/comcstfi.h	(revision 1017)
@@ -5,6 +5,4 @@
 
       REAL rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
-      real cpdet
-      external cpdet
 
 !-----------------------------------------------------------------------
Index: /trunk/LMDZ.VENUS/libf/phyvenus/diagphy.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/diagphy.F	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/diagphy.F	(revision 1017)
@@ -209,4 +209,5 @@
  
       use dimphy
+      use cpdet_mod, only: cpdet
       IMPLICIT NONE
 C
Index: /trunk/LMDZ.VENUS/libf/phyvenus/interface_surf.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/interface_surf.F90	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/interface_surf.F90	(revision 1017)
@@ -46,4 +46,5 @@
 
       use write_field_phy
+      use cpdet_mod, only: cpdet
 
       IMPLICIT none
@@ -196,4 +197,5 @@
 
   use write_field_phy
+  use cpdet_mod, only: t2tpot, tpot2t
 
   IMPLICIT none
Index: /trunk/LMDZ.VENUS/libf/phyvenus/lw_venus_ve.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/lw_venus_ve.F	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/lw_venus_ve.F	(revision 1017)
@@ -6,4 +6,5 @@
       
       use dimphy
+      use cpdet_mod, only: cpdet
       IMPLICIT none
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/lwi.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/lwi.F	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/lwi.F	(revision 1017)
@@ -2,4 +2,5 @@
 
       use dimphy
+      use cpdet_mod, only: cpdet
       implicit none
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1017)
@@ -64,4 +64,5 @@
       USE write_field_phy
       USE iophy
+      use cpdet_mod, only: cpdet, t2tpot
       IMPLICIT none
 c======================================================================
Index: /trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F	(revision 1017)
@@ -1,3 +1,4 @@
-
+#include "../dyn3d/cpdet_mod.F90"
+!***********************************************************************
       PROGRAM rcm1d
       
@@ -6,4 +7,5 @@
       use comgeomphy
       USE phys_state_var_mod
+      use cpdet_mod, only: ini_cpdet
       IMPLICIT NONE
 
@@ -466,4 +468,3 @@
 #include "../dyn3d/abort_gcm.F"
 #include "../dyn3d/dump2d.F"
-#include "../dyn3d/cpdet.F"
-
+
Index: /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_dc.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_dc.F	(revision 1016)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_dc.F	(revision 1017)
@@ -5,4 +5,5 @@
       
       use dimphy
+      use cpdet_mod, only: cpdet
       IMPLICIT none
 
