Index: codes/icosagcm/trunk/src/caldyn.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn.f90	(revision 133)
+++ codes/icosagcm/trunk/src/caldyn.f90	(revision 134)
@@ -20,6 +20,6 @@
       CASE('gcm')
         CALL init_caldyn_gcm
-      CASE('adv')
-        CALL init_caldyn_adv
+!      CASE('adv')
+!        CALL init_caldyn_adv
       CASE DEFAULT
         PRINT*,'Bad selector for variable caldyn : <', TRIM(caldyn_type),'> options are <gcm>, <adv>'
@@ -30,5 +30,6 @@
   END SUBROUTINE init_caldyn
   
-  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
+  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
+       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
   USE icosa
   USE caldyn_gcm_mod, ONLY : caldyn_gcm=>caldyn
@@ -41,4 +42,6 @@
   TYPE(t_field),POINTER :: f_u(:)
   TYPE(t_field),POINTER :: f_q(:)
+  TYPE(t_field),POINTER :: f_hflux(:)
+  TYPE(t_field),POINTER :: f_wflux(:)
   TYPE(t_field),POINTER :: f_dps(:)
   TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
@@ -47,7 +50,8 @@
     SELECT CASE (TRIM(caldyn_type))
       CASE('gcm')
-        CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
-      CASE('adv')
-        CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
+        CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
+             f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
+!      CASE('adv')
+!        CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
     END SELECT
   
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 133)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 134)
@@ -35,5 +35,5 @@
        STOP
     END SELECT
-     IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
+    IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
 
     def='direct'
@@ -84,5 +84,6 @@
   END SUBROUTINE allocate_caldyn
    
-  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
+  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
+       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
     USE icosa
     USE vorticity_mod
@@ -97,4 +98,5 @@
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
     TYPE(t_field),POINTER :: f_dps(:)
     TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
@@ -103,5 +105,5 @@
     REAL(rstd),POINTER :: phis(:), ps(:), dps(:)
     REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:), du(:,:)
+    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
     REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
     INTEGER :: ind,ij
@@ -133,4 +135,6 @@
           CALL swap_geometry(ind)
           phis=f_phis(ind)
+          hflux=f_hflux(ind)
+          wflux=f_wflux(ind)
           ps=f_ps(ind)
           dps=f_dps(ind)
@@ -144,5 +148,6 @@
           out_u=f_out_u(ind) 
           !$OMP PARALLEL DEFAULT(SHARED)
-          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
+               hflux, wflux, dps, dtheta_rhodz, du)
           !$OMP END PARALLEL
        ENDDO       
@@ -155,4 +160,6 @@
           ps=f_ps(ind)
           dps=f_dps(ind)
+          hflux=f_hflux(ind)
+          wflux=f_wflux(ind)
           theta_rhodz=f_theta_rhodz(ind)
           dtheta_rhodz=f_dtheta_rhodz(ind)
@@ -165,5 +172,6 @@
           !$OMP PARALLEL DEFAULT(SHARED)
           CALL compute_pvort(ps, u, p,rhodz,qu)
-          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
+               hflux, wflux, dps, dtheta_rhodz, du)
           !$OMP END PARALLEL
        ENDDO
@@ -276,5 +284,5 @@
   END SUBROUTINE compute_pvort
   
-  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du)
   USE icosa
   USE disvert_mod
@@ -289,7 +297,8 @@
     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
 
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm)
     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
     
     INTEGER :: i,j,ij,l
@@ -301,8 +310,7 @@
     REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function
     REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:)
-    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)              ! geopotential
-    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:), Ftheta(:,:)  ! mass flux, theta flux
+    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 :: w(:,:)      ! vertical mass flux
     REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
     
@@ -319,8 +327,6 @@
     ALLOCATE(beta(iim*jjm,llm))
     ALLOCATE(phi(iim*jjm,llm))      ! geopotential
-    ALLOCATE(Fe(3*iim*jjm,llm))     ! mass flux   
     ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
     ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
