Index: codes/icosagcm/trunk/src/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/advect.f90	(revision 137)
+++ codes/icosagcm/trunk/src/advect.f90	(revision 138)
@@ -3,4 +3,8 @@
   IMPLICIT NONE 
 
+  PRIVATE
+  
+  PUBLIC :: init_advect, compute_backward_traj, compute_gradq3d, compute_advect_horiz
+
 CONTAINS
 
@@ -8,9 +12,4 @@
 
   SUBROUTINE init_advect(normal,tangent)
-    USE domain_mod
-    USE dimensions
-    USE geometry
-    USE metric
-    USE vector
     IMPLICIT NONE
     REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
@@ -48,8 +47,4 @@
 
   SUBROUTINE compute_gradq3d(qi,gradq3d)
-    USE domain_mod
-    USE dimensions
-    USE geometry
-    USE metric
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm)
@@ -62,5 +57,10 @@
     REAL(rstd) :: ar
     INTEGER :: i,j,n,k,ind,l
+
+    ! TODO : precompute ar, drop arr as output argument of gradq ?
+
     !========================================================================================== GRADIENT 
+    ! Compute gradient at triangles solving a linear system
+    ! arr = area of triangle joining centroids of hexagons
     Do l = 1,llm 
        DO j=jj_begin-1,jj_end+1
@@ -114,10 +114,6 @@
   END SUBROUTINE compute_gradq3d
 
-  ! Backward trajectories, used in Miura approach
+  ! Backward trajectories, for use with Miura approach
   SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
-    USE domain_mod
-    USE dimensions
-    USE geometry
-    USE metric
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
@@ -127,9 +123,9 @@
     REAL(rstd),INTENT(IN)    :: tau
 
-    REAL(rstd) :: v_e(3), up_e
-    REAL(rstd) :: qe
-    REAL(rstd) :: ed(3)    
+    REAL(rstd) :: v_e(3), up_e, qe, ed(3)    
     INTEGER :: i,j,n,l
 
+    ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
+    
     ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
     DO l = 1,llm
@@ -189,174 +185,101 @@
   ! Horizontal transport (S. Dubey, T. Dubos)
   ! Slope-limited van Leer approach with hexagons
-  SUBROUTINE compute_advect_horiz(update_mass, normal,tangent,qi,gradq3d,him,ue,hfluxt,bigt)
-    USE domain_mod
-    USE dimensions
-    USE geometry
-    USE metric
-    IMPLICIT NONE
-    LOGICAL, INTENT(IN)      :: update_mass
-    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
-    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
-    REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
-    REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux
-    REAL(rstd),INTENT(IN)    :: bigt
-
-    REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm) 
-    REAL(rstd) :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
-    REAL(rstd) :: v_e(3*iim*jjm,llm,3)
-    REAL(rstd) :: up_e
-
-!    REAL(rstd) :: qe(3*iim*jjm,llm)
-    REAL(rstd) :: qe
-    REAL(rstd) :: ed(3)    
-    REAL(rstd) :: flxx(3*iim*jjm,llm)
+  SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 
+    IMPLICIT NONE
+    LOGICAL, INTENT(IN)       :: update_mass
+    REAL(rstd), INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
+    REAL(rstd), INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux
+    REAL(rstd), INTENT(IN)    :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
+    REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm)
+    REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm)
+
+    REAL(rstd) :: dq,dmass,qe,ed(3), newmass
+    REAL(rstd) :: qflux(3*iim*jjm,llm)
     INTEGER :: i,j,n,k,l
-    REAL(rstd):: him_old(llm) 
-
-
-    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*dt/2
-    ! TODO : this does not depend on q and should be done only once, not for each tracer
+
+    ! evaluate tracer field at point cc using piecewise linear reconstruction
+    ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
+    ! ne*hfluxt>0 iff outgoing
     DO l = 1,llm
        DO j=jj_begin-1,jj_end+1
           DO i=ii_begin-1,ii_end+1
              n=(j-1)*iim+i
