Index: codes/icosagcm/devel/src/diagnostics/check_conserve.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/check_conserve.f90	(revision 583)
+++ codes/icosagcm/devel/src/diagnostics/check_conserve.f90	(revision 584)
@@ -299,5 +299,5 @@
     REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv
     
-    REAL(rstd) :: ucenter(iim*jjm,3,llm)
+    REAL(rstd) :: ucenter(iim*jjm,llm,3)
     REAL(rstd) :: ulon(iim*jjm,llm)
     REAL(rstd) :: ulat(iim*jjm,llm)
Index: codes/icosagcm/devel/src/diagnostics/wind.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/wind.f90	(revision 583)
+++ codes/icosagcm/devel/src/diagnostics/wind.f90	(revision 584)
@@ -1,12 +1,17 @@
 MODULE wind_mod
-
-CONTAINS
-
-  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat)
+  USE omp_para
   USE icosa
   IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: compute_wind_centered, compute_wind_centered_lonlat_compound, &
+       compute_wind_centered_from_lonlat_compound, compute_wind_perp_from_lonlat_compound, &
+       un2ulonlat, ulonlat2un, compute_wind2d_perp_from_lonlat_compound
+       
+CONTAINS
+
+  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat)
     TYPE(t_field), POINTER :: f_u(:) ! IN  : normal velocity components on edges
-    TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
-    
+    TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons    
     REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
     INTEGER :: ind
@@ -25,6 +30,4 @@
 
   SUBROUTINE ulonlat2un(f_ulon, f_ulat,f_u)
-  USE icosa
-  IMPLICIT NONE
     TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! IN : velocity reconstructed at hexagons
     TYPE(t_field), POINTER :: f_u(:) ! OUT  : normal velocity components on edges
@@ -46,9 +49,6 @@
  
   SUBROUTINE compute_wind_centered(ue,ucenter)
-  USE icosa
-  USE omp_para
-  IMPLICIT NONE
   REAL(rstd) :: ue(3*iim*jjm,llm)
-  REAL(rstd) :: ucenter(iim*jjm,3,llm)
+  REAL(rstd) :: ucenter(iim*jjm,llm,3)
   INTEGER :: i,j,ij,l    
  
@@ -57,5 +57,7 @@
         DO i=ii_begin,ii_end
           ij=(j-1)*iim+i
-          ucenter(ij,:,l)=1/Ai(ij)*                                                                                                &
+          ! Perot reconstruction based on Gauss' theorem
+          ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
+          ucenter(ij,l,:)=1/Ai(ij)*                                                                                                &
                         ( ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:))  &
                          + ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:))          &
@@ -72,10 +74,6 @@
   
  SUBROUTINE compute_wind_on_edge(ue,uedge)
-  USE icosa
-  USE omp_para
-    
-  IMPLICIT NONE
   REAL(rstd) :: ue(3*iim*jjm,llm)
-  REAL(rstd) :: uedge(3*iim*jjm,3,llm)
+  REAL(rstd) :: uedge(3*iim*jjm,llm,3)
 
   REAL(rstd) :: ut(3*iim*jjm,llm)
@@ -88,7 +86,7 @@
         DO i=ii_begin,ii_end
           ij=(j-1)*iim+i
-          uedge(ij+u_right,:,l)=ue(ij+u_right,l)*ep_e(ij+u_right,:)*ne(ij,right) + ut(ij+u_right,l)*et_e(ij+u_right,:)*ne(ij,right) 
-          uedge(ij+u_lup,:,l)=ue(ij+u_lup,l)*ep_e(ij+u_lup,:)*ne(ij,lup) + ut(ij+u_lup,l)*et_e(ij+u_lup,:)*ne(ij,lup)
-          uedge(ij+u_ldown,:,l)=ue(ij+u_ldown,l)*ep_e(ij+u_ldown,:)*ne(ij,ldown) + ut(ij+u_ldown,l)*et_e(ij+u_ldown,:)*ne(ij,ldown)
+          uedge(ij+u_right,l,:)=ue(ij+u_right,l)*ep_e(ij+u_right,:)*ne(ij,right) + ut(ij+u_right,l)*et_e(ij+u_right,:)*ne(ij,right) 
+          uedge(ij+u_lup,l,:)=ue(ij+u_lup,l)*ep_e(ij+u_lup,:)*ne(ij,lup) + ut(ij+u_lup,l)*et_e(ij+u_lup,:)*ne(ij,lup)
+          uedge(ij+u_ldown,l,:)=ue(ij+u_ldown,l)*ep_e(ij+u_ldown,:)*ne(ij,ldown) + ut(ij+u_ldown,l)*et_e(ij+u_ldown,:)*ne(ij,ldown)
         ENDDO
       ENDDO
