Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 517)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 519)
@@ -2,4 +2,5 @@
   USE icosa
   USE transfert_mod
+  USE caldyn_kernels_hevi_mod
   USE caldyn_kernels_base_mod
   USE caldyn_kernels_mod
@@ -190,5 +191,5 @@
 
 ! temporary shared variable
-    REAL(rstd),POINTER  :: theta(:,:)  
+    REAL(rstd),POINTER  :: theta(:,:,:)  
     REAL(rstd),POINTER  :: pk(:,:)
     REAL(rstd),POINTER  :: geopot(:,:)
@@ -225,10 +226,18 @@
     IF(caldyn_eta==eta_mass) THEN
        CALL send_message(f_ps,req_ps) ! COM00
+       CALL wait_message(req_ps) ! COM00
     ELSE
        CALL send_message(f_mass,req_mass) ! COM00
+       CALL wait_message(req_mass) ! COM00
     END IF
     
     CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
+    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
     CALL send_message(f_u,req_u) ! COM02
+    CALL wait_message(req_u) ! COM02
+
+    IF(.NOT.hydrostatic) THEN
+       STOP 'caldyn_gcm may not be used yet when non-hydrostatic'
+    END IF
 
     SELECT CASE(caldyn_conserv)
@@ -245,11 +254,23 @@
           qu=f_qu(ind)
           qv=f_qv(ind)
+          pk = f_pk(ind)
+          geopot = f_geopot(ind)  
+          hflux=f_hflux(ind)
+          convm = f_dmass(ind)
+          dtheta_rhodz=f_dtheta_rhodz(ind)
+          du=f_du(ind)
           CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02
-       ENDDO
-!       CALL checksum(f_mass)
-!       CALL checksum(f_theta)
-
-       CALL send_message(f_qu,req_qu) ! COM03 wait_message in caldyn_horiz
-!       CALL wait_message(req_qu)
+!          CALL compute_theta(ps,theta_rhodz, mass,theta)
+!          CALL compute_pvort_only(u,mass,qu,qv)
+
+          CALL compute_geopot(mass,theta, ps,pk,geopot)
+!          du(:,:)=0.
+!          CALL compute_caldyn_fast(0.,u,mass,theta,pk,geopot,du)
+       ENDDO       
+
+       CALL send_message(f_u,req_u) ! COM02
+       CALL wait_message(req_u) ! COM02
+       CALL send_message(f_qu,req_qu) ! COM03
+       CALL wait_message(req_qu) ! COM03
 
        DO ind=1,ndomain
@@ -259,15 +280,19 @@
           ps=f_ps(ind)
           u=f_u(ind)
+          theta_rhodz = f_theta_rhodz(ind)
           mass=f_mass(ind)
           theta = f_theta(ind)
           qu=f_qu(ind)
+          qv=f_qv(ind)
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
-          CALL compute_geopot(ps,mass,theta, pk,geopot)
           hflux=f_hflux(ind)
           convm = f_dmass(ind)
           dtheta_rhodz=f_dtheta_rhodz(ind)
           du=f_du(ind)
+
           CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
+!          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .FALSE.) ! FIXME
+!          CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du) 
           IF(caldyn_eta==eta_mass) THEN
              wflux=f_wflux(ind)
@@ -278,9 +303,4 @@
        ENDDO       
    
-!       CALL checksum(f_geopot)
-!       CALL checksum(f_dmass)
-!       CALL checksum(f_pk)
-!       CALL checksum(f_pk)
-       
     CASE(enstrophy) ! enstrophy-conserving
        DO ind=1,ndomain
Index: codes/icosagcm/trunk/src/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 517)
+++ codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 519)
@@ -155,5 +155,5 @@
 
        IF(hydrostatic) THEN
-          CALL compute_caldyn_slow_hydro(u,mass,hflux,du)
+          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.)
        ELSE
           W = f_W(ind)