-
-             up_e =1/de(n+u_right)*(                                                   &
-                  wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        &
-                  wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        &
-                  wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      &
-                  wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    &
-                  wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
-                  wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    &
-                  wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    &
-                  wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    &
-                  wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        &
-                  wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         &
-                  )
-             v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:)
-
-             up_e=1/de(n+u_lup)*(                                                        &
-                  wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        &
-                  wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      &
-                  wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      &
-                  wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      &
-                  wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          & 
-                  wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          &
-                  wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              &
-                  wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              &
-                  wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            &
-                  wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           &
-                  )
-             v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:)
-
-             up_e=1/de(n+u_ldown)*(                                                    &
-                  wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &
-                  wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    &
-                  wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        &
-                  wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        &
-                  wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
-                  wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        &
-                  wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      &
-                  wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    &
-                  wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    &
-                  wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     &
-                  )
-             v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:)
-
-             cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt    !! redge is mid point of edge i 
-             cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt 
-             cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt 
-          ENDDO
-       ENDDO
-    END DO
-    ! end of tracer-independant computations
-
-    ! evaluate traver field at point cc usin piecewise linear reconstruction
-    ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
-    DO l = 1,llm
-       DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-             n=(j-1)*iim+i
-             if (ne(n,right)*ue(n+u_right,l)>0) then 
+             if (ne(n,right)*hfluxt(n+u_right,l)>0) then 
                 ed = cc(n+u_right,l,:) - xyz_i(n,:)
-                qe=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
+                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 
              else
                 ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:)
-                qe=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 
+                qe = qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 
              endif
-             flxx(n+u_right,l) = hfluxt(n+u_right,:)*qe*ne(n,right)
+             qflux(n+u_right,l) = hfluxt(n+u_right,l)*qe
              
-             if (ne(n,lup)*ue(n+u_lup,l)>0) then 
+             if (ne(n,lup)*hfluxt(n+u_lup,l)>0) then 
                 ed = cc(n+u_lup,l,:) - xyz_i(n,:)
-                qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)
+                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)
              else
                 ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:)
-                qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 
+                qe = qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 
              endif
-             flxx(n+u_lup,:) = hfluxt(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup) 
+             qflux(n+u_lup,l) = hfluxt(n+u_lup,l)*qe 
              
-             if (ne(n,ldown)*ue(n+u_ldown,l)>0) then 
+             if (ne(n,ldown)*hfluxt(n+u_ldown,l)>0) then 
                 ed = cc(n+u_ldown,l,:) - xyz_i(n,:)
-                qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed) 
+                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 
              else
                 ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:)
-                qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 
+                qe = qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 
              endif
-             flxx(n+u_ldown,:) = hfluxt(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown)
+             qflux(n+u_ldown,l) = hfluxt(n+u_ldown,l)*qe
           end do
        end do
     END DO
 
-    ! evaluate q flux as mass flux * qe
-    ! TODO : loop over k !
-    ! TODO : merge with previous loop and eliminate unnecessary temporary qe
-    DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-          n=(j-1)*iim+i 
-       ENDDO
-    ENDDO
-
-    ! update q and, if update_mass, mass
-    ! TODO : loop over k !
-    DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-          n=(j-1)*iim+i  
-
-          dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) &
-               + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) &
-               + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) ) 
-
-          dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       & 
-               +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    &
-               - flxx(n+u_left,:) - flxx(n+u_rdown,:) )   
-          him_old(:) = him(n,:) 
-          him(n,:) = him(n,:) + dhi(n,:) 
-          qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:)  
-
+    ! update q and, if update_mass, update mass
+    DO l = 1,llm
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             n=(j-1)*iim+i  
+             ! sign convention as Ringler et al. (2010) eq. 21
+             dmass =   hfluxt(n+u_right,l)*ne(n,right)   &
+                     + hfluxt(n+u_lup,l)  *ne(n,lup)     &
+                     + hfluxt(n+u_ldown,l)*ne(n,ldown)   &
+                     + hfluxt(n+u_rup,l)  *ne(n,rup)     &
+                     + hfluxt(n+u_left,l) *ne(n,left)    &
+                     + hfluxt(n+u_rdown,l)*ne(n,rdown)
+
+             dq    =   qflux(n+u_right,l) *ne(n,right)   &
+                     + qflux(n+u_lup,l)   *ne(n,lup)     &
+                     + qflux(n+u_ldown,l) *ne(n,ldown)   &
+                     + qflux(n+u_rup,l)   *ne(n,rup)     &
+                     + qflux(n+u_left,l)  *ne(n,left)    &
+                     + qflux(n+u_rdown,l) *ne(n,rdown)
+
+             newmass = mass(n,l) - dmass/Ai(n)
+             qi(n,l) = (qi(n,l)*mass(n,l) - dq/Ai(n) ) / newmass
+             IF(update_mass) mass(n,l) = newmass
+          END DO
        END DO
     END DO
 
   CONTAINS