@@ -100,7 +98,4 @@
  
  SUBROUTINE compute_tangential_compound(ue,ut)
-  USE icosa  
-  USE omp_para
-  IMPLICIT NONE
   REAL(rstd) :: ue(3*iim*jjm,llm)
   REAL(rstd) :: ut(3*iim*jjm,llm)
@@ -155,15 +150,36 @@
  END SUBROUTINE compute_tangential_compound
  
- SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat)
-  USE icosa  
-  USE omp_para
-    
-  IMPLICIT NONE
-  REAL(rstd) :: u(3*iim*jjm,3,llm)
-  REAL(rstd) :: ulon(3*iim*jjm,3,llm)
-  REAL(rstd) :: ulat(3*iim*jjm,3,llm)
-
-  INTEGER :: i,j,ij,l     
-    
+! SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat)
+!  REAL(rstd) :: u(3*iim*jjm,3,llm)
+!  REAL(rstd) :: ulon(3*iim*jjm,3,llm)
+!  REAL(rstd) :: ulat(3*iim*jjm,3,llm)
+!
+!  INTEGER :: i,j,ij,l     
+!    
+!  
+!    DO l=ll_begin,ll_end
+!      DO j=jj_begin-1,jj_end+1
+!        DO i=ii_begin-1,ii_end+1
+!          ij=(j-1)*iim+i
+!          ulon(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elon_e(ij+u_right,:))*elon_e(ij+u_right,:) 
+!          ulon(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elon_e(ij+u_lup,:))*elon_e(ij+u_lup,:)
+!          ulon(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elon_e(ij+u_ldown,:))*elon_e(ij+u_ldown,:)
+!          
+!          ulat(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elat_e(ij+u_right,:))*elat_e(ij+u_right,:) 
+!          ulat(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elat_e(ij+u_lup,:))*elat_e(ij+u_lup,:) 
+!          ulat(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elat_e(ij+u_ldown,:))*elat_e(ij+u_ldown,:) 
+!          
+!        ENDDO
+!      ENDDO
+!    ENDDO
+! 
+! END SUBROUTINE compute_wind_lonlat_compound
+ 
+  SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u)
+  REAL(rstd) :: u(3*iim*jjm,llm,3)
+  REAL(rstd) :: ulon(3*iim*jjm,llm)
+  REAL(rstd) :: ulat(3*iim*jjm,llm)
+
+  INTEGER :: i,j,ij,l     
   
     DO l=ll_begin,ll_end
@@ -171,36 +187,7 @@
         DO i=ii_begin-1,ii_end+1
           ij=(j-1)*iim+i