-    ALLOCATE(w(iim*jjm,llm))        ! vertical velocity      
     ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
 
@@ -342,11 +348,11 @@
       DO i=ii_begin-1,ii_end+1
         ij=(j-1)*iim+i
-        Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
-        Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
-        Fe(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))*Fe(ij+u_right,l)
-        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l)
-        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l)
+        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)
+        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
       ENDDO
     ENDDO
@@ -358,10 +364,10 @@
         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)*Fe(ij+u_right,l)   +  &
-                                ne(ij,rup)*Fe(ij+u_rup,l)       +  &  
-                                ne(ij,lup)*Fe(ij+u_lup,l)       +  &  
-                                ne(ij,left)*Fe(ij+u_left,l)     +  &  
-                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  &  
-                                ne(ij,rdown)*Fe(ij+u_rdown,l))    
+        convm(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)       +  &  
+                                ne(ij,left)*hflux(ij+u_left,l)     +  &  
+                                ne(ij,ldown)*hflux(ij+u_ldown,l)   +  &  
+                                ne(ij,rdown)*hflux(ij+u_rdown,l))    
 
         ! signe ? attention d (rho theta dz)
@@ -396,5 +402,5 @@
         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
         ! => w>0 for upward transport
-        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
+        wflux( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
       ENDDO
     ENDDO
@@ -406,5 +412,5 @@
     DO i=ii_begin,ii_end
       ij=(j-1)*iim+i
-      w(ij,1)  = 0.
+      wflux(ij,1)  = 0.
       ! dps/dt = -int(div flux)dz
       dps(ij)=-convm(ij,1) * g 
@@ -423,39 +429,39 @@
                 ij=(j-1)*iim+i
 
-                uu = wee(ij+u_right,1,1)*Fe(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
-                     wee(ij+u_right,2,1)*Fe(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
-                     wee(ij+u_right,3,1)*Fe(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
-                     wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
-                     wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
-                     wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
-                     wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
-                     wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
-                     wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
-                     wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
+                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
+                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
+                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
+                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
+                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
+                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
+                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
+                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
+                     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)
                 
-                uu = wee(ij+u_lup,1,1)*Fe(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
-                     wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
-                     wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
-                     wee(ij+u_lup,4,1)*Fe(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
-                     wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
-                     wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
-                     wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
-                     wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
-                     wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
-                     wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
+                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
+                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
+                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
+                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
+                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
+                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
+                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
+                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
+                     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)
 
                 
-                uu = wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
-                     wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
-                     wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
-                     wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
-                     wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
-                     wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
-                     wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
-                     wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
-                     wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
-                     wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
+                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
+                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
+                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
+                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
+                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
+                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
+                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
+                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
+                     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)
                 
@@ -472,39 +478,39 @@
                 ij=(j-1)*iim+i
 
-                uu = wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+                           &
-                     wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+                           &
-                     wee(ij+u_right,3,1)*Fe(ij+u_left,l)+                          &
-                     wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+                         &
-                     wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+                         & 
-                     wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+                 &
-                     wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+                 &
-                     wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+                 &
-                     wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+                   &
-                     wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l)
+                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
+                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
+                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
+                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
+                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
+                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
+                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
+                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
+                     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)
                 
                 
-                uu = wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+                        &
-                     wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+                       &
-                     wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+                       &
-                     wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+                       &
-                     wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+                         & 
-                     wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+                 &
-                     wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+                   &
-                     wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+                   &
-                     wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+                  &
-                     wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l)
+                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
+                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
+                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
+                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
+                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
+                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
+                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
+                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
+                     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)
 
-                uu = wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+                      &
-                     wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+                      &
-                     wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+                        &
-                     wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+                        &
-                     wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+                       & 
-                     wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+                &
-                     wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+               &
-                     wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+              &
-                     wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+              &
-                     wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l)
+                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
+                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
+                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
+                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
+                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
+                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
+                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
+                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
+                     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)
      