+
     !====================================================================================
-    REAL FUNCTION sum2(a1,a2)
+    PURE REAL(rstd) FUNCTION sum2(a1,a2)
       IMPLICIT NONE
-      REAL,INTENT(IN):: a1(3), a2(3)
+      REAL(rstd),INTENT(IN):: a1(3), a2(3)
       sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
+      !      sum2 = 0. ! Godunov scheme
     END FUNCTION sum2
 
   END SUBROUTINE compute_advect_horiz
+
   !==========================================================================
-  SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det)
-    USE geometry
-    USE metric
-    USE dimensions
+  PURE SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det)
     IMPLICIT NONE
     INTEGER, INTENT(IN) :: n0,n1,n2,n3
-    REAL,INTENT(IN)     :: q(iim*jjm)
-    REAL,INTENT(OUT)    :: dq(3)
-    REAL(rstd)  ::det,detx,dety,detz
+    REAL(rstd), INTENT(IN)     :: q(iim*jjm)
+    REAL(rstd), INTENT(OUT)    :: dq(3), det
+    
+    REAL(rstd)  ::detx,dety,detz
     INTEGER    :: info
     INTEGER    :: IPIV(3)
@@ -365,8 +288,8 @@
     REAL(rstd) :: B(3)
 
-    A(1,1)=xyz_i(n1,1) -xyz_i(n0,1); A(1,2)=xyz_i(n1,2)- xyz_i(n0,2); A(1,3)=xyz_i(n1,3)  - xyz_i(n0,3) 
-    A(2,1)=xyz_i(n2,1) - xyz_i(n0,1); A(2,2)=xyz_i(n2,2) - xyz_i(n0,2); A(2,3)=xyz_i(n2,3) - xyz_i(n0,3) 
-    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);             A(3,3)= xyz_v(n3,3)
-
+    ! TODO : replace A by A1,A2,A3
+    A(1,1)=xyz_i(n1,1)-xyz_i(n0,1);  A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3) 
+    A(2,1)=xyz_i(n2,1)-xyz_i(n0,1);  A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3) 
+    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);            A(3,3)=xyz_v(n3,3)
 
     dq(1) = q(n1)-q(n0)
@@ -382,9 +305,11 @@
     dq(3) = detz 
   END SUBROUTINE gradq
+
   !==========================================================================
-  SUBROUTINE determinant(a1,a2,a3,det)
-    IMPLICIT NONE
-    REAL, DIMENSION(3) :: a1, a2,a3
-    REAL ::  det,x1,x2,x3
+  PURE SUBROUTINE determinant(a1,a2,a3,det)
+    IMPLICIT NONE
+    REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3
+    REAL(rstd), INTENT(OUT) :: det
+    REAL(rstd) ::  x1,x2,x3
     x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
     x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
Index: codes/icosagcm/trunk/src/advect_tracer.f90
===================================================================
--- codes/icosagcm/trunk/src/advect_tracer.f90	(revision 137)
+++ codes/icosagcm/trunk/src/advect_tracer.f90	(revision 138)
@@ -1,4 +1,5 @@
 MODULE advect_tracer_mod
   USE icosa
