Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 132)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 133)
@@ -1,4 +1,10 @@
 MODULE timeloop_gcm_mod
-  
+
+  PRIVATE
+  
+  PUBLIC :: timeloop
+
+  INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
+
 CONTAINS
   
@@ -24,4 +30,5 @@
   TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
   TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
+  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
 
   REAL(rstd),POINTER :: phis(:)
@@ -29,13 +36,14 @@
   REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
   REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
-  REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
+  REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
   REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
   REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
   REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
+  REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
 
 !  REAL(rstd) :: dt, run_length
   INTEGER :: ind
-  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period
-  CHARACTER(len=255) :: scheme
+  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
+  CHARACTER(len=255) :: scheme_name
 !  INTEGER :: itaumax
 !  REAL(rstd) ::write_period
@@ -70,14 +78,21 @@
   CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
   CALL allocate_field(f_rhodz,field_t,type_real,llm)
-
-  scheme='runge_kutta'
-  CALL getin('scheme',scheme)
-
-  SELECT CASE (TRIM(scheme))
+  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
+  CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
+  CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
+
+  scheme_name='runge_kutta'
+  CALL getin('scheme',scheme_name)
+
+  SELECT CASE (TRIM(scheme_name))
   CASE('euler')
+     scheme=euler
      nb_stage=1
   CASE ('runge_kutta')
+     scheme=rk4
      nb_stage=4
   CASE ('leapfrog_matsuno')
+     scheme=mlf
      matsuno_period=5
      CALL getin('matsuno_period',matsuno_period)
@@ -88,5 +103,5 @@
      CALL allocate_field(f_um2,field_u,type_real,llm)
   CASE default
-     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
+     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
           ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
      STOP
@@ -120,5 +135,12 @@
   CALL transfert_request(f_q,req_i1) 
 
-  CALL advect_tracer_rhodz(f_ps, f_rhodz)  ! save rhodz for transport scheme before dynamics update ps
+  DO ind=1,ndomain
+     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
 
   DO it=0,itaumax
@@ -135,16 +157,14 @@
 !       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
 !            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
-!            f_Fe, f_W, f_dps, f_dtheta_rhodz, f_du)
+!            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)
-       SELECT CASE (TRIM(scheme))
-       CASE('euler')
-          CALL  euler_scheme(.TRUE.)
-          
-       CASE ('runge_kutta')
+       SELECT CASE (scheme)
+       CASE(euler)
+          CALL euler_scheme(.TRUE.)
+       CASE (rk4)
           CALL rk_scheme(stage)
-          
-       CASE ('leapfrog_matsuno')
+       CASE (mlf)
           CALL  leapfrog_matsuno_scheme(stage)
           
@@ -155,11 +175,6 @@
           !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
           !        CALL adam_bashforth_scheme
-       CASE default
-          !        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
-          !               ' > options are <euler>, <runge_kutta>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>'
-          PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             &
-               ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
+       CASE DEFAULT
           STOP
-          
        END SELECT
     END DO
@@ -168,5 +183,14 @@
     CALL euler_scheme(.FALSE.)
 
-!    CALL advect_tracer(f_Fe,f_W,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
+    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
+    END IF
+
 !    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
     ENDDO
@@ -183,4 +207,7 @@
           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)
        END IF
        u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
@@ -189,5 +216,5 @@
        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
     ENDDO
-      
+
     END SUBROUTINE Euler_scheme
 
@@ -212,6 +239,10 @@
          u(:,:)=um1(:,:)+tau*du(:,:)
          theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
+         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)
+         END IF
       END DO
-
     END SUBROUTINE RK_scheme
 
@@ -318,4 +349,29 @@
 
   END SUBROUTINE timeloop    
+
+  SUBROUTINE compute_rhodz(ps, rhodz)
+    USE icosa
+    USE disvert_mod
+    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
+    REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm)
+    INTEGER :: l,i,j,ij
+    DO l = 1, llm
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
+          ENDDO
+       ENDDO
+    ENDDO
+  END SUBROUTINE compute_rhodz
+
+  SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau)
+    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
+  END SUBROUTINE accumulate_fluxes
   
 END MODULE timeloop_gcm_mod
