Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/CREATE_amip/COMP/dynamico.card
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/CREATE_amip/COMP/dynamico.card	(revision 6361)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/CREATE_amip/COMP/dynamico.card	(revision 6362)
@@ -29,4 +29,5 @@
    	(${SUBMIT_DIR}/PARAM/run_dynamico.def_nbp${RESOL_NBP}, run_dynamico.def),\
    	(${SUBMIT_DIR}/PARAM/iodef.xml, .),\
+	(${SUBMIT_DIR}/PARAM/sponge_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_input_dynamico.xml, .),\
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/CREATE_hadisst/COMP/dynamico.card
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/CREATE_hadisst/COMP/dynamico.card	(revision 6361)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/CREATE_hadisst/COMP/dynamico.card	(revision 6362)
@@ -29,4 +29,5 @@
    	(${SUBMIT_DIR}/PARAM/run_dynamico.def_nbp${RESOL_NBP}, run_dynamico.def),\
    	(${SUBMIT_DIR}/PARAM/iodef.xml, .),\
+	(${SUBMIT_DIR}/PARAM/sponge_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_input_dynamico.xml, .),\
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/aquaplanet/CREATE_aquaplanet/COMP/dynamico.card
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/aquaplanet/CREATE_aquaplanet/COMP/dynamico.card	(revision 6361)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/aquaplanet/CREATE_aquaplanet/COMP/dynamico.card	(revision 6362)
@@ -29,4 +29,5 @@
    	(${SUBMIT_DIR}/PARAM/run_dynamico.def_nbp${RESOL_NBP}, run_dynamico.def),\
    	(${SUBMIT_DIR}/PARAM/iodef.xml, .),\
+	(${SUBMIT_DIR}/PARAM/sponge_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_input_dynamico.xml, .),\
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/aquaplanet/climAQUA_360d/COMP/dynamico.card
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/aquaplanet/climAQUA_360d/COMP/dynamico.card	(revision 6361)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/aquaplanet/climAQUA_360d/COMP/dynamico.card	(revision 6362)
@@ -29,4 +29,5 @@
    	(${SUBMIT_DIR}/PARAM/run_dynamico.def_nbp${RESOL_NBP}, run_dynamico.def),\
    	(${SUBMIT_DIR}/PARAM/iodef.xml, .),\
+	(${SUBMIT_DIR}/PARAM/sponge_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_input_dynamico.xml, .),\
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/clim_noleap/COMP/dynamico.card
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/clim_noleap/COMP/dynamico.card	(revision 6361)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/EXPERIMENTS/ICOLMDZ/clim_noleap/COMP/dynamico.card	(revision 6362)
@@ -38,4 +38,5 @@
    	(${SUBMIT_DIR}/PARAM/run_dynamico.def_nbp${RESOL_NBP}, run_dynamico.def),\
    	(${SUBMIT_DIR}/PARAM/iodef.xml, .),\
+	(${SUBMIT_DIR}/PARAM/sponge_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_input_dynamico.xml, .),\
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/DRIVER/dynamico.driver
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/DRIVER/dynamico.driver	(revision 6361)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/DRIVER/dynamico.driver	(revision 6362)
@@ -109,4 +109,9 @@
     # In iodef.xml add on the next line after "COMPONENT CONTEXT"
     echo '<context id="icosagcm" src="./context_dynamico.xml"/>' >> add.tmp
+
+   # Add inclusion of file sponge_dynamico.xml if this file exists
+   if [ -f sponge_dynamico.xml ] ; then
+        echo '<context id="icosagcm"        src="./sponge_dynamico.xml" />' >> add.tmp
+   fi
 
     # Add inclusion of file context_input_dynamico.xml if this file exists
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/PARAM/sponge_dynamico.xml
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/PARAM/sponge_dynamico.xml	(revision 6362)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/PARAM/sponge_dynamico.xml	(revision 6362)
@@ -0,0 +1,111 @@
+<!-- =========================================================================================================== -->
+<!-- nudging_dynamico.xml                                                                                           -->
+<!-- Definition of sponge layer zonal mean                                          -->
+<!-- =========================================================================================================== -->
+<context id="icosagcm" >
+
+
+   <axis_definition>
+     
+     <axis id="reduced_sponge" > <!-- n_glo is set from model -->
+       <reduce_domain  operation="average" direction="iDir"/>
+     </axis>
+
+     <axis id="lev_sponge" axis_type="Z" /> <!-- n_glo is set from model -->
+     
+    </axis_definition>
+    
+    <grid_definition>
+
+      <grid id="sponge_i">
+        <domain domain_ref="i" />
+        <axis axis_ref="lev_sponge"/>
+      </grid>
+   
+      <grid id="to_sponge">
+        <domain id="to_sponge" type="rectilinear" > <!-- ni_glo and nj_glo are set from model -->
+          <generate_rectilinear_domain />
+          <interpolate_domain order="1"/>
+        </domain> 
+        <axis axis_ref="lev_sponge"/>
+      </grid>
+         
+      <grid id="reduced_sponge">
+        <axis axis_ref="reduced_sponge"/>
+        <axis axis_ref="lev_sponge"/>
+      </grid>
+      
+      <grid id="from_sponge">
+        <domain domain_ref="i">
+          <interpolate_domain order="1"/>
+        </domain> 
+        <axis axis_ref="lev_sponge"/>
+      </grid>
+      
+    </grid_definition>
+    
+    <field_definition>
+      <field id="ulon_sponge"                                        grid_ref="sponge_i"/>
+      <field id="ulon_sponge_interp" field_ref="ulon_sponge"         grid_ref="to_sponge" />
+      <field id="ulon_sponge_reduced" field_ref="ulon_sponge_interp" grid_ref="reduced_sponge" read_access="false"/>
+      <field id="ulon_sponge_mean"                                   grid_ref="to_sponge" />
+      <field id="ulon_from_sponge" field_ref="ulon_sponge_mean"      grid_ref="from_sponge" read_access="false"/>
+
+      <field id="ulat_sponge"                                        grid_ref="sponge_i"/>
+      <field id="ulat_sponge_interp" field_ref="ulat_sponge"         grid_ref="to_sponge" />
+      <field id="ulat_sponge_reduced" field_ref="ulat_sponge_interp" grid_ref="reduced_sponge" read_access="false"/>
+      <field id="ulat_sponge_mean"                                   grid_ref="to_sponge" />
+      <field id="ulat_from_sponge" field_ref="ulat_sponge_mean"      grid_ref="from_sponge" read_access="false"/>
+
+      <field id="theta_sponge"                                        grid_ref="sponge_i"/>
+      <field id="theta_sponge_interp" field_ref="theta_sponge"         grid_ref="to_sponge" />
+      <field id="theta_sponge_reduced" field_ref="theta_sponge_interp" grid_ref="reduced_sponge" read_access="false"/>
+      <field id="theta_sponge_mean"                                   grid_ref="to_sponge" />
+      <field id="theta_from_sponge" field_ref="theta_sponge_mean"      grid_ref="from_sponge" read_access="false"/>
+
+      <field id="rhodz_sponge"                                        grid_ref="sponge_i"/>
+      <field id="rhodz_sponge_interp" field_ref="rhodz_sponge"         grid_ref="to_sponge" />
+      <field id="rhodz_sponge_reduced" field_ref="rhodz_sponge_interp" grid_ref="reduced_sponge" read_access="false"/>
+
+     <field id="dulon_sponge" domain_ref="i" axis_ref="lev" />
+     <field id="dulat_sponge" domain_ref="i" axis_ref="lev" />
+     <field id="dtheta_sponge" domain_ref="i" axis_ref="lev" />
+ 
+
+    </field_definition>
+
+    <file_definition>
+      <file id="check_sponge" output_freq="1ts" type="one_file" sync_freq="1ts" enabled="false" >
+        <field_group  prec="8">
+          <field field_ref="ulon_sponge" operation="instant" enabled="true"/>
+          <field field_ref="ulon_sponge_interp" operation="instant" enabled="true"/>
+          <field field_ref="ulon_sponge_reduced" operation="instant" enabled="true"/>
+          <field field_ref="ulon_sponge_mean" operation="instant" enabled="true"/>
+          <field field_ref="ulon_from_sponge" operation="instant" enabled="true"/>
+
+          <field field_ref="ulat_sponge" operation="instant" enabled="true"/>
+          <field field_ref="ulat_sponge_interp" operation="instant" enabled="true"/>
+          <field field_ref="ulat_sponge_reduced" operation="instant" enabled="true"/>
+          <field field_ref="ulat_sponge_mean" operation="instant" enabled="true"/>
+          <field field_ref="ulat_from_sponge" operation="instant" enabled="true"/>
+
+          <field field_ref="theta_sponge" operation="instant" enabled="true"/>
+          <field field_ref="theta_sponge_interp" operation="instant" enabled="true"/>
+          <field field_ref="theta_sponge_reduced" operation="instant" enabled="true"/>
+          <field field_ref="theta_sponge_mean" operation="instant" enabled="true"/>
+          <field field_ref="theta_from_sponge" operation="instant" enabled="true"/>
+          
+          <field field_ref="rhodz_sponge" operation="instant" enabled="true"/>
+          <field field_ref="rhodz_sponge_interp" operation="instant" enabled="true"/>
+          <field field_ref="rhodz_sponge_reduced" operation="instant" enabled="true"/>
+
+         <field field_ref="dulon_sponge" operation="instant" enabled="true"/>
+         <field field_ref="dulat_sponge" operation="instant" enabled="true"/>
+         <field field_ref="dtheta_sponge" operation="instant" enabled="true"/>
+
+        </field_group>
+      </file>
+    </file_definition>
+    
+</context>
+
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/PARAM/vert.def_L79
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/PARAM/vert.def_L79	(revision 6361)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/GENERAL/PARAM/vert.def_L79	(revision 6362)
@@ -29,7 +29,7 @@
 disvert=strato_custom
 # sponge parmeters
-iflag_sponge=0
+iflag_sponge=2
 tau_sponge=5.e-4