+  IMPLICIT NONE
   PRIVATE
 
@@ -8,7 +9,4 @@
   TYPE(t_field),POINTER :: f_cc(:)  ! starting point of backward-trajectory (Miura approach)
 
-!  TYPE(t_field),POINTER :: f_rhodzm1(:)
-!  TYPE(t_field),POINTER :: f_rhodz(:)
-
   REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz
 
@@ -19,16 +17,12 @@
   SUBROUTINE init_advect_tracer
     USE advect_mod
-    IMPLICIT NONE
     REAL(rstd),POINTER :: tangent(:,:)
     REAL(rstd),POINTER :: normal(:,:)
     INTEGER :: ind
 
-    CALL allocate_field(f_normal,field_u,type_real,3)
-    CALL allocate_field(f_tangent,field_u,type_real,3)
-    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3)
-    CALL allocate_field(f_cc,field_u,type_real,llm,3)
-
-!    CALL allocate_field(f_rhodz,field_t,type_real,llm) 
-!    CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 
+    CALL allocate_field(f_normal,field_u,type_real,3, name='normal')
+    CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent')
+    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d')
+    CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc')
 
     DO ind=1,ndomain
@@ -43,14 +37,5 @@
 
   SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz)
-    USE icosa
-    USE field_mod
-    USE domain_mod
-    USE dimensions
-    USE grid_param
-    USE geometry
-    USE metric
     USE advect_mod
-    USE advect_mod
-    USE disvert_mod
     USE mpipara
     IMPLICIT NONE
@@ -62,15 +47,97 @@
     TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step
 
-    REAL(rstd),POINTER :: q(:,:,:), normal(:,:,:), tangent(:,:,:), gradq3d(:,:,:)
+    REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), gradq3d(:,:,:), cc(:,:,:)
     REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:)
     REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 
-    INTEGER :: ind
+    INTEGER :: ind,k
 
     CALL transfert_request(f_u,req_e1)
-    CALL transfert_request(f_u,req_e1)
+!    CALL transfert_request(f_hfluxt,req_e1)  ! BUG : This (unnecessary) transfer makes the computation go wrong
+    CALL transfert_request(f_wfluxt,req_i1)
     CALL transfert_request(f_q,req_i1)
+    CALL transfert_request(f_rhodz,req_i1)
 	
     IF (is_mpi_root) PRINT *, 'Advection scheme'
 
+!    DO ind=1,ndomain
+!       CALL swap_dimensions(ind)
+!       CALL swap_geometry(ind)
+!       normal  = f_normal(ind)
+!       tangent = f_tangent(ind)
+!       cc      = f_cc(ind)
+!       u       = f_u(ind)
+!       q       = f_q(ind)
+!       rhodz   = f_rhodz(ind)
+!       hfluxt  = f_hfluxt(ind) 
+!       wfluxt  = f_wfluxt(ind) 
+!       gradq3d = f_gradq3d(ind)
+!
+!       ! 1/2 vertical transport
+!       DO k = 1, nqtot
+!          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k))
+!       END DO
+!
+!       ! horizontal transport
+!       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
+!       DO k = 1,nqtot
+!          CALL compute_gradq3d(q(:,:,k),gradq3d)
+!          CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k))
+!       END DO
+!
+!       ! 1/2 vertical transport
+!       DO k = 1,nqtot
+!          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))
+!       END DO
+!    END DO
+
+    ! 1/2 vertical transport + back-trajectories
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       normal  = f_normal(ind)
+       tangent = f_tangent(ind)
+       cc      = f_cc(ind)
+       u       = f_u(ind)
+       q       = f_q(ind)
+       rhodz   = f_rhodz(ind)
+       wfluxt  = f_wfluxt(ind) 
+       DO k = 1, nqtot
+          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k))
+       END DO
+       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
+    END DO
+
+    CALL transfert_request(f_q,req_i1)      ! necessary ?
+    CALL transfert_request(f_rhodz,req_i1)  ! necessary ?
+
+    ! horizontal transport - split in two to place transfer of gradq3d
+
+    DO k = 1, nqtot
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          q       = f_q(ind)
+          gradq3d = f_gradq3d(ind)
+          CALL compute_gradq3d(q(:,:,k),gradq3d)
+       END DO
+
+       CALL transfert_request(f_gradq3d,req_i1)
+
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          cc      = f_cc(ind)
+          q       = f_q(ind)
+          rhodz   = f_rhodz(ind)
+          hfluxt  = f_hfluxt(ind) 
+          gradq3d = f_gradq3d(ind)
+          CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k))
+       END DO
+    END DO 
+
+    CALL transfert_request(f_q,req_i1)      ! necessary ?
+    CALL transfert_request(f_rhodz,req_i1)  ! necessary ?
+
+    ! 1/2 vertical transport
     DO ind=1,ndomain
        CALL swap_dimensions(ind)