-          ulon(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elon_e(ij+u_right,:))*elon_e(ij+u_right,:) 
-          ulon(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elon_e(ij+u_lup,:))*elon_e(ij+u_lup,:)
-          ulon(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elon_e(ij+u_ldown,:))*elon_e(ij+u_ldown,:)
-          
-          ulat(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elat_e(ij+u_right,:))*elat_e(ij+u_right,:) 
-          ulat(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elat_e(ij+u_lup,:))*elat_e(ij+u_lup,:) 
-          ulat(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elat_e(ij+u_ldown,:))*elat_e(ij+u_ldown,:) 
-          
-        ENDDO
-      ENDDO
-    ENDDO
- 
- END SUBROUTINE compute_wind_lonlat_compound
- 
-  SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u)
-  USE icosa  
-  USE omp_para
-    
-  IMPLICIT NONE
-  REAL(rstd) :: u(3*iim*jjm,3,llm)
-  REAL(rstd) :: ulon(3*iim*jjm,llm)
-  REAL(rstd) :: ulat(3*iim*jjm,llm)
-
-  INTEGER :: i,j,ij,l     
-  
-    DO l=ll_begin,ll_end
-      DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-          ij=(j-1)*iim+i
-          u(ij+u_right,:,l)=ulon(ij+u_right,l)*elon_e(ij+u_right,:)+ ulat(ij+u_right,l)*elat_e(ij+u_right,:)
-          u(ij+u_lup,:,l)=ulon(ij+u_lup,l)*elon_e(ij+u_lup,:) + ulat(ij+u_lup,l)*elat_e(ij+u_lup,:)
-          u(ij+u_ldown,:,l)=ulon(ij+u_ldown,l)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown,l)*elat_e(ij+u_ldown,:)
+          u(ij+u_right,l,:)=ulon(ij+u_right,l)*elon_e(ij+u_right,:)+ ulat(ij+u_right,l)*elat_e(ij+u_right,:)
+          u(ij+u_lup,l,:)=ulon(ij+u_lup,l)*elon_e(ij+u_lup,:) + ulat(ij+u_lup,l)*elat_e(ij+u_lup,:)
+          u(ij+u_ldown,l,:)=ulon(ij+u_ldown,l)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown,l)*elat_e(ij+u_ldown,:)
         ENDDO
       ENDDO
@@ -210,12 +197,7 @@
  
   SUBROUTINE compute_wind_centered_from_lonlat_compound(ulon, ulat, u)
-  USE icosa  
-  USE omp_para
-    
-  IMPLICIT NONE
-  REAL(rstd) :: u(iim*jjm,3,llm)
+  REAL(rstd) :: u(iim*jjm,llm,3)
   REAL(rstd) :: ulon(iim*jjm,llm)
   REAL(rstd) :: ulat(iim*jjm,llm)
-
   INTEGER :: i,j,ij,l     
   DO l=ll_begin,ll_end
@@ -223,15 +205,11 @@
         DO i=ii_begin-1,ii_end+1
           ij=(j-1)*iim+i
-          u(ij,:,l)=ulon(ij,l)*elon_i(ij,:)+ ulat(ij,l)*elat_i(ij,:)
-        ENDDO
-      ENDDO
-    ENDDO
- 
+          u(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ ulat(ij,l)*elat_i(ij,:)
+        ENDDO
+      ENDDO
+    ENDDO 
   END SUBROUTINE compute_wind_centered_from_lonlat_compound
  
   SUBROUTINE compute_wind2D_from_lonlat_compound(ulon, ulat, u)
-  USE icosa  
-  
-  IMPLICIT NONE
   REAL(rstd) :: u(3*iim*jjm,3)
   REAL(rstd) :: ulon(3*iim*jjm)
@@ -252,12 +230,8 @@
  
   SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up)
-  USE icosa  
-  USE omp_para
-    
-  IMPLICIT NONE
   REAL(rstd) :: up(3*iim*jjm,llm)
   REAL(rstd) :: ulon(3*iim*jjm,llm)
   REAL(rstd) :: ulat(3*iim*jjm,llm)
-  REAL(rstd) :: u(3*iim*jjm,3,llm)
+  REAL(rstd) :: u(3*iim*jjm,llm,3)
 
   INTEGER :: i,j,ij,l     
@@ -269,7 +243,7 @@
         DO i=ii_begin-1,ii_end+1
           ij=(j-1)*iim+i
-          up(ij+u_right,l)=sum(u(ij+u_right,:,l)*ep_e(ij+u_right,:))
-          up(ij+u_lup,l)=sum(u(ij+u_lup,:,l)*ep_e(ij+u_lup,:))
-          up(ij+u_ldown,l)=sum(u(ij+u_ldown,:,l)*ep_e(ij+u_ldown,:))
+          up(ij+u_right,l)=sum(u(ij+u_right,l,:)*ep_e(ij+u_right,:))
+          up(ij+u_lup,l)=sum(u(ij+u_lup,l,:)*ep_e(ij+u_lup,:))
+          up(ij+u_ldown,l)=sum(u(ij+u_ldown,l,:)*ep_e(ij+u_ldown,:))
         ENDDO
       ENDDO
