Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/callkeys.h
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/callkeys.h	(revision 309)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/callkeys.h	(revision 310)
@@ -6,5 +6,5 @@
 
       COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj      &
-     &   , co2cond,callsoil                                             &
+     &   , co2cond,callsoil,callradsponge                               &
      &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
      &   , callstats,calleofdump                                        &
@@ -50,4 +50,5 @@
      &   , strictboundcorrk                                             
 
+      logical callradsponge
       logical enertest
       logical nonideal
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/inifis.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/inifis.F	(revision 309)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/inifis.F	(revision 310)
@@ -204,4 +204,10 @@
          if (is_master) write(*,*) "Rmean = ", Rmean
          
+         if (is_master) write(*,*) 
+     &      "Add sponge after radiative transfer?"
+         callradsponge=.false.
+         call getin_p("callradsponge", callradsponge)
+         if (is_master) write(*,*) "callradsponge = ", callradsponge
+
 ! Test of incompatibility:
 ! if tlocked, then diurnal should be false
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90	(revision 309)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90	(revision 310)
@@ -280,4 +280,5 @@
       real zdvmolvis(ngrid,nlayer)
       real zdqmoldiff(ngrid,nlayer,nq)
+      real dtradsponge(ngrid,nlayer) ! (K/s)
 
 !     Local variables for local calculations:
@@ -1088,4 +1089,9 @@
 !-------------------------
 