@@ -78,52 +145,12 @@
        q       = f_q(ind)
        rhodz   = f_rhodz(ind)
-       hfluxt  = f_hfluxt(ind) 
        wfluxt  = f_wfluxt(ind) 
-       normal  = f_normal(ind)
-       tangent = f_tangent(ind)
-       gradq3d = f_gradq3d(ind)
-       CALL compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d)
+       DO k = 1,nqtot
+          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))
+       END DO
     END DO
 
   END SUBROUTINE advect_tracer
 
-  SUBROUTINE compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d)
-    USE icosa
-    USE field_mod
-    USE domain_mod
-    USE dimensions
-    USE grid_param
-    USE geometry
-    USE metric
-    USE advect_mod
-    USE advect_mod
-    USE disvert_mod
-    USE mpipara
-    REAL(rstd), INTENT(IN) :: tangent(:,:,:), normal(:,:,:)
-    REAL(rstd), INTENT(IN) :: hfluxt(:,:), wfluxt(:,:), u(:,:)
-    REAL(rstd), INTENT(INOUT) :: rhodz(:,:)
-    REAL(rstd), INTENT(INOUT) :: q(:,:,:)
-    REAL(rstd), INTENT(OUT)   :: gradq3d(:,:,:)
-    INTEGER :: k
-    ! for mass/tracer consistency each transport step also updates the mass field rhodz
-    ! mass is updated after all tracers have been updated, when k==nqtot
-
-    ! 1/2 vertical transport
-    DO k = 1, nqtot
-       CALL vlz(k==nqtot, 0.5,wfluxt,rhodz,q(:,:,k))
-    END DO
-
-    ! horizontal transport
-    DO k = 1,nqtot
-       CALL compute_gradq3d(q(:,:,k),gradq3d)
-       CALL compute_advect_horiz(k==nqtot, tangent,normal,q(:,:,k),gradq3d,rhodz,u,hfluxt,dt*itau_adv) 
-    END DO
-
-    ! 1/2 vertical transport
-    DO k = 1,nqtot
-       CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))
-    END DO
-  END SUBROUTINE compute_adv_tracer
-  
   SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q)
     !
@@ -135,5 +162,4 @@
     !     wfluxt >0 for upward transport
     !    ********************************************************************
-    USE icosa
     IMPLICIT NONE
     LOGICAL, INTENT(IN)       :: update_mass
@@ -152,6 +178,6 @@
     ! finite difference of q
     DO l=2,llm
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
              ij=(j-1)*iim+i
              dzqw(ij,l)=q(ij,l)-q(ij,l-1)
@@ -164,6 +190,6 @@
     ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l
     DO l=2,llm-1
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
              ij=(j-1)*iim+i
              IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
@@ -179,6 +205,6 @@
 
     ! 0 slope in top and bottom layers
-    DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
+    DO j=jj_begin-1,jj_end+1
+       DO i=ii_begin-1,ii_end+1
           ij=(j-1)*iim+i
           dzq(ij,1)=0.