@@ -279,7 +253,4 @@
    
   SUBROUTINE compute_wind2D_perp_from_lonlat_compound(ulon, ulat, up)
-  USE icosa  
-    
-  IMPLICIT NONE
   REAL(rstd) :: up(3*iim*jjm)
   REAL(rstd) :: ulon(3*iim*jjm)
@@ -302,21 +273,16 @@
    
  SUBROUTINE compute_wind_centered_lonlat_compound(uc, ulon, ulat)
-  USE icosa  
-  USE omp_para
-    
-  IMPLICIT NONE
-  REAL(rstd) :: uc(iim*jjm,3,llm)
+  REAL(rstd) :: uc(iim*jjm,llm,3)
   REAL(rstd) :: ulon(iim*jjm,llm)
   REAL(rstd) :: ulat(iim*jjm,llm)
 
   INTEGER :: i,j,ij,l     
-    
-  
-    DO l=ll_begin,ll_end
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          ulon(ij,l)=sum(uc(ij,:,l)*elon_i(ij,:))
-          ulat(ij,l)=sum(uc(ij,:,l)*elat_i(ij,:)) 
+  
+    DO l=ll_begin,ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          ulon(ij,l)=sum(uc(ij,l,:)*elon_i(ij,:))
+          ulat(ij,l)=sum(uc(ij,l,:)*elat_i(ij,:)) 
         ENDDO
       ENDDO
@@ -326,20 +292,16 @@
 
  SUBROUTINE compute_wind_centered_from_wind_lonlat_centered(ulon, ulat,uc)
-  USE icosa  
-  USE omp_para
-    
-  IMPLICIT NONE
   REAL(rstd) :: ulon(iim*jjm,llm)
   REAL(rstd) :: ulat(iim*jjm,llm)
-  REAL(rstd) :: uc(iim*jjm,3,llm)
-
-  INTEGER :: i,j,ij,l     
-    
-  
-    DO l=ll_begin,ll_end
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          uc(ij,:,l)=ulon(ij,l)*elon_i(ij,:)+ulat(ij,l)*elat_i(ij,:)
+  REAL(rstd) :: uc(iim*jjm,llm,3)
+
+  INTEGER :: i,j,ij,l     
+    
+  
+    DO l=ll_begin,ll_end
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          uc(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ulat(ij,l)*elat_i(ij,:)
         ENDDO
       ENDDO
@@ -348,12 +310,8 @@
  END SUBROUTINE compute_wind_centered_from_wind_lonlat_centered
 
-
-
  SUBROUTINE compute_wind_perp_from_wind_centered(uc,un)
-  USE icosa  
-  USE omp_para
-    
+
   IMPLICIT NONE
-  REAL(rstd),INTENT(IN)   :: uc(iim*jjm,3,llm)
+  REAL(rstd),INTENT(IN)   :: uc(iim*jjm,llm,3)
   REAL(rstd),INTENT(OUT)  :: un(3*iim*jjm,llm)
 
@@ -365,7 +323,7 @@
         DO i=ii_begin,ii_end
           ij=(j-1)*iim+i
-          un(ij+u_right,l) = sum(0.5*(uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:))
-          un(ij+u_lup,l) = sum(0.5*(uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:))
-          un(ij+u_ldown,l) = sum(0.5*(uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:))
+          un(ij+u_right,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_right,l,:))*ep_e(ij+u_right,:))
+          un(ij+u_lup,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_lup,l,:))*ep_e(ij+u_lup,:))
+          un(ij+u_ldown,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:))
          ENDDO
       ENDDO
@@ -376,12 +334,9 @@
 
  SUBROUTINE compute_un2ulonlat(un, ulon, ulat)