+        if (callradsponge) then
+          call radsponge(ngrid,nlayer,pplay,pplev,pt,pdt,dtphys,dtradsponge)
+          pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtradsponge(1:ngrid,1:nlayer)
+        endif
+
       endif ! of if (callrad)
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/radsponge.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/radsponge.F90	(revision 310)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/radsponge.F90	(revision 310)
@@ -0,0 +1,141 @@
+subroutine radsponge(ngrid,nlayer,pplay,pplev,temp,pdt,dtphys,dtradsponge)
+! add a (vertical) sponge effect on temperature to damp oscillations
+! near the model top (topmost 
+! 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 ioipsl_getincom_p, only: getin_p
+use mod_phys_lmdz_para, only : is_master
+
+implicit none
+
+integer,intent(in) :: ngrid ! number of atmospheric columns
+integer,intent(in) :: nlayer ! number of atmospheric layers
+real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
+real,intent(in) :: pplev(ngrid,nlayer+1) ! interlayer pressure (Pa)
+real,intent(in) :: temp(ngrid,nlayer) ! temperature (K)
+real,intent(in) :: pdt(ngrid,nlayer) ! tendecies on temperature from previous processes (K/s)
+real,intent(in) :: dtphys ! physics time step (s)
+real,intent(out) :: dtradsponge(ngrid,nlayer) ! tendencies on temperature from sponge
+
+! local variables:
+logical,save :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall)
+integer,save :: iflag_sponge   !0 --> for no sponge
+			       !1 --> for sponge over 4 topmost layers
+			       !2 --> for sponge from top to ~1% of top layer pressure
+                               !3 --> for sponge over last nb_sponge_layers
+real,save :: tau_sponge !inverse of charactericstic relaxation time scale at the topmost layer (Hz)
+real,save :: nb_sponge_layers ! number of layers over which the sponge extends
+!$OMP THREADPRIVATE(iflag_sponge,tau_sponge,nb_sponge_layers)
+integer,save :: nlayer0 ! layer down to which sponge is applied
+real,allocatable,save :: rdamp(:) ! quenching coefficient
+real,allocatable,save :: lambda(:) ! inverse or quenching time scale (Hz)
+real,allocatable,save :: mass(:) ! g*mass (per m2) of atmospheric layer
+real,save :: spongemass ! (pseudo-) mass of atmosphere over the entire sponge layer
+
+real :: temp_cur(ngrid,nlayer) ! temperature (K)
+real :: mean_temp(ngrid) ! mean temperature towards which to relax
+
+integer :: i,k
+
+if (firstcall) then
+  ! load sponge parameters
+  iflag_sponge=0 ! default: no sponge
+  call getin_p("iflag_sponge",iflag_sponge)
+  tau_sponge=1.E-4 ! default value
+  call getin_p("tau_sponge",tau_sponge)
+  if (iflag_sponge==3) then
+    call getin_p("nb_sponge_layers",nb_sponge_layers)
+  endif
+
+!$OMP MASTER
+  allocate(rdamp(nlayer))
+  allocate(lambda(nlayer))
+  allocate(mass(nlayer))
+  
+  if (iflag_sponge == 1) then
+! sponge quenching over the topmost 4 atmospheric layers
+    nb_sponge_layers=4
+    nlayer0=nlayer-3
+    lambda(:)=0.
+    lambda(nlayer)=tau_sponge
+    lambda(nlayer-1)=tau_sponge/2.
+    lambda(nlayer-2)=tau_sponge/4.
+    lambda(nlayer-3)=tau_sponge/8.
+  else if (iflag_sponge == 2) then
+! sponge quenching over topmost layers down to pressures which are
+! higher than 100 times the topmost layer pressure
+    do k=nlayer,1,-1
+      lambda(k)=tau_sponge*max(pplay(1,nlayer)/pplay(1,k)-0.01,0.)
+      if (lambda(k).ne.0.) nlayer0=k
+    enddo
+  else if (iflag_sponge == 3) then
+    nlayer0=nlayer-nb_sponge_layers+1
+    lambda(:)=0
+    do k=nlayer,nlayer0,-1
+      lambda(k)=tau_sponge/(2.**(nlayer-k))
+    enddo
+  else
+    write(*,*) "radsponge: bad value for iflag_sponge:",iflag_sponge
+    stop
+  endif
+
+  rdamp(:)=1.-exp(-lambda(:)*dtphys)
+  
+  ! compute (pseudo-) atmospheric mass of each layer
+  ! (here we assume that they won't vay with time,
+  ! which is always the case near model top)
+  do k=1,nlayer
+    mass(k)=pplev(1,k)-pplev(1,k+1)
+  enddo
+  spongemass=0
+  do k=nlayer0,nlayer
+    spongemass=spongemass+mass(k)
+  enddo
+  
+  if (is_master) then
+    write(*,*) "radsponge: nlayer0=",nlayer0," spongemass=",spongemass
+    write(*,*) " layer    lambda       rdamp        mass"
+    do k=nlayer0,nlayer
+      write(*,'(i3,3x,3(1pe14.7,1x))') k,lambda(k),rdamp(k),mass(k)
+    enddo
+  endif
+!$OMP END MASTER
+!$OMP BARRIER
+  
+  firstcall=.false.
+endif ! of if (firstcall)
+
+
+! 1. compute current temperature and average teperature toward which to relax
+
+temp_cur(:,:)=temp(:,:)+dtphys*pdt(:,:)
+
+! compute mean temperature as a mass-weighed average
+mean_temp(1:ngrid)=0
+do k=nlayer0,nlayer
+  mean_temp(1:ngrid)=mean_temp(1:ngrid)+temp_cur(1:ngrid,k)*mass(k)
+enddo
+mean_temp(1:ngrid)=mean_temp(1:ngrid)/spongemass
+
+! 2. Temperature tendency due to sponge
+dtradsponge(:,:)=0
+do k=nlayer0,nlayer
+  ! A(t+dt)-A(t)=-(A(t)-Am)*(1-exp(-lambda*dt))
+  dtradsponge(1:ngrid,k)=-(temp_cur(1:ngrid,k)-mean_temp(1:ngrid))*rdamp(k)
+enddo
+! Ehouarn debug:
+!write(*,*)"radsponge: mean_temp(:)=",mean_temp(:)
+!write(*,*)" layer    temp     temp_cur     dtradsponge"
+!do k=nlayer0,nlayer
+!  write(*,'(i3,3x,3(1pe14.7,1x))') k,temp(1,k),temp_cur(1,k),dtradsponge(1,k)
+!enddo
+
+! return dtradsponge as a tendency in  K/s
+dtradsponge(1:ngrid,nlayer0:nlayer)=dtradsponge(1:ngrid,nlayer0:nlayer)/dtphys
+
+end subroutine radsponge
Index: codes/icosagcm/branches/SATURN_DYNAMICO/TEST/callphys.def
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/TEST/callphys.def	(revision 309)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/TEST/callphys.def	(revision 310)
@@ -46,4 +46,6 @@
 #iradia     = 160
 iradia     = 40
+# Add temperature sponge effect after radiative transfer?
+callradsponge=.true.
 # call multilayer correlated-k radiative transfer ?
 corrk      = .true.