@@ -190,28 +216,22 @@
     ! then amount of q leaving level l/l+1 = wq = w * qq
     DO l = 1,llm-1
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             IF(wfluxt(ij,l+1).gt.0.) THEN
-                ! downward transport, sigw<0
-                ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                
-                w          = fac*wfluxt(ij,l+1)
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             w = fac*wfluxt(ij,l+1)
+             IF(w>0.) THEN  ! upward transport, upwind side is at level l 
+                sigw       = w/mass(ij,l)
+                qq         = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
+             ELSE           ! downward transport, upwind side is at level l+1
                 sigw       = w/mass(ij,l+1)
-                qq         = q(ij,l+1) - 0.5*(1.+sigw)*dzq(ij,l+1)
-                wq(ij,l+1) = w*qq
-             ELSE
-                ! upward transport, sigw>0
-                ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
-                w          = fac*wfluxt(ij,l)
-                sigw       = w/mass(ij,l)
-                qq         = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l)
-                wq(ij,l+1) = w*qq
+                qq         = q(ij,l+1)-0.5*(1.+sigw)*dzq(ij,l+1) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                
              ENDIF
+             wq(ij,l+1) = w*qq
           ENDDO
        ENDDO
     END DO
     ! wq = 0 at top and bottom
-    DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
+    DO j=jj_begin-1,jj_end+1
+       DO i=ii_begin-1,ii_end+1
           ij=(j-1)*iim+i
           wq(ij,llm+1)=0.
@@ -222,6 +242,6 @@
     ! update q, mass is updated only after all q's have been updated
     DO l=1,llm
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
              ij=(j-1)*iim+i
              newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1))
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 137)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 138)
@@ -297,8 +297,8 @@
     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
 
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s
     REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
-    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux
+    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
     
     INTEGER :: i,j,ij,l
@@ -312,5 +312,5 @@
     REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
     REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux
-    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence
+    REAL(rstd),ALLOCATABLE,SAVE :: divm(:,:)  ! mass flux divergence
     REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
     
@@ -328,5 +328,5 @@
     ALLOCATE(phi(iim*jjm,llm))      ! geopotential
     ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
-    ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
+    ALLOCATE(divm(iim*jjm,llm))    ! mass flux convergence
     ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
 
@@ -363,6 +363,6 @@
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i  
-        ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
-        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  &
+        ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
+        divm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  &
                                 ne(ij,rup)*hflux(ij+u_rup,l)       +  &  
                                 ne(ij,lup)*hflux(ij+u_lup,l)       +  &  
@@ -383,5 +383,5 @@
   ENDDO
    
-!!! cumulate mass flux convergence from top to bottom
+!!! cumulate mass flux divergence from top to bottom
   DO  l = llm-1, 1, -1
 !$OMP DO
@@ -389,5 +389,5 @@
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
-        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+        divm(ij,l) = divm(ij,l) + divm(ij,l+1)
       ENDDO
     ENDDO
@@ -402,10 +402,10 @@
         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
         ! => w>0 for upward transport
-        wflux( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
+        wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 )
       ENDDO
     ENDDO
   ENDDO
 
-! compute dps, vertical mass flux at the surface = 0
+! compute dps, set vertical mass flux at the surface to 0
 !$OMP DO
   DO j=jj_begin,jj_end
@@ -414,5 +414,5 @@
       wflux(ij,1)  = 0.
       ! dps/dt = -int(div flux)dz
-      dps(ij)=-convm(ij,1) * g 
+      dps(ij)=-divm(ij,1) * g 
     ENDDO
   ENDDO
@@ -644,5 +644,5 @@
     DEALLOCATE(phi)     ! geopotential
     DEALLOCATE(Ftheta)  ! theta flux   
-    DEALLOCATE(convm)   ! mass flux convergence
+    DEALLOCATE(divm)    ! mass flux divergence
     DEALLOCATE(berni)   ! bernouilli term
 !!$OMP END MASTER