Index: codes/icosagcm/trunk/src/caldyn_kernels.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels.f90	(revision 517)
+++ codes/icosagcm/trunk/src/caldyn_kernels.f90	(revision 519)
@@ -32,5 +32,5 @@
   SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
     USE icosa
-    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
+    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop
     USE trace
     USE omp_para
@@ -48,18 +48,22 @@
     CALL trace_start("compute_pvort")  
 
-    IF(caldyn_eta==eta_mass) THEN
-       CALL wait_message(req_ps)  ! COM00
-    ELSE
-       CALL wait_message(req_mass) ! COM00
-    END IF
-    CALL wait_message(req_theta_rhodz) ! COM01
+!    IF(caldyn_eta==eta_mass) THEN
+!       CALL wait_message(req_ps)  ! COM00
+!    ELSE
+!       CALL wait_message(req_mass) ! COM00
+!    END IF
+!    CALL wait_message(req_theta_rhodz) ! COM01
 
     IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
        DO l = ll_begin,ll_end
-          CALL test_message(req_u) 
+!          CALL test_message(req_u) 
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext
-             m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
-             rhodz(ij,l) = m
+             IF(DEC) THEN  ! ps is actually Ms
+                m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l)
+             ELSE
+                m = mass_dak(l)+ps(ij)*mass_dbk(l)
+             END IF
+             rhodz(ij,l) = m/g
              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
           ENDDO
@@ -67,5 +71,5 @@
     ELSE ! Compute only theta
        DO l = ll_begin,ll_end
-          CALL test_message(req_u) 
+ !         CALL test_message(req_u) 
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext
@@ -75,5 +79,5 @@
     END IF
 
-    CALL wait_message(req_u) ! COM02
+ !   CALL wait_message(req_u) ! COM02
 
 !!! Compute shallow-water potential vorticity
@@ -81,24 +85,39 @@
        !DIR$ SIMD
        DO ij=ij_begin_ext,ij_end_ext
-          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
-               + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
-               - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
-
-          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
-               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
-               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
-
-          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
-
-          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
-               + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
-               - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
-
-          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
-               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
-               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
-
-          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
-
+          IF(DEC) THEN
+             etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)         &
+                                   + ne_left * u(ij+t_rup+u_left,l)  &
+                                   - ne_lup  * u(ij+u_lup,l) )                               
+
+             hv =   Riv2(ij,vup)          * rhodz(ij,l)           &
+                  + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
+                  + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+             qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
+
+             etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          &
+                                      + ne_right * u(ij+t_ldown+u_right,l)  &
+                                      - ne_rdown * u(ij+u_rdown,l) )
+             hv =   Riv2(ij,vdown)        * rhodz(ij,l)          &
+                  + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
+                  + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+             qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
+          ELSE
+             etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)        * de(ij+u_rup)         &
+                                   + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
+                                   - ne_lup  * u(ij+u_lup,l)        * de(ij+u_lup) )                               
+
+             hv =   Riv2(ij,vup)          * rhodz(ij,l)           &
+                  + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
+                  + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+             qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
+
+             etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
+                                      + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
+                                      - ne_rdown * u(ij+u_rdown,l)          * de(ij+u_rdown) )
+             hv =   Riv2(ij,vdown)        * rhodz(ij,l)          &
+                  + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
+                  + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+             qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
+          END IF
        ENDDO
 
@@ -115,5 +134,5 @@
   END SUBROUTINE compute_pvort
 
-  !************************* caldyn_horiz = caldyn_fast + caldyn_slow **********************
+  !************** caldyn_horiz = caldyn_fast + caldyn_slow + caldyn_coriolis ***************
 
   SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du)
@@ -139,4 +158,5 @@
     REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
+    REAL(rstd) :: uu_right, uu_lup, uu_ldown
 
     INTEGER :: i,j,ij,l
@@ -149,11 +169,19 @@
     DO l = ll_begin, ll_end
 !!!  Compute mass and theta fluxes