-  USE icosa  
-    
-  IMPLICIT NONE
   REAL(rstd),INTENT(IN)  :: un(3*iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm)
 
-  REAL(rstd)             :: uc(iim*jjm,3,llm)
+  REAL(rstd)             :: uc(iim*jjm,llm,3)
     
   CALL compute_wind_centered(un,uc)  
@@ -391,12 +346,9 @@
 
  SUBROUTINE compute_ulonlat2un(ulon, ulat,un)
-  USE icosa  
-    
-  IMPLICIT NONE
   REAL(rstd),INTENT(IN) :: ulon(iim*jjm,llm)
   REAL(rstd),INTENT(IN) :: ulat(iim*jjm,llm)
   REAL(rstd),INTENT(OUT)  :: un(3*iim*jjm,llm)
 
-  REAL(rstd)             :: uc(iim*jjm,3,llm)
+  REAL(rstd)             :: uc(iim*jjm,llm,3)
     
     CALL compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc)  
Index: codes/icosagcm/devel/src/physics/physics.f90
===================================================================
--- codes/icosagcm/devel/src/physics/physics.f90	(revision 583)
+++ codes/icosagcm/devel/src/physics/physics.f90	(revision 584)
@@ -247,5 +247,5 @@
 
     REAL(rstd) :: p(iim*jjm,llm+1)
-    REAL(rstd) :: uc(iim*jjm,3,llm)
+    REAL(rstd) :: uc(iim*jjm,llm,3)
     REAL(rstd) :: ulon(iim*jjm,llm)
     REAL(rstd) :: ulat(iim*jjm,llm)
@@ -311,5 +311,5 @@
     REAL(rstd) :: dulat(iim*jjm,llm)
     REAL(rstd) :: ue(3*iim*jjm,llm)
-    REAL(rstd) :: duc(iim*jjm,3,llm)
+    REAL(rstd) :: duc(iim*jjm,llm,3)
     REAL(rstd) :: dt2, due
     INTEGER :: i,j,ij,l
@@ -321,11 +321,11 @@
         DO i=ii_begin,ii_end
           ij=(j-1)*iim+i
-          due = sum( (duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
+          due = sum( (duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) )
           ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
 
-          due = sum( (duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
+          due = sum( (duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) )
           ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
 
-          due = sum( (duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
+          due = sum( (duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) )
           ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
         ENDDO
Index: codes/icosagcm/devel/src/physics/physics_lmdz_generic.F90
===================================================================
--- codes/icosagcm/devel/src/physics/physics_lmdz_generic.F90	(revision 583)
+++ codes/icosagcm/devel/src/physics/physics_lmdz_generic.F90	(revision 584)
@@ -78,5 +78,5 @@
     CALL allocate_field(f_dq,field_t,type_real,llm,nqtot)
     CALL allocate_field(f_dps,field_t,type_real)
-    CALL allocate_field(f_duc,field_t,type_real,3,llm)    
+    CALL allocate_field(f_duc,field_t,type_real,llm,3)
 !$OMP END PARALLEL    
 
@@ -362,5 +362,5 @@
           DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
-            duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:)
+            duc(ij,l,:)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:)
           ENDDO
         ENDDO
@@ -371,7 +371,7 @@
           DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
-            u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
-            u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
-            u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
+            u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) )
+            u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) )
+            u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) )
           ENDDO
         ENDDO
Index: codes/icosagcm/devel/src/time/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 583)
+++ codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 584)
@@ -168,4 +168,5 @@
     USE caldyn_mod
     USE advect_tracer_mod
+    USE diagflux_mod
     USE physics_mod
     USE mpipara
@@ -213,5 +214,5 @@
     CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)  
 
-    CALL trace_on
+    Call trace_on
 
     IF (xios_output) THEN ! we must call update_calendar before any XIOS output
@@ -298,5 +299,5 @@
        
        IF(MOD(it,itau_adv)==0) THEN
-          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
+          CALL advect_tracer(diagflux_on, f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz,f_qfluxt)  ! update q and rhodz after RK step
           fluxt_zero=.TRUE.
           ! FIXME : check that rhodz is consistent with ps