Index: codes/icosagcm/trunk/src/etat0_dcmip4.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip4.f90	(revision 137)
+++ codes/icosagcm/trunk/src/etat0_dcmip4.f90	(revision 138)
@@ -174,4 +174,5 @@
     
       q(:,:,1)=theta(:,:)
+      IF(nqtot>2) q(:,:,3)=1.
 
       DO l=1,llm
@@ -200,4 +201,5 @@
                                                             -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2)))
              q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta))
+             IF(nqtot>3) q(ij,l,4)=cos(lon)*cos(lat)
            ENDDO
          ENDDO
Index: codes/icosagcm/trunk/src/field.f90
===================================================================
--- codes/icosagcm/trunk/src/field.f90	(revision 137)
+++ codes/icosagcm/trunk/src/field.f90	(revision 138)
@@ -11,4 +11,5 @@
     
   TYPE t_field
+    CHARACTER(30)      :: name
     REAL(rstd),POINTER :: rval2d(:)
     REAL(rstd),POINTER :: rval3d(:,:)
@@ -45,5 +46,5 @@
 CONTAINS
 
-  SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2)
+  SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2,name)
   USE domain_mod
   IMPLICIT NONE
@@ -52,11 +53,17 @@
     INTEGER,INTENT(IN) :: data_type
     INTEGER,OPTIONAL :: dim1,dim2
+    CHARACTER(*), OPTIONAL :: name
     INTEGER :: ind
     INTEGER :: ii_size,jj_size
 
-    ALLOCATE(field(ndomain))    
+    ALLOCATE(field(ndomain))
 
     DO ind=1,ndomain
-  
+      IF(PRESENT(name)) THEN
+         field(ind)%name = name
+      ELSE
+         field(ind)%name = '(unkown)'
+      END IF
+
       IF (PRESENT(dim2)) THEN
         field(ind)%ndim=4 
@@ -224,5 +231,8 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=2 .OR. field%data_type/=type_real) STOP 'get_val_r2d : bad pointer assignation'        
+    IF (field%ndim/=2 .OR. field%data_type/=type_real) THEN
+       PRINT *, 'get_val_r2d : bad pointer assignation with ' // TRIM(field%name) 
+       STOP
+    END IF
     field_pt=>field%rval2d
   END SUBROUTINE  getval_r2d
@@ -233,5 +243,8 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=3 .OR. field%data_type/=type_real) STOP 'get_val_r3d : bad pointer assignation'        
+    IF (field%ndim/=3 .OR. field%data_type/=type_real) THEN
+       PRINT *, 'get_val_r3d : bad pointer assignation with ' // TRIM(field%name) 
+       STOP
+    END IF
     field_pt=>field%rval3d
   END SUBROUTINE  getval_r3d
@@ -242,8 +255,10 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=4 .OR. field%data_type/=type_real) STOP 'get_val_r4d : bad pointer assignation'        
+    IF (field%ndim/=4 .OR. field%data_type/=type_real) THEN
+       PRINT *, 'get_val_r4d : bad pointer assignation with ' // TRIM(field%name)
+       STOP
+    END IF
     field_pt=>field%rval4d
-  END SUBROUTINE  getval_r4d
-  
+  END SUBROUTINE  getval_r4d  
 
   
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 137)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 138)
@@ -46,5 +46,5 @@
   INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
   CHARACTER(len=255) :: scheme_name
-  LOGICAL :: fluxt_zero ! set to .TRUE. to start accumulating fluxes in time
+  LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
 !  INTEGER :: itaumax
 !  REAL(rstd) ::write_period
@@ -141,8 +141,18 @@
      CALL swap_geometry(ind)
      rhodz=f_rhodz(ind); ps=f_ps(ind)