-       IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+!       IF (caldyn_conserv==energy) CALL test_message(req_qu) 
        !DIR$ SIMD
        DO ij=ij_begin_ext,ij_end_ext         
-          hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
-          hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
-          hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
-
+          uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
+          uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
+          uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
+          uu_right= uu_right*le_de(ij+u_right)
+          uu_lup  = uu_lup  *le_de(ij+u_lup)
+          uu_ldown= uu_ldown*le_de(ij+u_ldown)
+          hflux(ij+u_right,l)=uu_right
+          hflux(ij+u_lup,l)  =uu_lup
+          hflux(ij+u_ldown,l)=uu_ldown
+!          hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
+!          hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
+!          hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
           Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
           Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
@@ -189,5 +217,5 @@
     CASE(energy) ! energy-conserving TRiSK
 
-       CALL wait_message(req_qu) ! COM03
+!       CALL wait_message(req_qu) ! COM03
 
        DO l=ll_begin,ll_end
@@ -205,5 +233,5 @@
                   wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
                   wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
-             du(ij+u_right,l) = .5*uu/de(ij+u_right)
+             du(ij+u_right,l) = .5*uu
 
              uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
@@ -217,6 +245,5 @@
                   wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
                   wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
-             du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
-
+             du(ij+u_lup,l) = .5*uu
 
              uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
@@ -230,5 +257,5 @@
                   wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
                   wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
-             du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
+             du(ij+u_ldown,l) = .5*uu
 
           ENDDO
@@ -251,5 +278,5 @@
                   wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
                   wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
-             du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
+             du(ij+u_right,l) = qu(ij+u_right,l)*uu
 
 
@@ -264,5 +291,5 @@
                   wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
                   wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
-             du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
+             du(ij+u_lup,l) = qu(ij+u_lup,l)*uu
 
              uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
@@ -276,5 +303,5 @@
                   wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
                   wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
-             du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
+             du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu
 
           ENDDO
@@ -290,12 +317,11 @@
           !DIR$ SIMD
           DO ij=ij_begin,ij_end         
-
-             berni(ij,l) = pk(ij,l) + &
-                  1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
-                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
-                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
-                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
-                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
-                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+             berni(ij,l) = pk(ij,l) + 1/(4*Ai(ij))*( &
+                  le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
+                  le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
+                  le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
+                  le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
+                  le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
+                  le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
              ! from now on pk contains the vertically-averaged geopotential
              pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
@@ -308,12 +334,12 @@
           !DIR$ SIMD
           DO ij=ij_begin,ij_end         
-
              berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                  + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
-                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
-                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
-                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
-                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
-                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+                  + 1/(4*Ai(ij))*( &
+                  le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
+                  le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
+                  le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
+                  le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
+                  le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
+                  le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
           ENDDO
        ENDDO
@@ -326,5 +352,5 @@
        DO ij=ij_begin,ij_end         
 
-          du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
+          du(ij+u_right,l) = du(ij+u_right,l) + (             &
                0.5*(theta(ij,l)+theta(ij+t_right,l))                &
                *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
@@ -332,10 +358,10 @@
 
 
-          du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
+          du(ij+u_lup,l) = du(ij+u_lup,l) +  (                  &
                0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
                *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
                + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
 
-          du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
+          du(ij+u_ldown,l) = du(ij+u_ldown,l) + (             &
                0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
                *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
Index: codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 517)
+++ codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 519)
@@ -649,9 +649,10 @@
   END SUBROUTINE compute_caldyn_Coriolis
 
-  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du)
+  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero)
+    LOGICAL, INTENT(IN) :: zero
     REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
-    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
     
     REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function
@@ -690,10 +691,22 @@
        ENDDO
        ! Compute du=-grad(Bernoulli)
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end
+       IF(zero) THEN
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end
              du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
              du(ij+u_lup,l)   = ne_lup*(berni(ij)-berni(ij+t_lup))
              du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
-       END DO
+          END DO
+       ELSE
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end
+             du(ij+u_right,l) = du(ij+u_right,l) + &
+                  ne_right*(berni(ij)-berni(ij+t_right))
+             du(ij+u_lup,l)   = du(ij+u_lup,l) + &
+                  ne_lup*(berni(ij)-berni(ij+t_lup))
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) + &
+                  ne_ldown*(berni(ij)-berni(ij+t_ldown))
+          END DO
+       END IF
     END DO
 