@@ -606,19 +612,19 @@
          ! ww>0 <=> upward transport
 
-        ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
+        ww            =   0.5 * wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
         dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
 
-        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1))
+        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_right,l+1))
         uu  = u(ij+u_right,l+1) - u(ij+u_right,l)  
         du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
         du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1)))
 
-        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_lup,l+1))
+        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_lup,l+1))
         uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l)  
         du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
         du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1)))
 
-        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_ldown,l+1))
+        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_ldown,l+1))
         uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l)  
         du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
@@ -631,14 +637,12 @@
 !!$OMP BARRIER
 !!$OMP MASTER  
-    DEALLOCATE(theta)  ! potential temperature
-    DEALLOCATE(pk)   ! Exner function
+    DEALLOCATE(theta)   ! potential temperature
+    DEALLOCATE(pk)      ! Exner function
     DEALLOCATE(pks)
     DEALLOCATE(alpha)
     DEALLOCATE(beta)
-    DEALLOCATE(phi)    ! geopotential
-    DEALLOCATE(Fe)   ! mass flux   
-    DEALLOCATE(Ftheta) ! theta flux   
-    DEALLOCATE(convm)    ! mass flux convergence
-    DEALLOCATE(w)       ! vertical velocity      
+    DEALLOCATE(phi)     ! geopotential
+    DEALLOCATE(Ftheta)  ! theta flux   
+    DEALLOCATE(convm)   ! mass flux convergence
     DEALLOCATE(berni)   ! bernouilli term
 !!$OMP END MASTER
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 133)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 134)
@@ -46,4 +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
 !  INTEGER :: itaumax
 !  REAL(rstd) ::write_period
@@ -78,4 +79,5 @@
   CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
   CALL allocate_field(f_rhodz,field_t,type_real,llm)
+! Mass fluxes
   CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
   CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
@@ -138,9 +140,8 @@
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
-     hfluxt=f_hfluxt(ind); hfluxt=0. ! clear accumulated fluxes
-     wfluxt=f_hfluxt(ind); wfluxt=0.
      rhodz=f_rhodz(ind); ps=f_ps(ind)
      CALL compute_rhodz(ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
   END DO
+  fluxt_zero=.FALSE.
 
   DO it=0,itaumax
@@ -155,10 +156,7 @@
 
     DO stage=1,nb_stage
-!       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
-!            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
-!            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
        CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
             f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
-            f_dps, f_dtheta_rhodz, f_du)
+            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
        SELECT CASE (scheme)
        CASE(euler)
@@ -185,10 +183,5 @@
     IF(MOD(it+1,itau_adv)==0) THEN
 !       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
-       DO ind=1,ndomain
-          hfluxt=f_hfluxt(ind)
-          wfluxt=f_hfluxt(ind)
-          hfluxt=0.               ! clear accumulated fluxes
-          wfluxt=0.
-       END DO
+       fluxt_zero=.TRUE.
     END IF
 
@@ -209,5 +202,5 @@
           hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
           wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
-          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt)
+          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero)
        END IF
        u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
@@ -242,5 +235,5 @@
             hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
             wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind)
-            CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt)
+            CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero)
          END IF
       END DO
@@ -366,11 +359,18 @@
   END SUBROUTINE compute_rhodz
 
-  SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau)
+  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)
     REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
     REAL(rstd), INTENT(IN) :: tau
-    hfluxt = hfluxt + tau*hflux
-    wfluxt = wfluxt + tau*wflux
+    LOGICAL, INTENT(INOUT) :: fluxt_zero
+    IF(fluxt_zero) THEN
+       fluxt_zero=.FALSE.
+       hfluxt = tau*hflux
+       wfluxt = tau*wflux
+    ELSE
+       hfluxt = hfluxt + tau*hflux
+       wfluxt = wfluxt + tau*wflux
+    END IF
   END SUBROUTINE accumulate_fluxes
   