-     CALL compute_rhodz(ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
+     CALL compute_rhodz(.TRUE., ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
   END DO
-  fluxt_zero=.FALSE.
-
+  fluxt_zero=.TRUE.
+
+  ! check that rhodz is consistent with ps
+  CALL transfert_request(f_rhodz,req_i1)
+  CALL transfert_request(f_ps,req_i1)
+  DO ind=1,ndomain
+     CALL swap_dimensions(ind)
+     CALL swap_geometry(ind)
+     rhodz=f_rhodz(ind); ps=f_ps(ind)
+     CALL compute_rhodz(.FALSE., ps, rhodz)   
+  END DO
+  
   DO it=0,itaumax
 
@@ -182,6 +192,23 @@
 
     IF(MOD(it+1,itau_adv)==0) THEN
+       CALL transfert_request(f_wfluxt,req_i1) ! FIXME
+!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME
+
        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
        fluxt_zero=.TRUE.
+
+       ! FIXME : check that rhodz is consistent with ps
+       CALL transfert_request(f_rhodz,req_i1)
+       CALL transfert_request(f_ps,req_i1)
+       CALL transfert_request(f_dps,req_i1) ! FIXME
+       CALL transfert_request(f_wflux,req_i1) ! FIXME
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); 
+          wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
+          CALL compute_rhodz(.FALSE., ps, rhodz)   
+       END DO
+
     END IF
 
@@ -195,12 +222,13 @@
     LOGICAL :: with_dps
     INTEGER :: ind
-      
     DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
        IF(with_dps) THEN
           ps=f_ps(ind) ; dps=f_dps(ind) ; 
           ps(:)=ps(:)+dt*dps(:)
           hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
-          wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
-          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero)
+          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
+          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
        END IF
        u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
@@ -221,4 +249,6 @@
       
       DO ind=1,ndomain
+         CALL swap_dimensions(ind)
+         CALL swap_geometry(ind)
          ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
          psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
@@ -234,6 +264,6 @@
          IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
             hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
-            wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
-            CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero)
+            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
+            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
          END IF
       END DO
@@ -245,4 +275,6 @@
 
       DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
@@ -277,4 +309,7 @@
     tau = dt/nb_stage
       DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+
         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
@@ -320,4 +355,6 @@
 
       DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
@@ -343,21 +380,46 @@
   END SUBROUTINE timeloop    
 
-  SUBROUTINE compute_rhodz(ps, rhodz)
+  SUBROUTINE compute_rhodz(comp, ps, rhodz)
     USE icosa
     USE disvert_mod
+    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
     REAL(rstd), INTENT(IN) :: ps(iim*jjm)
     REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm)
-    INTEGER :: l,i,j,ij
+    REAL(rstd) :: m, err
+    INTEGER :: l,i,j,ij,dd
+    err=0.
+    IF(comp) THEN 
+       dd=1
+    ELSE
+!       dd=-1
+       dd=0
+    END IF
+
     DO l = 1, llm
-       DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
+       DO j=jj_begin-dd,jj_end+dd
+          DO i=ii_begin-dd,ii_end+dd
              ij=(j-1)*iim+i
-             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
+             IF(comp) THEN
+                rhodz(ij,l) = m
+             ELSE
+                err = MAX(err,abs(m-rhodz(ij,l)))
+             END IF
           ENDDO
        ENDDO
     ENDDO
+
+    IF(.NOT. comp) THEN
+       IF(err>1e-10) THEN
+          PRINT *, 'Discrepancy between ps and rhodz detected', err
+          STOP
+       ELSE
+!          PRINT *, 'No discrepancy between ps and rhodz detected'
+       END IF
+    END IF
+
   END SUBROUTINE compute_rhodz
 
-  SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau,fluxt_zero)
+  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
     USE icosa
     REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
@@ -366,8 +428,10 @@
     LOGICAL, INTENT(INOUT) :: fluxt_zero
     IF(fluxt_zero) THEN
+!       PRINT *, 'Accumulating fluxes (first)'
        fluxt_zero=.FALSE.
        hfluxt = tau*hflux
        wfluxt = tau*wflux
     ELSE
+!       PRINT *, 'Accumulating fluxes (next)'
        hfluxt = hfluxt + tau*hflux
        wfluxt = wfluxt + tau*wflux