Index: codes/icosagcm/trunk/src/explicit_scheme.f90
===================================================================
--- codes/icosagcm/trunk/src/explicit_scheme.f90	(revision 517)
+++ codes/icosagcm/trunk/src/explicit_scheme.f90	(revision 519)
@@ -26,8 +26,10 @@
     REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:)
     REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:,:) , theta_rhodzm1(:,:,:), theta_rhodzm2(:,:,:), dtheta_rhodz(:,:,:)
     REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
     LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
     INTEGER :: it,stage
+    
+    CALL legacy_to_DEC(f_ps, f_u)
     DO stage=1,nb_stage
        !       CALL checksum(f_ps)
@@ -53,4 +55,5 @@
        END SELECT
     END DO
+    CALL DEC_to_legacy(f_ps, f_u)
 
   CONTAINS
@@ -142,5 +145,5 @@
                   um1(ij+u_lup,l)=u(ij+u_lup,l)
                   um1(ij+u_ldown,l)=u(ij+u_ldown,l)
-                  theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
+                  theta_rhodzm1(ij,l,1)=theta_rhodz(ij,l,1)
                ENDDO
             ENDDO
@@ -153,5 +156,5 @@
                u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
                u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
-               theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
+               theta_rhodz(ij,l,1)=theta_rhodzm1(ij,l,1)+tau*dtheta_rhodz(ij,l,1)
             ENDDO
          ENDDO
@@ -188,8 +191,8 @@
 
          IF (stage==1) THEN
-            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
+            psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:)
             ps(:)=psm1(:)+tau*dps(:)
             u(:,:)=um1(:,:)+tau*du(:,:)
-            theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
+            theta_rhodz(:,:,:)=theta_rhodzm1(:,:,:)+tau*dtheta_rhodz(:,:,:)
 
          ELSE IF (stage==2) THEN
@@ -197,8 +200,8 @@
             ps(:)=psm1(:)+tau*dps(:)
             u(:,:)=um1(:,:)+tau*du(:,:)
-            theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
-
-            psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
-            psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
+            theta_rhodz(:,:,:)=theta_rhodzm1(:,:,:)+tau*dtheta_rhodz(:,:,:)
+
+            psm2(:)=psm1(:) ; theta_rhodzm2(:,:,:)=theta_rhodzm1(:,:,:) ; um2(:,:)=um1(:,:) 
+            psm1(:)=ps(:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:) ; um1(:,:)=u(:,:) 
 
          ELSE 
@@ -206,8 +209,8 @@
             ps(:)=psm2(:)+2*tau*dps(:)
             u(:,:)=um2(:,:)+2*tau*du(:,:)
-            theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)
-
-            psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
-            psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
+            theta_rhodz(:,:,:)=theta_rhodzm2(:,:,:)+2*tau*dtheta_rhodz(:,:,:)
+
+            psm2(:)=psm1(:) ; theta_rhodzm2(:,:,:)=theta_rhodzm1(:,:,:) ; um2(:,:)=um1(:,:) 
+            psm1(:)=ps(:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:) ; um1(:,:)=u(:,:) 
 
          ENDIF
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 517)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 519)
@@ -211,20 +211,20 @@
     IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
 
+    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,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 (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  ! FIXME : write_start undefined
+   
+    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)
     CALL SYSTEM_CLOCK(count_rate=rate_clock)
     !$OMP END MASTER   
-
-    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,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 (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)  ! FIXME : write_start undefined
-   
-    CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
 
     DO it=itau0+1,itau0+itaumax