-mode_sponge=1
+mode_sponge=3
 
 ## Section 3: Following paramereters are shared between LMDZ/dynamics and DYNAMICO
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/dissip/sponge.f90
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/dissip/sponge.f90	(revision 6362)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/dissip/sponge.f90	(revision 6362)
@@ -0,0 +1,497 @@
+MODULE sponge_mod
+  USE icosa
+
+  PRIVATE
+  
+  REAL,SAVE :: tau_sponge !inverse of charactericstic relaxation time scale at the topmost layer (Hz)
+  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
+  INTEGER,SAVE :: mode_sponge !1 --> u and v relax towards 0
+                              !2 --> u and v relax towards their zonal mean
+                              !3 --> u,v and pot. temp. relax towards their zonal mean
+ 
+  !$OMP THREADPRIVATE(tau_sponge,iflag_sponge,mode_sponge)
+  REAL,ALLOCATABLE,SAVE :: rdamp(:) ! quenching coefficient
+  !$OMP THREADPRIVATE(rdamp)
+  REAL,ALLOCATABLE,SAVE:: lambda(:) ! inverse or quenching time scale (Hz)
+  !$OMP THREADPRIVATE(lambda)
+
+  PUBLIC sponge, init_sponge, iflag_sponge, pre_init_sponge
+
+  INTEGER,SAVE :: ni_glo
+  !$OMP THREADPRIVATE(ni_glo)
+  INTEGER,SAVE :: nj_glo
+  !$OMP THREADPRIVATE(nj_glo)
+  INTEGER,SAVE :: ni
+  !$OMP THREADPRIVATE(ni)
+  INTEGER,SAVE :: nj
+  !$OMP THREADPRIVATE(nj)
+  INTEGER,SAVE :: ibegin
+  !$OMP THREADPRIVATE(ibegin)
+  INTEGER,SAVE :: jbegin
+  !$OMP THREADPRIVATE(jbegin)
+ 
+  TYPE(t_field),POINTER,SAVE :: f_hflux(:)
+  TYPE(t_field),POINTER,SAVE :: f_ulon(:)
+  TYPE(t_field),POINTER,SAVE :: f_ulat(:)
+  TYPE(t_field),POINTER,SAVE :: f_ulon_sponge(:)
+  TYPE(t_field),POINTER,SAVE :: f_ulat_sponge(:)
+  TYPE(t_field),POINTER,SAVE :: f_rhodz_sponge(:)
+  TYPE(t_field),POINTER,SAVE :: f_theta_sponge(:)
+
+  TYPE(t_field),POINTER,SAVE :: f_dulon(:)
+  TYPE(t_field),POINTER,SAVE :: f_dulat(:)
+ 
+  INTEGER,SAVE  :: n_sponge
+  !$OMP THREADPRIVATE(n_sponge)
+  INTEGER,SAVE  :: begin_sponge
+  !$OMP THREADPRIVATE(begin_sponge)
+  
+CONTAINS
+
+  SUBROUTINE pre_init_sponge
+  USE icosa
+  USE disvert_mod
+  IMPLICIT NONE
+    INTEGER :: l
+
+    tau_sponge = 1.e-5
+    CALL getin("tau_sponge",tau_sponge)
+    PRINT*,'tau_sponge = ',tau_sponge
+
+    iflag_sponge = 0
+    CALL getin("iflag_sponge",iflag_sponge)
+    PRINT*,'iflag_sponge = ',iflag_sponge
+    
+    mode_sponge = 1
+    CALL getin("mode_sponge",mode_sponge)
+    PRINT*,'mode_sponge = ',mode_sponge
+
+    IF (iflag_sponge == 0) THEN
+      PRINT*,'init_sponge: no sponge'
+      RETURN
+    ENDIF
+    
+    IF (.NOT. (iflag_sponge == 1 .OR. iflag_sponge == 2 .OR. iflag_sponge == 3)) THEN
+      PRINT*,'Bad selector for variable iflag_sponge : <',iflag_sponge,'> options are 0,1,2'
+      STOP
+    ELSE
+      IF (.NOT. (mode_sponge == 1 .OR. mode_sponge == 2 .OR. mode_sponge == 3)) THEN
+        PRINT*,'Bad selector for variable mode_sponge : <',mode_sponge,'> options are 0,1,2'
+        STOP
+      ENDIF
+    ENDIF
+
+
+    IF (iflag_sponge == 0) RETURN
+    
+    ALLOCATE(rdamp(llm))
+    ALLOCATE(lambda(llm))
+    !$acc enter data create(rdamp(:),lambda(:)) async
+    
+    IF (iflag_sponge == 1) THEN
+! sponge quenching over the topmost 4 atmospheric layers
+        lambda(:)=0.
+        lambda(llm)=tau_sponge
+        lambda(llm-1)=tau_sponge/2.
+        lambda(llm-2)=tau_sponge/4.
+        lambda(llm-3)=tau_sponge/8.
+        n_sponge=4
+        begin_sponge=llm-3
+    ELSE IF (iflag_sponge == 2 ) THEN
+! sponge quenching over topmost layers down to pressures which are
+! higher than 100 times the topmost layer pressure
+        lambda(:)=tau_sponge*max(presnivs(llm)/presnivs(:)-0.01,0.)
+        DO l=llm,1,-1
+          IF (lambda(l)/=0.) begin_sponge=l
+        ENDDO
+        n_sponge=llm-begin_sponge+1
+    ENDIF
+     
+    IF (iflag_sponge>=1 .AND. (mode_sponge==2 .OR. mode_sponge==3)) THEN
+     CALL pre_init_sponge_zonal_mean_u
+    ENDIF
+    
+  END SUBROUTINE pre_init_sponge
+ 
+ 
+  SUBROUTINE init_sponge
+  USE icosa
+  USE disvert_mod
+  USE omp_para
+  USE mpipara, ONLY: is_mpi_master
+  IMPLICIT NONE
+    INTEGER :: l
+  
+    IF (iflag_sponge == 0) RETURN
+
+    rdamp(:)=itau_dissip*lambda(:)
+
+
+    IF (is_mpi_master) THEN
+      PRINT*,'lambda = '
+      DO l=ll_begin,ll_end
+        PRINT*,l,lambda(l)
+      ENDDO
+      PRINT*,'rdamp = '
+      DO l=ll_begin,ll_end
+        PRINT*,l,rdamp(l)
+      ENDDO
+    ENDIF
+    !$acc update device(rdamp(:), lambda(:)) async
+ 
+   IF (iflag_sponge>=1 .AND. (mode_sponge == 2 .OR. mode_sponge == 3)) THEN
+     CALL init_sponge_zonal_mean_u
+   ENDIF
+   
+  END SUBROUTINE init_sponge
+
+  
+  SUBROUTINE sponge(f_mass, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz)
+  USE icosa
+  USE theta2theta_rhodz_mod
+  USE exner_mod
+  USE geopotential_mod
+  USE trace
+  USE time_mod
+  USE omp_para
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_mass(:)
+    TYPE(t_field),POINTER :: f_due(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+
+    REAL(rstd),POINTER         :: due(:,:)!,due_sponge(:,:)
+    REAL(rstd),POINTER         :: ue(:,:)
+    REAL(rstd),POINTER         :: theta_rhodz(:,:,:)
+    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)!,dtheta_sponge(:,:)
+
+    INTEGER :: ind
+    INTEGER :: l,ij
+
+!$OMP BARRIER
+    
+    CALL trace_start("sponge")
+
+    IF (mode_sponge == 1) THEN
+      DO ind=1,ndomain
+        IF (.NOT. assigned_domain(ind)) CYCLE
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+        ue=f_ue(ind) 
+        due=f_due(ind) 
+        theta_rhodz=f_theta_rhodz(ind)
+        dtheta_rhodz=f_dtheta_rhodz(ind)
+
+       !$acc parallel loop collapse(2) present(ue(:,:), due(:,:), dtheta_rhodz)  async
+        DO l=ll_begin,ll_end
+         !$SIMD
+          DO ij=ij_begin,ij_end
+
+            due(ij+u_right,l) = ue(ij+u_right,l)
+            due(ij+u_lup,l)   = ue(ij+u_lup,l)
+            due(ij+u_ldown,l) = ue(ij+u_ldown,l)
+
+            dtheta_rhodz(ij,l,:) = 0.0
+          ENDDO
+        ENDDO
+       !$acc end parallel loop
+      END DO
+    ELSE IF (mode_sponge == 2 .OR. mode_sponge == 3) THEN
+      CALL sponge_zonal_mean_u(mode_sponge, f_mass, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz) 
+    ELSE
+      PRINT*,'Arnaud Bad selector for variable mode_sponge : <',mode_sponge,             &
+            '> options 2 and 3 not available for the moment!'
+      STOP
+    ENDIF
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      due=f_due(ind) 
+      dtheta_rhodz=f_dtheta_rhodz(ind)
+
+      !$acc parallel loop collapse(2) present(ue(:,:), due(:,:), rdamp(:), dtheta_rhodz)  async
+      DO l=ll_begin,ll_end
+      !$SIMD
+        DO ij=ij_begin,ij_end
+          due(ij+u_right,l)        = -rdamp(l)*due(ij+u_right,l)
+          due(ij+u_lup,l)          = -rdamp(l)*due(ij+u_lup,l)
+          due(ij+u_ldown,l)        = -rdamp(l)*due(ij+u_ldown,l)
+          dtheta_rhodz(ij,l,:)     = -rdamp(l)*dtheta_rhodz(ij,l,:)
+        ENDDO
+      ENDDO
+     !$acc end parallel loop
+    END DO
+
+   CALL trace_end("sponge")
+
+!$OMP BARRIER
+  END SUBROUTINE sponge
+
+  SUBROUTINE pre_init_sponge_zonal_mean_u()
+  USE xios_mod
+  USE grid_param, ONLY : iim_glo
+  USE disvert_mod
+  IMPLICIT NONE
+  INTEGER :: nlon
+  INTEGER :: nlat  
+
+!$OMP MASTER
+    nlon=360.*iim_glo/80.
+    nlat=180.*iim_glo/80.
+
+    CALL xios_set_axis_attr("lev_sponge", n_glo=n_sponge, value=presnivs(begin_sponge:))
+    CALL xios_set_domain_attr("to_sponge", ni_glo=nlon, nj_glo=nlat)
+    CALL xios_set_axis_attr("reduced_sponge", n_glo=nlat)
+    CALL xios_set_field_attr("ulon_sponge_reduced", read_access=.TRUE.)
+    CALL xios_set_field_attr("ulon_from_sponge", read_access=.TRUE.)
+    CALL xios_set_field_attr("ulat_sponge_reduced", read_access=.TRUE.)
+    CALL xios_set_field_attr("ulat_from_sponge", read_access=.TRUE.)
+    CALL xios_set_field_attr("rhodz_sponge_reduced", read_access=.TRUE.)
+    IF (mode_sponge==3) CALL xios_set_field_attr("theta_sponge_reduced", read_access=.TRUE.)
+    IF (mode_sponge==3) CALL xios_set_field_attr("theta_from_sponge", read_access=.TRUE.)
+!$OMP END MASTER
+
+  END SUBROUTINE  pre_init_sponge_zonal_mean_u
+
+
+  SUBROUTINE init_sponge_zonal_mean_u()
+  USE xios_mod
+  IMPLICIT NONE
+    
+!$OMP MASTER
+    CALL xios_get_domain_attr("to_sponge", ni_glo=ni_glo, nj_glo=nj_glo,ni=ni, nj=nj, ibegin=ibegin,jbegin=jbegin)
+!$OMP END MASTER
+    CALL allocate_field(f_hflux, field_u, type_real, llm, name='hflux') 
+    CALL allocate_field(f_ulon, field_t, type_real, llm, name='ulon') 
+    CALL allocate_field(f_ulat, field_t, type_real, llm, name='ulat') 
+    CALL allocate_field(f_ulon_sponge, field_t, type_real, n_sponge, name='ulon_sponge') 
+    CALL allocate_field(f_ulat_sponge, field_t, type_real, n_sponge, name='ulat_sponge') 
+    CALL allocate_field(f_rhodz_sponge, field_t, type_real, n_sponge, name='rhodz_sponge') 
+    CALL allocate_field(f_theta_sponge, field_t, type_real, n_sponge, name='theta_sponge') 
+    CALL allocate_field(f_dulon, field_t, type_real, llm,name='dulon_sponge') 
+    CALL allocate_field(f_dulat, field_t, type_real, llm,name='dulat_sponge') 
+
+  END SUBROUTINE init_sponge_zonal_mean_u
+
+  
+  SUBROUTINE sponge_zonal_mean_u(mode_sponge, f_rhodz, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz)
+  USE icosa
+  USE theta2theta_rhodz_mod
+  USE exner_mod
+  USE geopotential_mod
+  USE trace
+  USE time_mod
+  USE omp_para
+  USE wind_mod
+  USE xios_mod
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)    :: mode_sponge
+    TYPE(t_field),POINTER :: f_rhodz(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_due(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+
+    REAL(rstd),POINTER         :: due(:,:)!,due_sponge(:,:)
+    REAL(rstd),POINTER         :: ue(:,:)
+    REAL(rstd),POINTER         :: rhodz(:,:)
+    REAL(rstd),POINTER         :: hflux(:,:)
+    REAL(rstd),POINTER         :: theta_rhodz(:,:,:)
+    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)!,dtheta_sponge(:,:)
+    REAL(rstd),POINTER         :: ulon(:,:)
+    REAL(rstd),POINTER         :: ulat(:,:)
+    REAL(rstd),POINTER         :: ulon_sponge(:,:)
+    REAL(rstd),POINTER         :: ulat_sponge(:,:)
+    REAL(rstd),POINTER         :: theta_sponge(:,:)
+    REAL(rstd),POINTER         :: rhodz_sponge(:,:)
+    
+    REAL(rstd) :: ulon_zonal_mean(nj_glo,n_sponge)
+    REAL(rstd) :: ulat_zonal_mean(nj_glo,n_sponge)
+    REAL(rstd) :: theta_zonal_mean(nj_glo,n_sponge)
+    REAL(rstd) :: rhodz_zonal_mean(nj_glo,n_sponge)
+    REAL(rstd) :: domain_zonal_mean(ni,nj,n_sponge)
+   
+    INTEGER :: i,j,ind,ij,l,ll 
+
+    DO ind = 1 , ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ue=f_ue(ind)
+      rhodz=f_rhodz(ind)
+      hflux = f_hflux(ind)
+         
+      DO l = ll_begin, ll_end
+         !DIR$ IVDEP
+         DO ij=ij_begin_ext,ij_end_ext         
+            hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*ue(ij+u_right,l)
+            hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*ue(ij+u_lup,l)
+            hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*ue(ij+u_ldown,l)
+         ENDDO
+       ENDDO
+    ENDDO
+
+              
+    CALL transfert_request(f_hflux,req_e1_vect)
+    CALL un2ulonlat(f_hflux, f_ulon, f_ulat)
+
+    DO ind = 1 , ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      
+      ulon=f_ulon(ind)
+      ulat=f_ulat(ind)
+      ulon_sponge = f_ulon_sponge(ind)
+      ulat_sponge = f_ulat_sponge(ind)
+      rhodz_sponge = f_rhodz_sponge(ind)
+      theta_sponge = f_theta_sponge(ind)
+      theta_rhodz = f_theta_rhodz(ind)
+              
+      DO l = 1, llm
+         IF (l<begin_sponge) CYCLE
+         ll=l-begin_sponge+1
+         !DIR$ IVDEP
+         DO ij=ij_begin_ext,ij_end_ext         
+            ulon_sponge(ij,ll)= ulon(ij,l) 
+            ulat_sponge(ij,ll)= ulat(ij,l) 
+            rhodz_sponge(ij,ll) = rhodz(ij,l)
+            IF (mode_sponge==3) theta_sponge(ij,ll) = theta_rhodz(ij,l,1)
+         ENDDO
+       ENDDO
+    ENDDO
+
+!$OMP BARRIER
+    
+    CALL xios_write_field("ulon_sponge",f_ulon_sponge)
+    CALL xios_write_field("ulat_sponge",f_ulat_sponge)
+    CALL xios_write_field("rhodz_sponge",f_rhodz_sponge)
+!$OMP MASTER
+    CALL xios_recv_field("ulon_sponge_reduced",ulon_zonal_mean)
+    CALL xios_recv_field("ulat_sponge_reduced",ulat_zonal_mean)
+    CALL xios_recv_field("rhodz_sponge_reduced",rhodz_zonal_mean)
+!$OMP END MASTER
+
+    IF (mode_sponge==3) THEN
+      CALL xios_write_field("theta_sponge",f_theta_sponge)
+!$OMP MASTER
+      CALL xios_recv_field("theta_sponge_reduced",theta_zonal_mean)
+!$OMP END MASTER
+    ENDIF
+    
+!$OMP MASTER
+    DO j=1,nj
+      DO i=1,ni
+        domain_zonal_mean(i,j,:)=ulon_zonal_mean(jbegin+j,:)/rhodz_zonal_mean(jbegin+j,:)
+      ENDDO
+    ENDDO
+    
+    CALL xios_send_field("ulon_sponge_mean",domain_zonal_mean)
+!$OMP END MASTER
+    CALL xios_read_field("ulon_from_sponge",f_ulon_sponge)
+  
+!$OMP MASTER
+    DO j=1,nj
+      DO i=1,ni
+        domain_zonal_mean(i,j,:)=ulat_zonal_mean(jbegin+j,:)/rhodz_zonal_mean(jbegin+j,:)
+      ENDDO
+    ENDDO
+    
+    CALL xios_send_field("ulat_sponge_mean",domain_zonal_mean)
+!$OMP END MASTER
+    CALL xios_read_field("ulat_from_sponge",f_ulat_sponge)
+    
+    IF (mode_sponge==3) THEN
+!$OMP MASTER
+      DO j=1,nj
+        DO i=1,ni
+          domain_zonal_mean(i,j,:)=theta_zonal_mean(jbegin+j,:)
+        ENDDO
+      ENDDO
+      CALL xios_send_field("theta_sponge_mean",domain_zonal_mean)
+!$OMP END MASTER
+      CALL xios_read_field("theta_from_sponge",f_theta_sponge)
+    ENDIF
+    
+
+    DO ind = 1 , ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      
+      ulon=f_ulon(ind)
+      ulat=f_ulat(ind)
+      ulon_sponge = f_ulon_sponge(ind)
+      ulat_sponge = f_ulat_sponge(ind)
+      theta_sponge = f_theta_sponge(ind)
+      dtheta_rhodz = f_dtheta_rhodz(ind)    
+      theta_rhodz = f_theta_rhodz(ind)
+      DO l = ll_begin, ll_end
+        IF (l<begin_sponge) CYCLE
+        ll=l-begin_sponge+1
+        IF (mode_sponge==2) THEN
+          DO ij=ij_begin_ext,ij_end_ext         
+            ulon(ij,l) = ulon_sponge(ij,ll) 
+            ulat(ij,l) = ulat_sponge(ij,ll) 
+            dtheta_rhodz(ij,l,1)=theta_rhodz(ij,ll,1) 
+          ENDDO
+        ELSE  IF (mode_sponge==3) THEN
+          DO ij=ij_begin_ext,ij_end_ext         
+            ulon(ij,l) = ulon_sponge(ij,ll) 
+            ulat(ij,l) = ulat_sponge(ij,ll) 
+            dtheta_rhodz(ij,l,1) = theta_sponge(ij,ll)
+          ENDDO
+        ENDIF
+      ENDDO
+    ENDDO
+    
+    CALL transfert_request(f_ulon,req_i0)
+    CALL transfert_request(f_ulat,req_i0)
+    IF (mode_sponge==3) THEN
+      CALL transfert_request(f_dtheta_rhodz,req_i0)
+    ENDIF 
+    CALL ulonlat2un(f_ulon,f_ulat,f_due)
+!$OMP BARRIER
+   
+    DO ind = 1 , ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      due=f_due(ind)
+      ue=f_ue(ind)
+      dtheta_rhodz=f_dtheta_rhodz(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+ 
+       DO l = ll_begin, ll_end
+         IF (l<begin_sponge) THEN
+           !DIR$ IVDEP
+           DO ij=ij_begin,ij_end         
+              due(ij+u_right,l) = 0.
+              due(ij+u_lup,l)   = 0.
+              due(ij+u_ldown,l) = 0.
+              dtheta_rhodz(ij,l,:) = 0.
+           ENDDO         
+         ELSE
+           !DIR$ IVDEP
+           DO ij=ij_begin,ij_end         
+             due(ij+u_right,l) = ue(ij+u_right,l) - due(ij+u_right,l)
+             due(ij+u_lup,l)   = ue(ij+u_lup,l)   - due(ij+u_lup,l)
+             due(ij+u_ldown,l) = ue(ij+u_ldown,l) - due(ij+u_ldown,l)
+             dtheta_rhodz(ij,l,1) = theta_rhodz(ij,l,1) - dtheta_rhodz(ij,l,1)
+          ENDDO
+        ENDIF
+      ENDDO
+    ENDDO  
+   
+    CALL transfert_request(f_due,req_e1_vect)
+!ym to check results
+!    CALL un2ulonlat(f_due, f_dulon, f_dulat)
+!    CALL xios_write_field("dulon_sponge",f_dulon)
+!    CALL xios_write_field("dulat_sponge",f_dulat)
+!    CALL xios_write_field("dtheta_sponge",f_dtheta_rhodz)
+    
+  END SUBROUTINE sponge_zonal_mean_u
+    
+
+END MODULE sponge_mod
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/icosa_init.f90
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/icosa_init.f90	(revision 6362)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/icosa_init.f90	(revision 6362)
@@ -0,0 +1,144 @@
+MODULE icosa_init_mod
+  
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: icosa_init
+  
+CONTAINS
+  
+  SUBROUTINE icosa_init
+    USE profiling_mod,        ONLY : init_profiling
+    USE mpipara,              ONLY : is_mpi_master, is_mpi_root, mpi_rank, init_mpipara, finalize_mpipara
+    USE signature_mod,        ONLY : print_signature
+    USE trace,                ONLY : trace_off, init_trace
+    USE xios_mod,             ONLY : xios_init
+    USE geometry,             ONLY : geom
+    USE earth_const,          ONLY : init_earth_const
+    USE grid_param    ,       ONLY : grid_type, grid_ico
+    USE init_grid_param_mod,  ONLY : init_grid_param
+    USE omp_para,             ONLY : init_omp_para, switch_omp_no_distrib_level
+    USE metric,               ONLY : init_metric
+    USE domain_mod,           ONLY : compute_domain
+    USE transfert_mod,        ONLY : init_transfert
+    USE write_field_mod,      ONLY : init_writefield, writefield, close_files
+    USE trace,                ONLY : init_trace
+    USE read_metric_mod,      ONLY : read_metric
+    USE sponge_mod,           ONLY : pre_init_sponge
+    USE geometry,             ONLY : compute_geometry
+    USE disvert_mod,          ONLY : init_disvert, write_apbp
+    USE vertical_interp_mod,  ONLY : init_vertical_interp
+    USE tracer_icosa_mod,     ONLY : init_tracer
+    USE time_mod,             ONLY : init_time, close_time_counter
+    USE restart_mod,          ONLY : init_restart
+    USE etat0_mod,            ONLY : init_etat0
+    USE output_field_mod,     ONLY : output_field_init, output_field_finalize
+    USE timeloop_gcm_mod,     ONLY : init_timeloop, timeloop
+    USE physics_mod,          ONLY : init_physics, zero_du_phys
+    USE diagflux_mod,         ONLY : init_diagflux
+    USE guided_mod,           ONLY : pre_init_guided
+    
+    CALL init_profiling
+    CALL init_mpipara
+    IF(is_mpi_master) CALL print_signature
+    CALL trace_off
+    CALL xios_init
+    CALL init_earth_const 
+    CALL init_grid_param(is_mpi_master, mpi_rank) ! sets max_primal_deg, etc.
+    CALL init_omp_para(is_mpi_master)
+    IF(grid_type == grid_ico) CALL init_metric
+    CALL compute_domain ! allocates domain_glo(ndomain_glo), domain(ndomain)
+    CALL init_transfert
+    
+    CALL init_writefield
+    CALL init_trace
+    
+    !$OMP PARALLEL  
+    CALL switch_omp_no_distrib_level
+
+    IF(grid_type == grid_ico) CALL read_metric
+    CALL compute_geometry
+    
+    IF(grid_type == grid_ico) CALL check_total_area
+    
+    CALL init_disvert 
+    CALL init_vertical_interp
+    CALL init_tracer
+    
+    CALL WriteField("Ai",geom%Ai)
+    
+    IF (is_mpi_root) CALL write_apbp
+    CALL init_time
+    CALL init_etat0
+    CALL init_restart
+    CALL pre_init_sponge
+    CALL pre_init_guided ! must be done before closing xios context, i.e. before output_field_init
+    CALL output_field_init
+    CALL init_timeloop
+    
+    !$OMP END PARALLEL
+    
+    CALL init_physics
+    
+    !$OMP PARALLEL
+    
+    CALL init_diagflux
+    CALL zero_du_phys
+    
+    CALL timeloop
+    
+    CALL switch_omp_no_distrib_level
+    !$OMP END PARALLEL
+    
+    CALL output_field_finalize
+    CALL close_files
+    CALL close_time_counter
+    CALL finalize_mpipara
+    
+  END SUBROUTINE icosa_init
+  
+  SUBROUTINE check_total_area
+    USE prec,       ONLY : rstd
+    USE domain_mod, ONLY : ndomain, domain
+    USE geometry,   ONLY : swap_geometry, Ai, ne
+    USE mpipara,    ONLY : is_mpi_root
+    USE dimensions
+    USE field_mod
+    USE domain_param
+    USE omp_para
+    USE icosa
+    
+    TYPE(t_field), POINTER, SAVE :: sum_ne(:)
+    REAL(rstd), POINTER          :: pt_sum_ne(:)
+    INTEGER :: ind,i,j,k,n 
+    REAL(rstd) :: tot_sum=0
+    
+    CALL allocate_field(sum_ne,field_T,type_real)
+    
+    !$OMP BARRIER
+    !$OMP MASTER
+    DO ind=1,ndomain
+       
+       pt_sum_ne=sum_ne(ind)
+       
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             n=(j-1)*iim+i
+             pt_sum_ne(n)=0
+             DO k=1,6
+                pt_sum_ne(n)=pt_sum_ne(n)+ne(n,k)
+             ENDDO
+             IF (domain(ind)%own(i,j)) tot_sum=tot_sum+Ai(n)
+          ENDDO
+       ENDDO
+    ENDDO
+    
+    IF (is_mpi_root) PRINT *," Diff surf",1-tot_sum/(4*Pi*radius*radius)
+    !$OMP END MASTER
+    
+  END SUBROUTINE check_total_area
+  
+END MODULE icosa_init_mod
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/time/timeloop_gcm.F90
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/time/timeloop_gcm.F90	(revision 6362)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/DYNAMICO/trunk/time/timeloop_gcm.F90	(revision 6362)
@@ -0,0 +1,535 @@
+MODULE timeloop_gcm_mod
+  USE profiling_mod
+  USE icosa
+  USE disvert_mod
+  USE trace
+  USE omp_para
+  USE euler_scheme_mod
+  USE explicit_scheme_mod
+  USE hevi_scheme_mod
+  IMPLICIT NONE
+  PRIVATE
+  
+  INTEGER, PARAMETER :: sync_it=10
+  TYPE(t_message),SAVE :: req_ps0, req_rhodz0, req_mass0, req_theta_rhodz0, req_u0, req_q0, req_W0, req_geopot0
+  LOGICAL, SAVE :: positive_theta
+  INTEGER :: itau_prof, id_timeloop, id_dyn, id_phys, id_dissip, id_adv, id_diags
+  PUBLIC :: init_timeloop, timeloop
+
+CONTAINS
+  
+  SUBROUTINE init_timeloop
+    USE dissip_gcm_mod
+    USE observable_mod
+    USE caldyn_mod
+    USE etat0_mod
+    USE guided_mod
+    USE advect_tracer_mod
+    USE check_conserve_mod
+    USE output_field_mod
+    USE theta2theta_rhodz_mod
+    USE sponge_mod
+
+    CHARACTER(len=255) :: def
+
+    CALL register_id('timeloop',id_timeloop)
+    CALL register_id('dyn',id_dyn)
+    CALL register_id('dissip',id_dissip)
+    CALL register_id('phys',id_phys)
+    CALL register_id('adv',id_adv)
+    CALL register_id('diags',id_diags)
+
+    CALL init_caldyn
+    
+!    IF (xios_output) itau_out=1
+    IF (.NOT. enable_io) itau_out=HUGE(itau_out)
+
+    itau_prof=1000
+    CALL getin('itau_profiling',itau_prof)
+
+    positive_theta=.FALSE.
+    CALL getin('positive_theta',positive_theta)
+    IF(positive_theta .AND. nqtot<1) THEN
+       PRINT *, 'nqtot must be >0 if positive_theta is .TRUE.'
+       STOP
+    END IF
+
+    def='ARK2.3'
+    CALL getin('time_scheme',def)
+    SELECT CASE (TRIM(def))
+    CASE('euler')
+       scheme_family=explicit
+       scheme=euler
+       nb_stage=1
+    CASE ('runge_kutta')
+       scheme_family=explicit
+       scheme=rk4
+       nb_stage=4
+    CASE ('RK2.5')
+       scheme_family=explicit
+       scheme=rk25
+       nb_stage=5
+    CASE ('leapfrog_matsuno')
+       scheme_family=explicit
+       scheme=mlf
+       matsuno_period=5
+       CALL getin('matsuno_period',matsuno_period)
+       nb_stage=matsuno_period+1
+    CASE('ARK2.3')
+       scheme_family=hevi
+       scheme=ark23
+       nb_stage=3
+       CALL set_coefs_ark23(dt)
+    CASE('ARK3.3')
+       scheme_family=hevi
+       scheme=ark33
+       nb_stage=3
+       CALL set_coefs_ark33(dt)
+    CASE ('none')
+       nb_stage=0
+    CASE default
+       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             &
+            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>,<RK2.5>,<ARK2.3>'
+       STOP
+    END SELECT
+    
+    IF (scheme_family /= hevi) THEN
+       CALL abort_acc("scheme_family /= hevi")
+    END IF
+
+    ! Time-independant orography
+    CALL allocate_field(f_phis,field_t,type_real,name='phis')
+    ! Model state at current time step
+    CALL allocate_field(f_ps,field_t,type_real, name='ps')
+    CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
+    CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
+    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,nqdyn,name='theta_rhodz')
+    CALL allocate_field(f_u,field_u,type_real,llm,name='u')
+    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')
+    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W') ! used only if .not. hydrostatic
+    CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q')
+    ! Mass fluxes
+    CALL allocate_field(f_hflux,field_u,type_real,llm, ondevice=.TRUE.)    ! instantaneous mass fluxes
+    CALL allocate_field(f_hfluxt,field_u,type_real,llm,ondevice=.TRUE.)   ! mass "fluxes" accumulated in time
+    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes
+    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt',ondevice=.TRUE.)
+    
+    SELECT CASE(scheme_family)
+    CASE(explicit)
+       ! Trends
+       CALL allocate_field(f_dps,field_t,type_real,name='dps')
+       CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
+       CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,nqdyn,name='dtheta_rhodz')
+       CALL allocate_field(f_du,field_u,type_real,llm,name='du')
+       ! Model state at previous time step (RK/MLF)
+       CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
+       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
+       CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,nqdyn,name='theta_rhodzm1')
+       CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')
+    CASE(hevi)
+       ! Trends
+       CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_dW_slow, field_t,type_real,llm+1,name='dW_slow')
+       CALL allocate_fields(nb_stage,f_dW_fast, field_t,type_real,llm+1,name='dW_fast')
+       CALL allocate_fields(nb_stage,f_dPhi_slow, field_t,type_real,llm+1,name='dPhi_slow')
+       CALL allocate_fields(nb_stage,f_dPhi_fast, field_t,type_real,llm+1,name='dPhi_fast')
+       f_dps => f_dps_slow(:,1)
+       f_du => f_du_slow(:,1)
+       f_dtheta_rhodz => f_dtheta_rhodz_slow(:,1)
+    END SELECT
+
+    SELECT CASE(scheme)
+    CASE(mlf)
+       ! Model state 2 time steps ago (MLF)
+       CALL allocate_field(f_psm2,field_t,type_real)
+       CALL allocate_field(f_massm2,field_t,type_real,llm)
+       CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm,nqdyn)
+       CALL allocate_field(f_um2,field_u,type_real,llm)
+    END SELECT
+
+    CALL init_theta2theta_rhodz
+    CALL init_dissip
+    CALL init_sponge
+    CALL init_observable
+    CALL init_guided
+    CALL init_advect_tracer
+    CALL init_check_conserve
+
+    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
+
+    CALL transfert_request(f_phis,req_i0) 
+    CALL transfert_request(f_phis,req_i1) 
+
+    CALL init_message(f_ps,req_i0,req_ps0)
+    CALL init_message(f_mass,req_i0,req_mass0)
+    CALL init_message(f_rhodz,req_i0,req_rhodz0)
+    CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
+    CALL init_message(f_u,req_e0_vect,req_u0)
+    CALL init_message(f_q,req_i0,req_q0)
+    CALL init_message(f_geopot,req_i0,req_geopot0)
+    CALL init_message(f_W,req_i0,req_W0)
+
+  END SUBROUTINE init_timeloop
+  
+  SUBROUTINE timeloop
+    USE abort_mod
+    USE dissip_gcm_mod
+    USE sponge_mod
+    USE observable_mod
+    USE etat0_mod
+    USE guided_mod
+    USE caldyn_mod
+    USE advect_tracer_mod
+    USE diagflux_mod
+    USE physics_mod
+    USE mpipara
+    USE transfert_mod
+    USE check_conserve_mod
+    USE xios_mod
+    USE output_field_mod
+    USE write_etat0_mod
+    USE restart_mod
+    USE checksum_mod
+    USE explicit_scheme_mod
+    USE hevi_scheme_mod
+    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:)
+
+    REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out
+    INTEGER :: ind, it,l
+    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time
+    LOGICAL, PARAMETER :: check_rhodz=.FALSE.
+    INTEGER(kind=8) :: start_clock, stop_clock, rate_clock
+    INTEGER :: itau_sync   ! best iteration for synchronisation and ensure 1+1=2
+    INTEGER :: i
+    
+    LOGICAL,SAVE :: first_physic=.TRUE.
+    !$OMP THREADPRIVATE(first_physic)    
+
+    itau_sync=1
+    DO i=2,3*sync_it
+      IF (MOD(86400,INT(i*dt))==0 .AND. ABS((sync_it-itau_sync)*1./sync_it )/sync_it < (sync_it-itau_sync)*1./sync_it) itau_sync=i
+    ENDDO
+    CALL getin("itau_sync",itau_sync)
+    IF (is_master) PRINT*,"Synchronize frontier every itau_sync =",itau_sync
+      
+    CALL switch_omp_distrib_level
+    CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces
+
+    !$OMP BARRIER
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
+       IF(caldyn_eta==eta_mass) THEN
+          CALL compute_rhodz(.TRUE., ps, rhodz, ondevice=.FALSE.) ! save rhodz for transport scheme before dynamics update ps
+       ELSE
+          DO l=ll_begin,ll_end
+             rhodz(:,l)=mass(:,l)
+          ENDDO
+       END IF
+    END DO
+    !$OMP BARRIER
+
+    fluxt_zero=.TRUE.
+
+    IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
+    IF(diagflux_on) THEN
+       adv_over_out = itau_adv*(1./itau_out)
+    ELSE
+       adv_over_out = 0.
+    END IF
+
+    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,itau0)  
+
+    Call trace_on
+
+    IF (xios_output) THEN ! we must call update_calendar before any XIOS output
+      IF (is_omp_master) CALL xios_update_calendar(1)
+    END IF
+    
+!    IF (is_write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
+    IF (is_write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W)
+   
+    CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
+
+    !$OMP MASTER
+    CALL SYSTEM_CLOCK(start_clock, rate_clock)
+    !$OMP END MASTER   
+    call update_device_field(f_ps)
+    call update_device_field(f_mass)
+    CALL update_device_field(f_theta_rhodz)
+    CALL update_device_field(f_u)
+    CALL update_device_field(f_q)
+    CALL update_device_field(f_geopot)
+    CALL update_device_field(f_wflux)
+    CALL update_device_field(f_rhodz)
+    call reset_profiling()  
+
+
+    DO it=itau0+1,itau0+itaumax
+
+       CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)
+
+       CALL enter_profile(id_timeloop)
+
+       IF (xios_output) THEN
+
+          IF(it>itau0+1 .AND. is_omp_master) CALL xios_update_calendar(it-itau0)
+       ELSE
+          CALL update_time_counter(dt*it)
+       END IF
+
+       IF (it==itau0+1 .OR. MOD(it-1,itau_sync)==0) THEN
+          CALL send_message(f_ps,req_ps0)
+          CALL wait_message(req_ps0)
+          CALL send_message(f_mass,req_mass0)
+          CALL wait_message(req_mass0)
+          CALL send_message(f_rhodz,req_rhodz0)
+          CALL wait_message(req_rhodz0)
+          CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
+          CALL wait_message(req_theta_rhodz0)
+          CALL send_message(f_u,req_u0)
+          CALL wait_message(req_u0)
+          CALL send_message(f_q,req_q0) 
+          CALL wait_message(req_q0)
+          IF(.NOT. hydrostatic) THEN
+             CALL send_message(f_geopot,req_geopot0)
+             CALL wait_message(req_geopot0)
+             CALL send_message(f_W,req_W0)
+             CALL wait_message(req_W0)
+          END IF
+       ENDIF
+
+!       CALL checksum(f_ps)
+!       CALL checksum(f_theta_rhodz)
+!       CALL checksum(f_u)
+!       CALL checksum(f_q)
+!       CALL checksum(f_mass)
+!       CALL checksum(f_geopot)
+!       CALL checksum(f_rhodz)
+!       CALL checksum(f_wflux)
+
+       CALL guided(it*dt, f_ps, f_rhodz, f_theta_rhodz,f_u,f_q)
+       
+       CALL enter_profile(id_dyn)
+       SELECT CASE(scheme_family)
+       CASE(explicit)
+          CALL abort_acc("explicit_scheme")
+          CALL explicit_scheme(it, fluxt_zero)
+       CASE(hevi)
+          CALL HEVI_scheme(it, fluxt_zero)
+       END SELECT
+       CALL exit_profile(id_dyn)
+
+!       CALL checksum(f_ps)
+!       CALL checksum(f_theta_rhodz)
+!       CALL checksum(f_u)
+!       CALL checksum(f_q)
+!       CALL checksum(f_mass)
+!       CALL checksum(f_geopot)
+!       CALL checksum(f_rhodz)
+!       CALL checksum(f_wflux)
+
+       CALL enter_profile(id_dissip)
+       IF (MOD(it,itau_dissip)==0) THEN
+          
+          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+          
+          CALL enter_profile(id_diags)
+          CALL check_conserve_detailed(it, AAM_dyn, &
+               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL exit_profile(id_diags)
+
+          CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du)
+          
+          CALL euler_scheme(.FALSE.)  ! update only u, theta
+          IF (iflag_sponge > 0) THEN
+             CALL sponge(f_mass,f_u,f_du,f_theta_rhodz,f_dtheta_rhodz)
+             CALL euler_scheme(.FALSE.)  ! update only u, theta
+          END IF
+
+          CALL update_border(it*dt,f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+
+          CALL enter_profile(id_diags)
+          CALL check_conserve_detailed(it, AAM_dissip, &
+               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL exit_profile(id_diags)
+       END IF
+       CALL exit_profile(id_dissip)
+!       CALL checksum(f_ps)
+!       CALL checksum(f_theta_rhodz)
+!       CALL checksum(f_u)
+!       CALL checksum(f_q)
+!       CALL checksum(f_hfluxt)
+!       CALL checksum(f_wfluxt)
+!       CALL checksum(f_u)
+!       CALL checksum(f_rhodz)
+
+       CALL enter_profile(id_adv)
+       IF(MOD(it,itau_adv)==0) THEN
+          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
+               adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt)  ! accumulate mass and fluxes if diagflux_on
+          fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
+          ! At this point advect_tracer has obtained the halos of u and rhodz,
+          ! needed for correct computation of kinetic energy
+          IF(diagflux_on) CALL abort_acc("diagflux_on")
+          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta,f_buf_i, f_hfluxt)
+
+          DO ind=1,ndomain
+            IF (.NOT. assigned_domain(ind)) CYCLE
+            CALL swap_dimensions(ind)
+            CALL swap_geometry(ind)
+            rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
+            IF(caldyn_eta==eta_mass) THEN
+              CALL compute_rhodz(.TRUE., ps, rhodz, ondevice=.TRUE.) ! save rhodz for transport scheme before dynamics update ps
+            ELSE
+              DO l=ll_begin,ll_end
+                 rhodz(:,l)=mass(:,l)
+              ENDDO
+            END IF
+          END DO
+
+          IF(positive_theta) CALL abort_acc("positive_theta")
+          IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
+          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+       END IF
+       CALL exit_profile(id_adv)
+!       CALL checksum(f_ps)
+!       CALL checksum(f_theta_rhodz)
+!       CALL checksum(f_u)
+!       CALL checksum(f_q)
+
+       CALL enter_profile(id_diags)
+!       IF (MOD(it,itau_physics)==0) THEN
+          CALL check_conserve_detailed(it, AAM_dyn, &
+            f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL enter_profile(id_phys)
+          CALL physics(it,f_phis, f_geopot, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
+          CALL exit_profile(id_phys)
+          CALL check_conserve_detailed(it, AAM_phys, &
+               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          !$OMP MASTER
+          IF (first_physic) CALL SYSTEM_CLOCK(start_clock)
+          !$OMP END MASTER   
+          first_physic=.FALSE.
+!       END IF
+
+       IF (MOD(it,itau_check_conserv)==0) THEN
+          CALL update_host_field(f_ps)
+          CALL update_host_field(f_theta_rhodz)
+          CALL update_host_field(f_u)
+          CALL update_host_field(f_dps)
+          CALL update_host_field(f_q)
+          CALL check_conserve_detailed(it, AAM_dyn, &
+               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,it) 
+       ENDIF
+
+       IF (mod(it,itau_out)==0 ) THEN
+ 
+          CALL update_border(it*dt,f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+          
+          CALL transfert_request(f_u,req_e1_vect)
+          CALL update_host_field(f_ps)              
+          CALL update_host_field(f_mass)
+          CALL update_host_field(f_theta_rhodz)
+          CALL update_host_field(f_geopot)
+          CALL update_host_field(f_u)
+          CALL update_host_field(f_q)
+          CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
+       ENDIF
+       CALL exit_profile(id_diags)
+
+       CALL exit_profile(id_timeloop)
+    END DO
+    
+    CALL update_host_field(f_ps)
+    CALL update_host_field(f_theta_rhodz)
+    CALL update_host_field(f_u)
+    CALL update_host_field(f_q)
+    CALL update_host_field(f_geopot)
+
+!    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q) 
+    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W) 
+
+    CALL update_host_field(f_dps)    
+    CALL check_conserve_detailed(it, AAM_dyn, &
+         f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,it)  
+    
+    !$OMP MASTER
+    CALL SYSTEM_CLOCK(stop_clock)
+    
+    IF (mpi_rank==0) THEN 
+       PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock 
+    ENDIF
+    !$OMP END MASTER 
+    
+    ! CONTAINS
+  END SUBROUTINE timeloop
+
+  SUBROUTINE print_iteration(it,itau0,itaumax,start_clock,rate_clock)
+    INTEGER :: it, itau0, itaumax, throughput
+    INTEGER(kind=8) :: start_clock, stop_clock, rate_clock
+    REAL :: per_step,total, elapsed
+    
+    IF(is_master) THEN
+       WRITE(*,'(A,I7,A,F14.1)') "It No :",it,"   t :",dt*it
+       IF(MOD(it,10)==0) THEN
+          CALL SYSTEM_CLOCK(stop_clock)
+          elapsed = (stop_clock-start_clock)*1./rate_clock
+          per_step = elapsed/(it-itau0)
+          throughput = INT(dt/per_step)
+          total = per_step*itaumax
+          WRITE(*,'(A,I5,A,F6.2,A,I6)') 'Time spent (s):',INT(elapsed), &
+               '  -- ms/step : ', 1000*per_step, &
+               '  -- Throughput :', throughput
+          WRITE(*,'(A,I5,A,I5)') 'Whole job (min) :', INT(total/60.), &
+               '  -- Completion in (min) : ', INT((total-elapsed)/60.)
+       END IF
+    END IF
+    IF(MOD(it,itau_prof)==0) CALL print_profile
+
+  END SUBROUTINE print_iteration
+
+  SUBROUTINE copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
+    INTEGER :: ind, iq
+    DO ind=1, ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       rhodz=f_rhodz(ind)
+       q=f_q(ind)
+       DO iq=1, nqdyn
+          q(:,:,iq)  = theta_rhodz(:,:,iq)/rhodz(:,:)
+       END DO
+    END DO
+  END SUBROUTINE copy_theta_to_q
+
+  SUBROUTINE copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
+    INTEGER :: ind, iq
+    DO ind=1, ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       rhodz=f_rhodz(ind)
+       q=f_q(ind)
+       DO iq=1,nqdyn
+          theta_rhodz(:,:,iq) = rhodz(:,:)*q(:,:,iq)
+       END DO
+    END DO
+  END SUBROUTINE copy_q_to_theta
+
+END MODULE timeloop_gcm_mod
Index: /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/compile_icolmdzor.sh
===================================================================
--- /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/compile_icolmdzor.sh	(revision 6361)
+++ /CONFIG/UNIFORM/v7/ICOLMDZOR_v7/compile_icolmdzor.sh	(revision 6362)
@@ -290,5 +290,15 @@
     echo; echo "NOW COMPILE DYNAMICO "
     echo >> $outfile ; echo " NOW COMPILE DYNAMICO"   >> $outfile 
-    
+
+	if [ -d $modipsl/config/ICOLMDZOR_v7/SOURCES/DYNAMICO ] ; then 
+         cd $modipsl/config/ICOLMDZOR_v7/SOURCES/DYNAMICO 
+         for dir in $( find . -mindepth 1 -maxdepth 1 -type d ); do 
+             for file in $( find $dir -type f ! -path '*.svn*' ); do 
+                 echo cp $file $modipsl/modeles/DYNAMICO/$file >> $outfile  
+                 cp $file $modipsl/modeles/DYNAMICO/$file  >> $outfile 2>&1 
+             done 
+         done 
+         cd $modipsl/modeles/DYNAMICO 
+        fi 
 
     echo ./make_icosa -$optmode -with_fcm1 -parallel $parallel -external_ioipsl -with_xios -arch $fcm_arch -arch_path $arch_path -job 8 $full_dyna    >> $outfile 
Index: /CONFIG/UNIFORM/v7/IPSLCM7/EXPERIMENTS/ICOLMDZOR/clim_noleap/COMP/dynamico.card
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/EXPERIMENTS/ICOLMDZOR/clim_noleap/COMP/dynamico.card	(revision 6361)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/EXPERIMENTS/ICOLMDZOR/clim_noleap/COMP/dynamico.card	(revision 6362)
@@ -30,4 +30,5 @@
    	(${SUBMIT_DIR}/PARAM/run_dynamico.def_nbp${RESOL_NBP}, run_dynamico.def),\
    	(${SUBMIT_DIR}/PARAM/iodef.xml, .),\
+	(${SUBMIT_DIR}/PARAM/sponge_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_input_dynamico.xml, .),\
Index: /CONFIG/UNIFORM/v7/IPSLCM7/EXPERIMENTS/IPSLCM/piControl_TEST/COMP/dynamico.card
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/EXPERIMENTS/IPSLCM/piControl_TEST/COMP/dynamico.card	(revision 6361)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/EXPERIMENTS/IPSLCM/piControl_TEST/COMP/dynamico.card	(revision 6362)
@@ -30,4 +30,5 @@
    	(${SUBMIT_DIR}/PARAM/run_dynamico.def_nbp${RESOL_NBP}, run_dynamico.def),\
    	(${SUBMIT_DIR}/PARAM/iodef.xml, .),\
+	(${SUBMIT_DIR}/PARAM/sponge_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_dynamico.xml, .),\
    	(${MODIPSL}/modeles/DYNAMICO/xml/context_input_dynamico.xml, .),\
Index: /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/DRIVER/dynamico.driver
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/DRIVER/dynamico.driver	(revision 6361)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/DRIVER/dynamico.driver	(revision 6362)
@@ -86,10 +86,14 @@
     fi 
 
-
     # Add include of context xml files for DYNAMICO in iodef.xml
     # In iodef.xml add on the next line after "COMPONENT CONTEXT"
     echo '<context id="icosagcm" src="./context_dynamico.xml"/>' >> add.tmp
 
-    # Add inclusion of file context_input_dynamico.xml if this file exists
+   # Add inclusion of file sponge_dynamico.xml if this file exists
+   if [ -f sponge_dynamico.xml ] ; then
+	echo '<context id="icosagcm"        src="./sponge_dynamico.xml" />' >> add.tmp
+   fi
+
+   # Add inclusion of file context_input_dynamico.xml if this file exists
     if [ -f context_input_dynamico.xml ] ; then
         echo '<context id="icosagcm_input" src="./context_input_dynamico.xml"/>' >> add.tmp
Index: /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/PARAM/sponge_dynamico.xml
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/PARAM/sponge_dynamico.xml	(revision 6362)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/PARAM/sponge_dynamico.xml	(revision 6362)
@@ -0,0 +1,111 @@
+<!-- =========================================================================================================== -->
+<!-- nudging_dynamico.xml                                                                                           -->
+<!-- Definition of sponge layer zonal mean                                          -->
+<!-- =========================================================================================================== -->
+<context id="icosagcm" >
+
+
+   <axis_definition>
+     
+     <axis id="reduced_sponge" > <!-- n_glo is set from model -->
+       <reduce_domain  operation="average" direction="iDir"/>
+     </axis>
+
+     <axis id="lev_sponge" axis_type="Z" /> <!-- n_glo is set from model -->
+     
+    </axis_definition>
+    
+    <grid_definition>
+
+      <grid id="sponge_i">
+        <domain domain_ref="i" />
+        <axis axis_ref="lev_sponge"/>
+      </grid>
+   
+      <grid id="to_sponge">
+        <domain id="to_sponge" type="rectilinear" > <!-- ni_glo and nj_glo are set from model -->
+          <generate_rectilinear_domain />
+          <interpolate_domain order="1"/>
+        </domain> 
+        <axis axis_ref="lev_sponge"/>
+      </grid>
+         
+      <grid id="reduced_sponge">
+        <axis axis_ref="reduced_sponge"/>
+        <axis axis_ref="lev_sponge"/>
+      </grid>
+      
+      <grid id="from_sponge">
+        <domain domain_ref="i">
+          <interpolate_domain order="1"/>
+        </domain> 
+        <axis axis_ref="lev_sponge"/>
+      </grid>
+      
+    </grid_definition>
+    
+    <field_definition>
+      <field id="ulon_sponge"                                        grid_ref="sponge_i"/>
+      <field id="ulon_sponge_interp" field_ref="ulon_sponge"         grid_ref="to_sponge" />
+      <field id="ulon_sponge_reduced" field_ref="ulon_sponge_interp" grid_ref="reduced_sponge" read_access="false"/>
+      <field id="ulon_sponge_mean"                                   grid_ref="to_sponge" />
+      <field id="ulon_from_sponge" field_ref="ulon_sponge_mean"      grid_ref="from_sponge" read_access="false"/>
+
+      <field id="ulat_sponge"                                        grid_ref="sponge_i"/>
+      <field id="ulat_sponge_interp" field_ref="ulat_sponge"         grid_ref="to_sponge" />
+      <field id="ulat_sponge_reduced" field_ref="ulat_sponge_interp" grid_ref="reduced_sponge" read_access="false"/>
+      <field id="ulat_sponge_mean"                                   grid_ref="to_sponge" />
+      <field id="ulat_from_sponge" field_ref="ulat_sponge_mean"      grid_ref="from_sponge" read_access="false"/>
+
+      <field id="theta_sponge"                                        grid_ref="sponge_i"/>
+      <field id="theta_sponge_interp" field_ref="theta_sponge"         grid_ref="to_sponge" />
+      <field id="theta_sponge_reduced" field_ref="theta_sponge_interp" grid_ref="reduced_sponge" read_access="false"/>
+      <field id="theta_sponge_mean"                                   grid_ref="to_sponge" />
+      <field id="theta_from_sponge" field_ref="theta_sponge_mean"      grid_ref="from_sponge" read_access="false"/>
+
+      <field id="rhodz_sponge"                                        grid_ref="sponge_i"/>
+      <field id="rhodz_sponge_interp" field_ref="rhodz_sponge"         grid_ref="to_sponge" />
+      <field id="rhodz_sponge_reduced" field_ref="rhodz_sponge_interp" grid_ref="reduced_sponge" read_access="false"/>
+
+     <field id="dulon_sponge" domain_ref="i" axis_ref="lev" />
+     <field id="dulat_sponge" domain_ref="i" axis_ref="lev" />
+     <field id="dtheta_sponge" domain_ref="i" axis_ref="lev" />
+ 
+
+    </field_definition>
+
+    <file_definition>
+      <file id="check_sponge" output_freq="1ts" type="one_file" sync_freq="1ts" enabled="false" >
+        <field_group  prec="8">
+          <field field_ref="ulon_sponge" operation="instant" enabled="true"/>
+          <field field_ref="ulon_sponge_interp" operation="instant" enabled="true"/>
+          <field field_ref="ulon_sponge_reduced" operation="instant" enabled="true"/>
+          <field field_ref="ulon_sponge_mean" operation="instant" enabled="true"/>
+          <field field_ref="ulon_from_sponge" operation="instant" enabled="true"/>
+
+          <field field_ref="ulat_sponge" operation="instant" enabled="true"/>
+          <field field_ref="ulat_sponge_interp" operation="instant" enabled="true"/>
+          <field field_ref="ulat_sponge_reduced" operation="instant" enabled="true"/>
+          <field field_ref="ulat_sponge_mean" operation="instant" enabled="true"/>
+          <field field_ref="ulat_from_sponge" operation="instant" enabled="true"/>
+
+          <field field_ref="theta_sponge" operation="instant" enabled="true"/>
+          <field field_ref="theta_sponge_interp" operation="instant" enabled="true"/>
+          <field field_ref="theta_sponge_reduced" operation="instant" enabled="true"/>
+          <field field_ref="theta_sponge_mean" operation="instant" enabled="true"/>
+          <field field_ref="theta_from_sponge" operation="instant" enabled="true"/>
+          
+          <field field_ref="rhodz_sponge" operation="instant" enabled="true"/>
+          <field field_ref="rhodz_sponge_interp" operation="instant" enabled="true"/>
+          <field field_ref="rhodz_sponge_reduced" operation="instant" enabled="true"/>
+
+         <field field_ref="dulon_sponge" operation="instant" enabled="true"/>
+         <field field_ref="dulat_sponge" operation="instant" enabled="true"/>
+         <field field_ref="dtheta_sponge" operation="instant" enabled="true"/>
+
+        </field_group>
+      </file>
+    </file_definition>
+    
+</context>
+
Index: /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/PARAM/vert.def_L79
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/PARAM/vert.def_L79	(revision 6361)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/GENERAL/PARAM/vert.def_L79	(revision 6362)
@@ -31,5 +31,5 @@
 iflag_sponge=2
 tau_sponge=5.e-4
-mode_sponge=1
+mode_sponge=3
 
 ## Section 3: Following paramereters are shared between LMDZ/dynamics and DYNAMICO
Index: /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/dissip/sponge.f90
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/dissip/sponge.f90	(revision 6362)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/dissip/sponge.f90	(revision 6362)
@@ -0,0 +1,497 @@
+MODULE sponge_mod
+  USE icosa
+
+  PRIVATE
+  
+  REAL,SAVE :: tau_sponge !inverse of charactericstic relaxation time scale at the topmost layer (Hz)
+  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
+  INTEGER,SAVE :: mode_sponge !1 --> u and v relax towards 0
+                              !2 --> u and v relax towards their zonal mean
+                              !3 --> u,v and pot. temp. relax towards their zonal mean
+ 
+  !$OMP THREADPRIVATE(tau_sponge,iflag_sponge,mode_sponge)
+  REAL,ALLOCATABLE,SAVE :: rdamp(:) ! quenching coefficient
+  !$OMP THREADPRIVATE(rdamp)
+  REAL,ALLOCATABLE,SAVE:: lambda(:) ! inverse or quenching time scale (Hz)
+  !$OMP THREADPRIVATE(lambda)
+
+  PUBLIC sponge, init_sponge, iflag_sponge, pre_init_sponge
+
+  INTEGER,SAVE :: ni_glo
+  !$OMP THREADPRIVATE(ni_glo)
+  INTEGER,SAVE :: nj_glo
+  !$OMP THREADPRIVATE(nj_glo)
+  INTEGER,SAVE :: ni
+  !$OMP THREADPRIVATE(ni)
+  INTEGER,SAVE :: nj
+  !$OMP THREADPRIVATE(nj)
+  INTEGER,SAVE :: ibegin
+  !$OMP THREADPRIVATE(ibegin)
+  INTEGER,SAVE :: jbegin
+  !$OMP THREADPRIVATE(jbegin)
+ 
+  TYPE(t_field),POINTER,SAVE :: f_hflux(:)
+  TYPE(t_field),POINTER,SAVE :: f_ulon(:)
+  TYPE(t_field),POINTER,SAVE :: f_ulat(:)
+  TYPE(t_field),POINTER,SAVE :: f_ulon_sponge(:)
+  TYPE(t_field),POINTER,SAVE :: f_ulat_sponge(:)
+  TYPE(t_field),POINTER,SAVE :: f_rhodz_sponge(:)
+  TYPE(t_field),POINTER,SAVE :: f_theta_sponge(:)
+
+  TYPE(t_field),POINTER,SAVE :: f_dulon(:)
+  TYPE(t_field),POINTER,SAVE :: f_dulat(:)
+ 
+  INTEGER,SAVE  :: n_sponge
+  !$OMP THREADPRIVATE(n_sponge)
+  INTEGER,SAVE  :: begin_sponge
+  !$OMP THREADPRIVATE(begin_sponge)
+  
+CONTAINS
+
+  SUBROUTINE pre_init_sponge
+  USE icosa
+  USE disvert_mod
+  IMPLICIT NONE
+    INTEGER :: l
+
+    tau_sponge = 1.e-5
+    CALL getin("tau_sponge",tau_sponge)
+    PRINT*,'tau_sponge = ',tau_sponge
+
+    iflag_sponge = 0
+    CALL getin("iflag_sponge",iflag_sponge)
+    PRINT*,'iflag_sponge = ',iflag_sponge
+    
+    mode_sponge = 1
+    CALL getin("mode_sponge",mode_sponge)
+    PRINT*,'mode_sponge = ',mode_sponge
+
+    IF (iflag_sponge == 0) THEN
+      PRINT*,'init_sponge: no sponge'
+      RETURN
+    ENDIF
+    
+    IF (.NOT. (iflag_sponge == 1 .OR. iflag_sponge == 2 .OR. iflag_sponge == 3)) THEN
+      PRINT*,'Bad selector for variable iflag_sponge : <',iflag_sponge,'> options are 0,1,2'
+      STOP
+    ELSE
+      IF (.NOT. (mode_sponge == 1 .OR. mode_sponge == 2 .OR. mode_sponge == 3)) THEN
+        PRINT*,'Bad selector for variable mode_sponge : <',mode_sponge,'> options are 0,1,2'
+        STOP
+      ENDIF
+    ENDIF
+
+
+    IF (iflag_sponge == 0) RETURN
+    
+    ALLOCATE(rdamp(llm))
+    ALLOCATE(lambda(llm))
+    !$acc enter data create(rdamp(:),lambda(:)) async
+    
+    IF (iflag_sponge == 1) THEN
+! sponge quenching over the topmost 4 atmospheric layers
+        lambda(:)=0.
+        lambda(llm)=tau_sponge
+        lambda(llm-1)=tau_sponge/2.
+        lambda(llm-2)=tau_sponge/4.
+        lambda(llm-3)=tau_sponge/8.
+        n_sponge=4
+        begin_sponge=llm-3
+    ELSE IF (iflag_sponge == 2 ) THEN
+! sponge quenching over topmost layers down to pressures which are
+! higher than 100 times the topmost layer pressure
+        lambda(:)=tau_sponge*max(presnivs(llm)/presnivs(:)-0.01,0.)
+        DO l=llm,1,-1
+          IF (lambda(l)/=0.) begin_sponge=l
+        ENDDO
+        n_sponge=llm-begin_sponge+1
+    ENDIF
+     
+    IF (iflag_sponge>=1 .AND. (mode_sponge==2 .OR. mode_sponge==3)) THEN
+     CALL pre_init_sponge_zonal_mean_u
+    ENDIF
+    
+  END SUBROUTINE pre_init_sponge
+ 
+ 
+  SUBROUTINE init_sponge
+  USE icosa
+  USE disvert_mod
+  USE omp_para
+  USE mpipara, ONLY: is_mpi_master
+  IMPLICIT NONE
+    INTEGER :: l
+  
+    IF (iflag_sponge == 0) RETURN
+
+    rdamp(:)=itau_dissip*lambda(:)
+
+
+    IF (is_mpi_master) THEN
+      PRINT*,'lambda = '
+      DO l=ll_begin,ll_end
+        PRINT*,l,lambda(l)
+      ENDDO
+      PRINT*,'rdamp = '
+      DO l=ll_begin,ll_end
+        PRINT*,l,rdamp(l)
+      ENDDO
+    ENDIF
+    !$acc update device(rdamp(:), lambda(:)) async
+ 
+   IF (iflag_sponge>=1 .AND. (mode_sponge == 2 .OR. mode_sponge == 3)) THEN
+     CALL init_sponge_zonal_mean_u
+   ENDIF
+   
+  END SUBROUTINE init_sponge
+
+  
+  SUBROUTINE sponge(f_mass, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz)
+  USE icosa
+  USE theta2theta_rhodz_mod
+  USE exner_mod
+  USE geopotential_mod
+  USE trace
+  USE time_mod
+  USE omp_para
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_mass(:)
+    TYPE(t_field),POINTER :: f_due(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+
+    REAL(rstd),POINTER         :: due(:,:)!,due_sponge(:,:)
+    REAL(rstd),POINTER         :: ue(:,:)
+    REAL(rstd),POINTER         :: theta_rhodz(:,:,:)
+    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)!,dtheta_sponge(:,:)
+
+    INTEGER :: ind
+    INTEGER :: l,ij
+
+!$OMP BARRIER
+    
+    CALL trace_start("sponge")
+
+    IF (mode_sponge == 1) THEN
+      DO ind=1,ndomain
+        IF (.NOT. assigned_domain(ind)) CYCLE
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+        ue=f_ue(ind) 
+        due=f_due(ind) 
+        theta_rhodz=f_theta_rhodz(ind)
+        dtheta_rhodz=f_dtheta_rhodz(ind)
+
+       !$acc parallel loop collapse(2) present(ue(:,:), due(:,:), dtheta_rhodz)  async
+        DO l=ll_begin,ll_end
+         !$SIMD
+          DO ij=ij_begin,ij_end
+
+            due(ij+u_right,l) = ue(ij+u_right,l)
+            due(ij+u_lup,l)   = ue(ij+u_lup,l)
+            due(ij+u_ldown,l) = ue(ij+u_ldown,l)
+
+            dtheta_rhodz(ij,l,:) = 0.0
+          ENDDO
+        ENDDO
+       !$acc end parallel loop
+      END DO
+    ELSE IF (mode_sponge == 2 .OR. mode_sponge == 3) THEN
+      CALL sponge_zonal_mean_u(mode_sponge, f_mass, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz) 
+    ELSE
+      PRINT*,'Arnaud Bad selector for variable mode_sponge : <',mode_sponge,             &
+            '> options 2 and 3 not available for the moment!'
+      STOP
+    ENDIF
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      due=f_due(ind) 
+      dtheta_rhodz=f_dtheta_rhodz(ind)
+
+      !$acc parallel loop collapse(2) present(ue(:,:), due(:,:), rdamp(:), dtheta_rhodz)  async
+      DO l=ll_begin,ll_end
+      !$SIMD
+        DO ij=ij_begin,ij_end
+          due(ij+u_right,l)        = -rdamp(l)*due(ij+u_right,l)
+          due(ij+u_lup,l)          = -rdamp(l)*due(ij+u_lup,l)
+          due(ij+u_ldown,l)        = -rdamp(l)*due(ij+u_ldown,l)
+          dtheta_rhodz(ij,l,:)     = -rdamp(l)*dtheta_rhodz(ij,l,:)
+        ENDDO
+      ENDDO
+     !$acc end parallel loop
+    END DO
+
+   CALL trace_end("sponge")
+
+!$OMP BARRIER
+  END SUBROUTINE sponge
+
+  SUBROUTINE pre_init_sponge_zonal_mean_u()
+  USE xios_mod
+  USE grid_param, ONLY : iim_glo
+  USE disvert_mod
+  IMPLICIT NONE
+  INTEGER :: nlon
+  INTEGER :: nlat  
+
+!$OMP MASTER
+    nlon=360.*iim_glo/80.
+    nlat=180.*iim_glo/80.
+
+    CALL xios_set_axis_attr("lev_sponge", n_glo=n_sponge, value=presnivs(begin_sponge:))
+    CALL xios_set_domain_attr("to_sponge", ni_glo=nlon, nj_glo=nlat)
+    CALL xios_set_axis_attr("reduced_sponge", n_glo=nlat)
+    CALL xios_set_field_attr("ulon_sponge_reduced", read_access=.TRUE.)
+    CALL xios_set_field_attr("ulon_from_sponge", read_access=.TRUE.)
+    CALL xios_set_field_attr("ulat_sponge_reduced", read_access=.TRUE.)
+    CALL xios_set_field_attr("ulat_from_sponge", read_access=.TRUE.)
+    CALL xios_set_field_attr("rhodz_sponge_reduced", read_access=.TRUE.)
+    IF (mode_sponge==3) CALL xios_set_field_attr("theta_sponge_reduced", read_access=.TRUE.)
+    IF (mode_sponge==3) CALL xios_set_field_attr("theta_from_sponge", read_access=.TRUE.)
+!$OMP END MASTER
+
+  END SUBROUTINE  pre_init_sponge_zonal_mean_u
+
+
+  SUBROUTINE init_sponge_zonal_mean_u()
+  USE xios_mod
+  IMPLICIT NONE
+    
+!$OMP MASTER
+    CALL xios_get_domain_attr("to_sponge", ni_glo=ni_glo, nj_glo=nj_glo,ni=ni, nj=nj, ibegin=ibegin,jbegin=jbegin)
+!$OMP END MASTER
+    CALL allocate_field(f_hflux, field_u, type_real, llm, name='hflux') 
+    CALL allocate_field(f_ulon, field_t, type_real, llm, name='ulon') 
+    CALL allocate_field(f_ulat, field_t, type_real, llm, name='ulat') 
+    CALL allocate_field(f_ulon_sponge, field_t, type_real, n_sponge, name='ulon_sponge') 
+    CALL allocate_field(f_ulat_sponge, field_t, type_real, n_sponge, name='ulat_sponge') 
+    CALL allocate_field(f_rhodz_sponge, field_t, type_real, n_sponge, name='rhodz_sponge') 
+    CALL allocate_field(f_theta_sponge, field_t, type_real, n_sponge, name='theta_sponge') 
+    CALL allocate_field(f_dulon, field_t, type_real, llm,name='dulon_sponge') 
+    CALL allocate_field(f_dulat, field_t, type_real, llm,name='dulat_sponge') 
+
+  END SUBROUTINE init_sponge_zonal_mean_u
+
+  
+  SUBROUTINE sponge_zonal_mean_u(mode_sponge, f_rhodz, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz)
+  USE icosa
+  USE theta2theta_rhodz_mod
+  USE exner_mod
+  USE geopotential_mod
+  USE trace
+  USE time_mod
+  USE omp_para
+  USE wind_mod
+  USE xios_mod
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)    :: mode_sponge
+    TYPE(t_field),POINTER :: f_rhodz(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_due(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+
+    REAL(rstd),POINTER         :: due(:,:)!,due_sponge(:,:)
+    REAL(rstd),POINTER         :: ue(:,:)
+    REAL(rstd),POINTER         :: rhodz(:,:)
+    REAL(rstd),POINTER         :: hflux(:,:)
+    REAL(rstd),POINTER         :: theta_rhodz(:,:,:)
+    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)!,dtheta_sponge(:,:)
+    REAL(rstd),POINTER         :: ulon(:,:)
+    REAL(rstd),POINTER         :: ulat(:,:)
+    REAL(rstd),POINTER         :: ulon_sponge(:,:)
+    REAL(rstd),POINTER         :: ulat_sponge(:,:)
+    REAL(rstd),POINTER         :: theta_sponge(:,:)
+    REAL(rstd),POINTER         :: rhodz_sponge(:,:)
+    
+    REAL(rstd) :: ulon_zonal_mean(nj_glo,n_sponge)
+    REAL(rstd) :: ulat_zonal_mean(nj_glo,n_sponge)
+    REAL(rstd) :: theta_zonal_mean(nj_glo,n_sponge)
+    REAL(rstd) :: rhodz_zonal_mean(nj_glo,n_sponge)
+    REAL(rstd) :: domain_zonal_mean(ni,nj,n_sponge)
+   
+    INTEGER :: i,j,ind,ij,l,ll 
+
+    DO ind = 1 , ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ue=f_ue(ind)
+      rhodz=f_rhodz(ind)
+      hflux = f_hflux(ind)
+         
+      DO l = ll_begin, ll_end
+         !DIR$ IVDEP
+         DO ij=ij_begin_ext,ij_end_ext         
+            hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*ue(ij+u_right,l)
+            hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*ue(ij+u_lup,l)
+            hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*ue(ij+u_ldown,l)
+         ENDDO
+       ENDDO
+    ENDDO
+
+              
+    CALL transfert_request(f_hflux,req_e1_vect)
+    CALL un2ulonlat(f_hflux, f_ulon, f_ulat)
+
+    DO ind = 1 , ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      
+      ulon=f_ulon(ind)
+      ulat=f_ulat(ind)
+      ulon_sponge = f_ulon_sponge(ind)
+      ulat_sponge = f_ulat_sponge(ind)
+      rhodz_sponge = f_rhodz_sponge(ind)
+      theta_sponge = f_theta_sponge(ind)
+      theta_rhodz = f_theta_rhodz(ind)
+              
+      DO l = 1, llm
+         IF (l<begin_sponge) CYCLE
+         ll=l-begin_sponge+1
+         !DIR$ IVDEP
+         DO ij=ij_begin_ext,ij_end_ext         
+            ulon_sponge(ij,ll)= ulon(ij,l) 
+            ulat_sponge(ij,ll)= ulat(ij,l) 
+            rhodz_sponge(ij,ll) = rhodz(ij,l)
+            IF (mode_sponge==3) theta_sponge(ij,ll) = theta_rhodz(ij,l,1)
+         ENDDO
+       ENDDO
+    ENDDO
+
+!$OMP BARRIER
+    
+    CALL xios_write_field("ulon_sponge",f_ulon_sponge)
+    CALL xios_write_field("ulat_sponge",f_ulat_sponge)
+    CALL xios_write_field("rhodz_sponge",f_rhodz_sponge)
+!$OMP MASTER
+    CALL xios_recv_field("ulon_sponge_reduced",ulon_zonal_mean)
+    CALL xios_recv_field("ulat_sponge_reduced",ulat_zonal_mean)
+    CALL xios_recv_field("rhodz_sponge_reduced",rhodz_zonal_mean)
+!$OMP END MASTER
+
+    IF (mode_sponge==3) THEN
+      CALL xios_write_field("theta_sponge",f_theta_sponge)
+!$OMP MASTER
+      CALL xios_recv_field("theta_sponge_reduced",theta_zonal_mean)
+!$OMP END MASTER
+    ENDIF
+    
+!$OMP MASTER
+    DO j=1,nj
+      DO i=1,ni
+        domain_zonal_mean(i,j,:)=ulon_zonal_mean(jbegin+j,:)/rhodz_zonal_mean(jbegin+j,:)
+      ENDDO
+    ENDDO
+    
+    CALL xios_send_field("ulon_sponge_mean",domain_zonal_mean)
+!$OMP END MASTER
+    CALL xios_read_field("ulon_from_sponge",f_ulon_sponge)
+  
+!$OMP MASTER
+    DO j=1,nj
+      DO i=1,ni
+        domain_zonal_mean(i,j,:)=ulat_zonal_mean(jbegin+j,:)/rhodz_zonal_mean(jbegin+j,:)
+      ENDDO
+    ENDDO
+    
+    CALL xios_send_field("ulat_sponge_mean",domain_zonal_mean)
+!$OMP END MASTER
+    CALL xios_read_field("ulat_from_sponge",f_ulat_sponge)
+    
+    IF (mode_sponge==3) THEN
+!$OMP MASTER
+      DO j=1,nj
+        DO i=1,ni
+          domain_zonal_mean(i,j,:)=theta_zonal_mean(jbegin+j,:)
+        ENDDO
+      ENDDO
+      CALL xios_send_field("theta_sponge_mean",domain_zonal_mean)
+!$OMP END MASTER
+      CALL xios_read_field("theta_from_sponge",f_theta_sponge)
+    ENDIF
+    
+
+    DO ind = 1 , ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      
+      ulon=f_ulon(ind)
+      ulat=f_ulat(ind)
+      ulon_sponge = f_ulon_sponge(ind)
+      ulat_sponge = f_ulat_sponge(ind)
+      theta_sponge = f_theta_sponge(ind)
+      dtheta_rhodz = f_dtheta_rhodz(ind)    
+      theta_rhodz = f_theta_rhodz(ind)
+      DO l = ll_begin, ll_end
+        IF (l<begin_sponge) CYCLE
+        ll=l-begin_sponge+1
+        IF (mode_sponge==2) THEN
+          DO ij=ij_begin_ext,ij_end_ext         
+            ulon(ij,l) = ulon_sponge(ij,ll) 
+            ulat(ij,l) = ulat_sponge(ij,ll) 
+            dtheta_rhodz(ij,l,1)=theta_rhodz(ij,ll,1) 
+          ENDDO
+        ELSE  IF (mode_sponge==3) THEN
+          DO ij=ij_begin_ext,ij_end_ext         
+            ulon(ij,l) = ulon_sponge(ij,ll) 
+            ulat(ij,l) = ulat_sponge(ij,ll) 
+            dtheta_rhodz(ij,l,1) = theta_sponge(ij,ll)
+          ENDDO
+        ENDIF
+      ENDDO
+    ENDDO
+    
+    CALL transfert_request(f_ulon,req_i0)
+    CALL transfert_request(f_ulat,req_i0)
+    IF (mode_sponge==3) THEN
+      CALL transfert_request(f_dtheta_rhodz,req_i0)
+    ENDIF 
+    CALL ulonlat2un(f_ulon,f_ulat,f_due)
+!$OMP BARRIER
+   
+    DO ind = 1 , ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      due=f_due(ind)
+      ue=f_ue(ind)
+      dtheta_rhodz=f_dtheta_rhodz(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+ 
+       DO l = ll_begin, ll_end
+         IF (l<begin_sponge) THEN
+           !DIR$ IVDEP
+           DO ij=ij_begin,ij_end         
+              due(ij+u_right,l) = 0.
+              due(ij+u_lup,l)   = 0.
+              due(ij+u_ldown,l) = 0.
+              dtheta_rhodz(ij,l,:) = 0.
+           ENDDO         
+         ELSE
+           !DIR$ IVDEP
+           DO ij=ij_begin,ij_end         
+             due(ij+u_right,l) = ue(ij+u_right,l) - due(ij+u_right,l)
+             due(ij+u_lup,l)   = ue(ij+u_lup,l)   - due(ij+u_lup,l)
+             due(ij+u_ldown,l) = ue(ij+u_ldown,l) - due(ij+u_ldown,l)
+             dtheta_rhodz(ij,l,1) = theta_rhodz(ij,l,1) - dtheta_rhodz(ij,l,1)
+          ENDDO
+        ENDIF
+      ENDDO
+    ENDDO  
+   
+    CALL transfert_request(f_due,req_e1_vect)
+!ym to check results
+!    CALL un2ulonlat(f_due, f_dulon, f_dulat)
+!    CALL xios_write_field("dulon_sponge",f_dulon)
+!    CALL xios_write_field("dulat_sponge",f_dulat)
+!    CALL xios_write_field("dtheta_sponge",f_dtheta_rhodz)
+    
+  END SUBROUTINE sponge_zonal_mean_u
+    
+
+END MODULE sponge_mod
Index: /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/icosa_init.f90
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/icosa_init.f90	(revision 6362)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/icosa_init.f90	(revision 6362)
@@ -0,0 +1,144 @@
+MODULE icosa_init_mod
+  
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: icosa_init
+  
+CONTAINS
+  
+  SUBROUTINE icosa_init
+    USE profiling_mod,        ONLY : init_profiling
+    USE mpipara,              ONLY : is_mpi_master, is_mpi_root, mpi_rank, init_mpipara, finalize_mpipara
+    USE signature_mod,        ONLY : print_signature
+    USE trace,                ONLY : trace_off, init_trace
+    USE xios_mod,             ONLY : xios_init
+    USE geometry,             ONLY : geom
+    USE earth_const,          ONLY : init_earth_const
+    USE grid_param    ,       ONLY : grid_type, grid_ico
+    USE init_grid_param_mod,  ONLY : init_grid_param
+    USE omp_para,             ONLY : init_omp_para, switch_omp_no_distrib_level
+    USE metric,               ONLY : init_metric
+    USE domain_mod,           ONLY : compute_domain
+    USE transfert_mod,        ONLY : init_transfert
+    USE write_field_mod,      ONLY : init_writefield, writefield, close_files
+    USE trace,                ONLY : init_trace
+    USE read_metric_mod,      ONLY : read_metric
+    USE sponge_mod,           ONLY : pre_init_sponge
+    USE geometry,             ONLY : compute_geometry
+    USE disvert_mod,          ONLY : init_disvert, write_apbp
+    USE vertical_interp_mod,  ONLY : init_vertical_interp
+    USE tracer_icosa_mod,     ONLY : init_tracer
+    USE time_mod,             ONLY : init_time, close_time_counter
+    USE restart_mod,          ONLY : init_restart
+    USE etat0_mod,            ONLY : init_etat0
+    USE output_field_mod,     ONLY : output_field_init, output_field_finalize
+    USE timeloop_gcm_mod,     ONLY : init_timeloop, timeloop
+    USE physics_mod,          ONLY : init_physics, zero_du_phys
+    USE diagflux_mod,         ONLY : init_diagflux
+    USE guided_mod,           ONLY : pre_init_guided
+    
+    CALL init_profiling
+    CALL init_mpipara
+    IF(is_mpi_master) CALL print_signature
+    CALL trace_off
+    CALL xios_init
+    CALL init_earth_const 
+    CALL init_grid_param(is_mpi_master, mpi_rank) ! sets max_primal_deg, etc.
+    CALL init_omp_para(is_mpi_master)
+    IF(grid_type == grid_ico) CALL init_metric
+    CALL compute_domain ! allocates domain_glo(ndomain_glo), domain(ndomain)
+    CALL init_transfert
+    
+    CALL init_writefield
+    CALL init_trace
+    
+    !$OMP PARALLEL  
+    CALL switch_omp_no_distrib_level
+
+    IF(grid_type == grid_ico) CALL read_metric
+    CALL compute_geometry
+    
+    IF(grid_type == grid_ico) CALL check_total_area
+    
+    CALL init_disvert 
+    CALL init_vertical_interp
+    CALL init_tracer
+    
+    CALL WriteField("Ai",geom%Ai)
+    
+    IF (is_mpi_root) CALL write_apbp
+    CALL init_time
+    CALL init_etat0
+    CALL init_restart
+    CALL pre_init_sponge
+    CALL pre_init_guided ! must be done before closing xios context, i.e. before output_field_init
+    CALL output_field_init
+    CALL init_timeloop
+    
+    !$OMP END PARALLEL
+    
+    CALL init_physics
+    
+    !$OMP PARALLEL
+    
+    CALL init_diagflux
+    CALL zero_du_phys
+    
+    CALL timeloop
+    
+    CALL switch_omp_no_distrib_level
+    !$OMP END PARALLEL
+    
+    CALL output_field_finalize
+    CALL close_files
+    CALL close_time_counter
+    CALL finalize_mpipara
+    
+  END SUBROUTINE icosa_init
+  
+  SUBROUTINE check_total_area
+    USE prec,       ONLY : rstd
+    USE domain_mod, ONLY : ndomain, domain
+    USE geometry,   ONLY : swap_geometry, Ai, ne
+    USE mpipara,    ONLY : is_mpi_root
+    USE dimensions
+    USE field_mod
+    USE domain_param
+    USE omp_para
+    USE icosa
+    
+    TYPE(t_field), POINTER, SAVE :: sum_ne(:)
+    REAL(rstd), POINTER          :: pt_sum_ne(:)
+    INTEGER :: ind,i,j,k,n 
+    REAL(rstd) :: tot_sum=0
+    
+    CALL allocate_field(sum_ne,field_T,type_real)
+    
+    !$OMP BARRIER
+    !$OMP MASTER
+    DO ind=1,ndomain
+       
+       pt_sum_ne=sum_ne(ind)
+       
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             n=(j-1)*iim+i
+             pt_sum_ne(n)=0
+             DO k=1,6
+                pt_sum_ne(n)=pt_sum_ne(n)+ne(n,k)
+             ENDDO
+             IF (domain(ind)%own(i,j)) tot_sum=tot_sum+Ai(n)
+          ENDDO
+       ENDDO
+    ENDDO
+    
+    IF (is_mpi_root) PRINT *," Diff surf",1-tot_sum/(4*Pi*radius*radius)
+    !$OMP END MASTER
+    
+  END SUBROUTINE check_total_area
+  
+END MODULE icosa_init_mod
Index: /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/time/timeloop_gcm.F90
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/time/timeloop_gcm.F90	(revision 6362)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/time/timeloop_gcm.F90	(revision 6362)
@@ -0,0 +1,535 @@
+MODULE timeloop_gcm_mod
+  USE profiling_mod
+  USE icosa
+  USE disvert_mod
+  USE trace
+  USE omp_para
+  USE euler_scheme_mod
+  USE explicit_scheme_mod
+  USE hevi_scheme_mod
+  IMPLICIT NONE
+  PRIVATE
+  
+  INTEGER, PARAMETER :: sync_it=10
+  TYPE(t_message),SAVE :: req_ps0, req_rhodz0, req_mass0, req_theta_rhodz0, req_u0, req_q0, req_W0, req_geopot0
+  LOGICAL, SAVE :: positive_theta
+  INTEGER :: itau_prof, id_timeloop, id_dyn, id_phys, id_dissip, id_adv, id_diags
+  PUBLIC :: init_timeloop, timeloop
+
+CONTAINS
+  
+  SUBROUTINE init_timeloop
+    USE dissip_gcm_mod
+    USE observable_mod
+    USE caldyn_mod
+    USE etat0_mod
+    USE guided_mod
+    USE advect_tracer_mod
+    USE check_conserve_mod
+    USE output_field_mod
+    USE theta2theta_rhodz_mod
+    USE sponge_mod
+
+    CHARACTER(len=255) :: def
+
+    CALL register_id('timeloop',id_timeloop)
+    CALL register_id('dyn',id_dyn)
+    CALL register_id('dissip',id_dissip)
+    CALL register_id('phys',id_phys)
+    CALL register_id('adv',id_adv)
+    CALL register_id('diags',id_diags)
+
+    CALL init_caldyn
+    
+!    IF (xios_output) itau_out=1
+    IF (.NOT. enable_io) itau_out=HUGE(itau_out)
+
+    itau_prof=1000
+    CALL getin('itau_profiling',itau_prof)
+
+    positive_theta=.FALSE.
+    CALL getin('positive_theta',positive_theta)
+    IF(positive_theta .AND. nqtot<1) THEN
+       PRINT *, 'nqtot must be >0 if positive_theta is .TRUE.'
+       STOP
+    END IF
+
+    def='ARK2.3'
+    CALL getin('time_scheme',def)
+    SELECT CASE (TRIM(def))
+    CASE('euler')
+       scheme_family=explicit
+       scheme=euler
+       nb_stage=1
+    CASE ('runge_kutta')
+       scheme_family=explicit
+       scheme=rk4
+       nb_stage=4
+    CASE ('RK2.5')
+       scheme_family=explicit
+       scheme=rk25
+       nb_stage=5
+    CASE ('leapfrog_matsuno')
+       scheme_family=explicit
+       scheme=mlf
+       matsuno_period=5
+       CALL getin('matsuno_period',matsuno_period)
+       nb_stage=matsuno_period+1
+    CASE('ARK2.3')
+       scheme_family=hevi
+       scheme=ark23
+       nb_stage=3
+       CALL set_coefs_ark23(dt)
+    CASE('ARK3.3')
+       scheme_family=hevi
+       scheme=ark33
+       nb_stage=3
+       CALL set_coefs_ark33(dt)
+    CASE ('none')
+       nb_stage=0
+    CASE default
+       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             &
+            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>,<RK2.5>,<ARK2.3>'
+       STOP
+    END SELECT
+    
+    IF (scheme_family /= hevi) THEN
+       CALL abort_acc("scheme_family /= hevi")
+    END IF
+
+    ! Time-independant orography
+    CALL allocate_field(f_phis,field_t,type_real,name='phis')
+    ! Model state at current time step
+    CALL allocate_field(f_ps,field_t,type_real, name='ps')
+    CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
+    CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
+    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,nqdyn,name='theta_rhodz')
+    CALL allocate_field(f_u,field_u,type_real,llm,name='u')
+    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')
+    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W') ! used only if .not. hydrostatic
+    CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q')
+    ! Mass fluxes
+    CALL allocate_field(f_hflux,field_u,type_real,llm, ondevice=.TRUE.)    ! instantaneous mass fluxes
+    CALL allocate_field(f_hfluxt,field_u,type_real,llm,ondevice=.TRUE.)   ! mass "fluxes" accumulated in time
+    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes
+    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt',ondevice=.TRUE.)
+    
+    SELECT CASE(scheme_family)
+    CASE(explicit)
+       ! Trends
+       CALL allocate_field(f_dps,field_t,type_real,name='dps')
+       CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
+       CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,nqdyn,name='dtheta_rhodz')
+       CALL allocate_field(f_du,field_u,type_real,llm,name='du')
+       ! Model state at previous time step (RK/MLF)
+       CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
+       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
+       CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,nqdyn,name='theta_rhodzm1')
+       CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')
+    CASE(hevi)
+       ! Trends
+       CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_dW_slow, field_t,type_real,llm+1,name='dW_slow')
+       CALL allocate_fields(nb_stage,f_dW_fast, field_t,type_real,llm+1,name='dW_fast')
+       CALL allocate_fields(nb_stage,f_dPhi_slow, field_t,type_real,llm+1,name='dPhi_slow')
+       CALL allocate_fields(nb_stage,f_dPhi_fast, field_t,type_real,llm+1,name='dPhi_fast')
+       f_dps => f_dps_slow(:,1)
+       f_du => f_du_slow(:,1)
+       f_dtheta_rhodz => f_dtheta_rhodz_slow(:,1)
+    END SELECT
+
+    SELECT CASE(scheme)
+    CASE(mlf)
+       ! Model state 2 time steps ago (MLF)
+       CALL allocate_field(f_psm2,field_t,type_real)
+       CALL allocate_field(f_massm2,field_t,type_real,llm)
+       CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm,nqdyn)
+       CALL allocate_field(f_um2,field_u,type_real,llm)
+    END SELECT
+
+    CALL init_theta2theta_rhodz
+    CALL init_dissip
+    CALL init_sponge
+    CALL init_observable
+    CALL init_guided
+    CALL init_advect_tracer
+    CALL init_check_conserve
+
+    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
+
+    CALL transfert_request(f_phis,req_i0) 
+    CALL transfert_request(f_phis,req_i1) 
+
+    CALL init_message(f_ps,req_i0,req_ps0)
+    CALL init_message(f_mass,req_i0,req_mass0)
+    CALL init_message(f_rhodz,req_i0,req_rhodz0)
+    CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
+    CALL init_message(f_u,req_e0_vect,req_u0)
+    CALL init_message(f_q,req_i0,req_q0)
+    CALL init_message(f_geopot,req_i0,req_geopot0)
+    CALL init_message(f_W,req_i0,req_W0)
+
+  END SUBROUTINE init_timeloop
+  
+  SUBROUTINE timeloop
+    USE abort_mod
+    USE dissip_gcm_mod
+    USE sponge_mod
+    USE observable_mod
+    USE etat0_mod
+    USE guided_mod
+    USE caldyn_mod
+    USE advect_tracer_mod
+    USE diagflux_mod
+    USE physics_mod
+    USE mpipara
+    USE transfert_mod
+    USE check_conserve_mod
+    USE xios_mod
+    USE output_field_mod
+    USE write_etat0_mod
+    USE restart_mod
+    USE checksum_mod
+    USE explicit_scheme_mod
+    USE hevi_scheme_mod
+    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:)
+
+    REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out
+    INTEGER :: ind, it,l
+    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time
+    LOGICAL, PARAMETER :: check_rhodz=.FALSE.
+    INTEGER(kind=8) :: start_clock, stop_clock, rate_clock
+    INTEGER :: itau_sync   ! best iteration for synchronisation and ensure 1+1=2
+    INTEGER :: i
+    
+    LOGICAL,SAVE :: first_physic=.TRUE.
+    !$OMP THREADPRIVATE(first_physic)    
+
+    itau_sync=1
+    DO i=2,3*sync_it
+      IF (MOD(86400,INT(i*dt))==0 .AND. ABS((sync_it-itau_sync)*1./sync_it )/sync_it < (sync_it-itau_sync)*1./sync_it) itau_sync=i
+    ENDDO
+    CALL getin("itau_sync",itau_sync)
+    IF (is_master) PRINT*,"Synchronize frontier every itau_sync =",itau_sync
+      
+    CALL switch_omp_distrib_level
+    CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces
+
+    !$OMP BARRIER
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
+       IF(caldyn_eta==eta_mass) THEN
+          CALL compute_rhodz(.TRUE., ps, rhodz, ondevice=.FALSE.) ! save rhodz for transport scheme before dynamics update ps
+       ELSE
+          DO l=ll_begin,ll_end
+             rhodz(:,l)=mass(:,l)
+          ENDDO
+       END IF
+    END DO
+    !$OMP BARRIER
+
+    fluxt_zero=.TRUE.
+
+    IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
+    IF(diagflux_on) THEN
+       adv_over_out = itau_adv*(1./itau_out)
+    ELSE
+       adv_over_out = 0.
+    END IF
+
+    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,itau0)  
+
+    Call trace_on
+
+    IF (xios_output) THEN ! we must call update_calendar before any XIOS output
+      IF (is_omp_master) CALL xios_update_calendar(1)
+    END IF
+    
+!    IF (is_write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
+    IF (is_write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W)
+   
+    CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
+
+    !$OMP MASTER
+    CALL SYSTEM_CLOCK(start_clock, rate_clock)
+    !$OMP END MASTER   
+    call update_device_field(f_ps)
+    call update_device_field(f_mass)
+    CALL update_device_field(f_theta_rhodz)
+    CALL update_device_field(f_u)
+    CALL update_device_field(f_q)
+    CALL update_device_field(f_geopot)
+    CALL update_device_field(f_wflux)
+    CALL update_device_field(f_rhodz)
+    call reset_profiling()  
+
+
+    DO it=itau0+1,itau0+itaumax
+
+       CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)
+
+       CALL enter_profile(id_timeloop)
+
+       IF (xios_output) THEN
+
+          IF(it>itau0+1 .AND. is_omp_master) CALL xios_update_calendar(it-itau0)
+       ELSE
+          CALL update_time_counter(dt*it)
+       END IF
+
+       IF (it==itau0+1 .OR. MOD(it-1,itau_sync)==0) THEN
+          CALL send_message(f_ps,req_ps0)
+          CALL wait_message(req_ps0)
+          CALL send_message(f_mass,req_mass0)
+          CALL wait_message(req_mass0)
+          CALL send_message(f_rhodz,req_rhodz0)
+          CALL wait_message(req_rhodz0)
+          CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
+          CALL wait_message(req_theta_rhodz0)
+          CALL send_message(f_u,req_u0)
+          CALL wait_message(req_u0)
+          CALL send_message(f_q,req_q0) 
+          CALL wait_message(req_q0)
+          IF(.NOT. hydrostatic) THEN
+             CALL send_message(f_geopot,req_geopot0)
+             CALL wait_message(req_geopot0)
+             CALL send_message(f_W,req_W0)
+             CALL wait_message(req_W0)
+          END IF
+       ENDIF
+
+!       CALL checksum(f_ps)
+!       CALL checksum(f_theta_rhodz)
+!       CALL checksum(f_u)
+!       CALL checksum(f_q)
+!       CALL checksum(f_mass)
+!       CALL checksum(f_geopot)
+!       CALL checksum(f_rhodz)
+!       CALL checksum(f_wflux)
+
+       CALL guided(it*dt, f_ps, f_rhodz, f_theta_rhodz,f_u,f_q)
+       
+       CALL enter_profile(id_dyn)
+       SELECT CASE(scheme_family)
+       CASE(explicit)
+          CALL abort_acc("explicit_scheme")
+          CALL explicit_scheme(it, fluxt_zero)
+       CASE(hevi)
+          CALL HEVI_scheme(it, fluxt_zero)
+       END SELECT
+       CALL exit_profile(id_dyn)
+
+!       CALL checksum(f_ps)
+!       CALL checksum(f_theta_rhodz)
+!       CALL checksum(f_u)
+!       CALL checksum(f_q)
+!       CALL checksum(f_mass)
+!       CALL checksum(f_geopot)
+!       CALL checksum(f_rhodz)
+!       CALL checksum(f_wflux)
+
+       CALL enter_profile(id_dissip)
+       IF (MOD(it,itau_dissip)==0) THEN
+          
+          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+          
+          CALL enter_profile(id_diags)
+          CALL check_conserve_detailed(it, AAM_dyn, &
+               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL exit_profile(id_diags)
+
+          CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du)
+          
+          CALL euler_scheme(.FALSE.)  ! update only u, theta
+          IF (iflag_sponge > 0) THEN
+             CALL sponge(f_mass,f_u,f_du,f_theta_rhodz,f_dtheta_rhodz)
+             CALL euler_scheme(.FALSE.)  ! update only u, theta
+          END IF
+
+          CALL update_border(it*dt,f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+
+          CALL enter_profile(id_diags)
+          CALL check_conserve_detailed(it, AAM_dissip, &
+               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL exit_profile(id_diags)
+       END IF
+       CALL exit_profile(id_dissip)
+!       CALL checksum(f_ps)
+!       CALL checksum(f_theta_rhodz)
+!       CALL checksum(f_u)
+!       CALL checksum(f_q)
+!       CALL checksum(f_hfluxt)
+!       CALL checksum(f_wfluxt)
+!       CALL checksum(f_u)
+!       CALL checksum(f_rhodz)
+
+       CALL enter_profile(id_adv)
+       IF(MOD(it,itau_adv)==0) THEN
+          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
+               adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt)  ! accumulate mass and fluxes if diagflux_on
+          fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
+          ! At this point advect_tracer has obtained the halos of u and rhodz,
+          ! needed for correct computation of kinetic energy
+          IF(diagflux_on) CALL abort_acc("diagflux_on")
+          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta,f_buf_i, f_hfluxt)
+
+          DO ind=1,ndomain
+            IF (.NOT. assigned_domain(ind)) CYCLE
+            CALL swap_dimensions(ind)
+            CALL swap_geometry(ind)
+            rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
+            IF(caldyn_eta==eta_mass) THEN
+              CALL compute_rhodz(.TRUE., ps, rhodz, ondevice=.TRUE.) ! save rhodz for transport scheme before dynamics update ps
+            ELSE
+              DO l=ll_begin,ll_end
+                 rhodz(:,l)=mass(:,l)
+              ENDDO
+            END IF
+          END DO
+
+          IF(positive_theta) CALL abort_acc("positive_theta")
+          IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
+          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+       END IF
+       CALL exit_profile(id_adv)
+!       CALL checksum(f_ps)
+!       CALL checksum(f_theta_rhodz)
+!       CALL checksum(f_u)
+!       CALL checksum(f_q)
+
+       CALL enter_profile(id_diags)
+!       IF (MOD(it,itau_physics)==0) THEN
+          CALL check_conserve_detailed(it, AAM_dyn, &
+            f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL enter_profile(id_phys)
+          CALL physics(it,f_phis, f_geopot, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
+          CALL exit_profile(id_phys)
+          CALL check_conserve_detailed(it, AAM_phys, &
+               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          !$OMP MASTER
+          IF (first_physic) CALL SYSTEM_CLOCK(start_clock)
+          !$OMP END MASTER   
+          first_physic=.FALSE.
+!       END IF
+
+       IF (MOD(it,itau_check_conserv)==0) THEN
+          CALL update_host_field(f_ps)
+          CALL update_host_field(f_theta_rhodz)
+          CALL update_host_field(f_u)
+          CALL update_host_field(f_dps)
+          CALL update_host_field(f_q)
+          CALL check_conserve_detailed(it, AAM_dyn, &
+               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,it) 
+       ENDIF
+
+       IF (mod(it,itau_out)==0 ) THEN
+ 
+          CALL update_border(it*dt,f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
+          
+          CALL transfert_request(f_u,req_e1_vect)
+          CALL update_host_field(f_ps)              
+          CALL update_host_field(f_mass)
+          CALL update_host_field(f_theta_rhodz)
+          CALL update_host_field(f_geopot)
+          CALL update_host_field(f_u)
+          CALL update_host_field(f_q)
+          CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
+       ENDIF
+       CALL exit_profile(id_diags)
+
+       CALL exit_profile(id_timeloop)
+    END DO
+    
+    CALL update_host_field(f_ps)
+    CALL update_host_field(f_theta_rhodz)
+    CALL update_host_field(f_u)
+    CALL update_host_field(f_q)
+    CALL update_host_field(f_geopot)
+
+!    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q) 
+    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W) 
+
+    CALL update_host_field(f_dps)    
+    CALL check_conserve_detailed(it, AAM_dyn, &
+         f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,it)  
+    
+    !$OMP MASTER
+    CALL SYSTEM_CLOCK(stop_clock)
+    
+    IF (mpi_rank==0) THEN 
+       PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock 
+    ENDIF
+    !$OMP END MASTER 
+    
+    ! CONTAINS
+  END SUBROUTINE timeloop
+
+  SUBROUTINE print_iteration(it,itau0,itaumax,start_clock,rate_clock)
+    INTEGER :: it, itau0, itaumax, throughput
+    INTEGER(kind=8) :: start_clock, stop_clock, rate_clock
+    REAL :: per_step,total, elapsed
+    
+    IF(is_master) THEN
+       WRITE(*,'(A,I7,A,F14.1)') "It No :",it,"   t :",dt*it
+       IF(MOD(it,10)==0) THEN
+          CALL SYSTEM_CLOCK(stop_clock)
+          elapsed = (stop_clock-start_clock)*1./rate_clock
+          per_step = elapsed/(it-itau0)
+          throughput = INT(dt/per_step)
+          total = per_step*itaumax
+          WRITE(*,'(A,I5,A,F6.2,A,I6)') 'Time spent (s):',INT(elapsed), &
+               '  -- ms/step : ', 1000*per_step, &
+               '  -- Throughput :', throughput
+          WRITE(*,'(A,I5,A,I5)') 'Whole job (min) :', INT(total/60.), &
+               '  -- Completion in (min) : ', INT((total-elapsed)/60.)
+       END IF
+    END IF
+    IF(MOD(it,itau_prof)==0) CALL print_profile
+
+  END SUBROUTINE print_iteration
+
+  SUBROUTINE copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
+    INTEGER :: ind, iq
+    DO ind=1, ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       rhodz=f_rhodz(ind)
+       q=f_q(ind)
+       DO iq=1, nqdyn
+          q(:,:,iq)  = theta_rhodz(:,:,iq)/rhodz(:,:)
+       END DO
+    END DO
+  END SUBROUTINE copy_theta_to_q
+
+  SUBROUTINE copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
+    INTEGER :: ind, iq
+    DO ind=1, ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       rhodz=f_rhodz(ind)
+       q=f_q(ind)
+       DO iq=1,nqdyn
+          theta_rhodz(:,:,iq) = rhodz(:,:)*q(:,:,iq)
+       END DO
+    END DO
+  END SUBROUTINE copy_q_to_theta
+
+END MODULE timeloop_gcm_mod
Index: /CONFIG/UNIFORM/v7/IPSLCM7/compile_ipslcm7.sh
===================================================================
--- /CONFIG/UNIFORM/v7/IPSLCM7/compile_ipslcm7.sh	(revision 6361)
+++ /CONFIG/UNIFORM/v7/IPSLCM7/compile_ipslcm7.sh	(revision 6362)
@@ -352,4 +352,16 @@
 echo >> $outfile ; echo " NOW COMPILE DYNAMICO"   >> $outfile 
 
+        if [ -d $modipsl/config/IPSLCM7/SOURCES/DYNAMICO ] ; then
+         cd $modipsl/config/IPSLCM7/SOURCES/DYNAMICO
+         for dir in $( find . -mindepth 1 -maxdepth 1 -type d ); do
+             for file in $( find $dir -type f ! -path '*.svn*' ); do
+                 echo cp $file $modipsl/modeles/DYNAMICO/$file >> $outfile
+                 cp $file $modipsl/modeles/DYNAMICO/$file  >> $outfile 2>&1
+             done
+         done
+         cd $modipsl/modeles/DYNAMICO
+        fi
+
+
 echo ./make_icosa -$optmode -with_fcm1 -parallel $parallel -external_ioipsl -with_oasis -with_xios -arch $fcm_arch -arch_path $arch_path -job 8 $full_dyna    >> $outfile 
      ./make_icosa -$optmode -with_fcm1 -parallel $parallel -external_ioipsl -with_oasis -with_xios -arch $fcm_arch -arch_path $arch_path -job 8 $full_dyna    >> $outfile 2>&1
