Index: codes/icosagcm/trunk/src/diagnostics/check_conserve.f90
===================================================================
--- codes/icosagcm/trunk/src/diagnostics/check_conserve.f90	(revision 581)
+++ codes/icosagcm/trunk/src/diagnostics/check_conserve.f90	(revision 599)
@@ -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/trunk/src/diagnostics/diagflux.F90
===================================================================
--- codes/icosagcm/trunk/src/diagnostics/diagflux.F90	(revision 599)
+++ codes/icosagcm/trunk/src/diagnostics/diagflux.F90	(revision 599)
@@ -0,0 +1,160 @@
+MODULE diagflux_mod
+  USE icosa
+  USE omp_para
+  IMPLICIT NONE
+  SAVE
+  PRIVATE
+
+  TYPE(t_field), POINTER, PUBLIC :: &
+       f_masst(:), f_qmasst(:), & ! time-averaged mass, tracer mass, 
+       f_massfluxt(:), f_qfluxt(:), & ! time-integrated mass flux and tracer flux
+       f_qfluxt_lon(:), f_qfluxt_lat(:), & ! scalar flux reconstructed at cell centers
+       f_epot(:), f_ekin(:), f_enthalpy(:), & ! time-averaged potential E, kinetic E and enthalpy
+       f_epotfluxt(:), f_ekinfluxt(:), f_enthalpyfluxt(:) ! time averaged 'fluxes' of epot, ekin and enthalpy
+  LOGICAL :: diagflux_on
+  !$OMP THREADPRIVATE(diagflux_on)
+
+  PUBLIC :: diagflux_on, init_diagflux, zero_qfluxt, qflux_centered_lonlat, flux_centered_lonlat, diagflux_energy
+
+CONTAINS
+
+  SUBROUTINE init_diagflux
+    USE getin_mod
+    INTEGER :: ll
+    diagflux_on = .FALSE.
+    CALL getin("diagflux", diagflux_on)
+    ll = MERGE(llm,1,diagflux_on)
+    CALL allocate_field(f_masst,         field_t,type_real,ll,       name="masst")
+    CALL allocate_field(f_epot,          field_t,type_real,ll,       name="epot")
+    CALL allocate_field(f_ekin,          field_t,type_real,ll,       name="ekin")
+    CALL allocate_field(f_enthalpy,      field_t,type_real,ll,       name="enthalpy")
+    CALL allocate_field(f_qmasst,        field_t,type_real,ll,nqtot, name="qmasst")
+    CALL allocate_field(f_massfluxt,     field_u,type_real,ll,       name="massfluxt")
+    CALL allocate_field(f_epotfluxt,     field_u,type_real,ll,       name="epotfluxt")
+    CALL allocate_field(f_ekinfluxt,     field_u,type_real,ll,       name="ekinfluxt")
+    CALL allocate_field(f_enthalpyfluxt, field_u,type_real,ll,       name="enthalpyfluxt")
+    CALL allocate_field(f_qfluxt,        field_u,type_real,ll,nqtot, name="qfluxt")
+    CALL allocate_field(f_qfluxt_lon,    field_t,type_real,ll,nqtot, name="qfluxt_lon")
+    CALL allocate_field(f_qfluxt_lat,    field_t,type_real,ll,nqtot, name="qfluxt_lat")
+    IF(diagflux_on) CALL zero_qfluxt
+  END SUBROUTINE init_diagflux
+
+#define ZERO2(field) buf2=field(ind) ; buf2(:,ll_begin:ll_end)=0.
+#define ZERO3(field) buf3=field(ind) ; buf3(:,ll_begin:ll_end,:)=0.
+
+  SUBROUTINE zero_qfluxt
+    INTEGER :: ind
+    REAL(rstd), POINTER :: buf2(:,:),buf3(:,:,:)
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       ZERO2(f_masst)
+       ZERO2(f_epot)
+       ZERO2(f_ekin)
+       ZERO2(f_enthalpy)
+       ZERO3(f_qmasst)
+       ZERO2(f_massfluxt)
+       ZERO2(f_epotfluxt)
+       ZERO2(f_ekinfluxt)
+       ZERO2(f_enthalpyfluxt)
+       ZERO3(f_qfluxt)
+    END DO
+  END SUBROUTINE zero_qfluxt
+
+!------------------------------------ Reconstruct fluxes at cell centers ---------------------------------------
+
+  SUBROUTINE qflux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
+    REAL(rstd), INTENT(IN) :: scale
+    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
+    REAL(rstd), POINTER :: flux(:,:,:), flux_lon(:,:,:), flux_lat(:,:,:)
+    INTEGER :: ind, itrac
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       flux=f_flux(ind)
+       flux_lon=f_flux_lon(ind)
+       flux_lat=f_flux_lat(ind)
+       DO itrac=1,nqtot
+          CALL compute_flux_centered_lonlat(scale, flux(:,:,itrac), flux_lon(:,:,itrac), flux_lat(:,:,itrac))
+       END DO
+    END DO
+  END SUBROUTINE qflux_centered_lonlat
+  
+  SUBROUTINE flux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
+    REAL(rstd), INTENT(IN) :: scale
+    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
+    REAL(rstd), POINTER :: flux(:,:), flux_lon(:,:), flux_lat(:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       flux=f_flux(ind)
+       flux_lon=f_flux_lon(ind)
+       flux_lat=f_flux_lat(ind)
+       CALL compute_flux_centered_lonlat(scale, flux, flux_lon, flux_lat)
+    END DO
+  END SUBROUTINE flux_centered_lonlat
+  
+  SUBROUTINE compute_flux_centered_lonlat(scale, flux, flux_lon, flux_lat)
+    USE wind_mod
+    REAL(rstd), INTENT(IN) :: scale
+    REAL(rstd), INTENT(IN) :: flux(3*iim*jjm,llm)
+    REAL(rstd), INTENT(OUT) :: flux_lon(iim*jjm,llm), flux_lat(iim*jjm,llm)
+    REAL(rstd) :: flux_3d(iim*jjm,llm,3)
+    CALL compute_flux_centered(scale, flux, flux_3d)
+    CALL compute_wind_centered_lonlat_compound(flux_3d, flux_lon, flux_lat)
+  END SUBROUTINE compute_flux_centered_lonlat
+
+!------------------------------------ Compute energy fluxes ---------------------------------------
+
+  SUBROUTINE diagflux_energy(frac, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt)
+    REAL(rstd), INTENT(IN) :: frac
+    TYPE(t_field),POINTER :: f_phis(:),f_rhodz(:),f_theta_rhodz(:),f_u(:), f_geopot(:), f_theta(:), f_hfluxt(:)
+    REAL(rstd), POINTER :: phis(:), rhodz(:,:), theta_rhodz(:,:,:), u(:,:), &
+         geopot(:,:), pk(:,:,:), hfluxt(:,:), &
+         epot(:,:), ekin(:,:), enthalpy(:,:), &
+         epotflux(:,:), ekinflux(:,:), enthalpyflux(:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       hfluxt = f_hfluxt(ind)
+       phis = f_phis(ind)
+       rhodz = f_rhodz(ind)
+       theta_rhodz = f_theta_rhodz(ind)
+       u = f_u(ind)
+       geopot = f_geopot(ind)
+       pk = f_theta(ind) ! buffer
+       epot = f_epot(ind)
+       ekin = f_ekin(ind)
+       enthalpy = f_enthalpy(ind)
+       epotflux = f_epotfluxt(ind)
+       ekinflux = f_ekinfluxt(ind)
+       enthalpyflux = f_enthalpyfluxt(ind)
+       CALL compute_diagflux_energy(frac,hfluxt, phis,rhodz,theta_rhodz,u, geopot,pk, epot,ekin,enthalpy, epotflux, ekinflux, enthalpyflux)
+    END DO
+  END SUBROUTINE diagflux_energy
+
+  SUBROUTINE compute_diagflux_energy(frac, massflux, phis,rhodz,theta_rhodz,u, geopot,pk, epot,ekin,enthalpy, epot_flux, ekin_flux, enthalpy_flux)
+    USE disvert_mod, ONLY : ptop
+    REAL(rstd), INTENT(IN) :: frac
+    REAL(rstd), INTENT(IN) :: massflux(3*iim*jjm,llm), u(3*iim*jjm,llm),&
+                              phis(iim*jjm), rhodz(iim*jjm,llm), theta_rhodz(iim*jjm,llm,nqtot)
+    REAL(rstd), INTENT(INOUT) :: geopot(iim*jjm,llm+1), pk(iim*jjm,llm) ! pk = buffer
+    REAL(rstd), INTENT(INOUT), DIMENSION(iim*jjm, llm)   ::  epot, ekin, enthalpy
+    REAL(rstd), INTENT(INOUT), DIMENSION(3*iim*jjm, llm) ::  epot_flux, ekin_flux, enthalpy_flux    
+    REAL(rstd) :: energy, p_ik, theta_ik, temp_ik, gv, Rd 
+    INTEGER :: ij, l, ij_omp_begin_ext, ij_omp_end_ext
+    Rd = kappa*cpp
+    ! even if loops are of the _ext variant, we still need halo exchanges before reconstructing fluxes at cell centers
+    ! => loop over interior region
+    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin_ext,ij_omp_end_ext)
+    ij_omp_begin_ext = ij_omp_begin_ext+ij_begin-1
+    ij_omp_end_ext = ij_omp_end_ext+ij_begin-1
+#include "../kernels/energy_fluxes.k90"
+  END SUBROUTINE compute_diagflux_energy
+
+END MODULE diagflux_mod
Index: codes/icosagcm/trunk/src/diagnostics/observable.f90
===================================================================
--- codes/icosagcm/trunk/src/diagnostics/observable.f90	(revision 581)
+++ codes/icosagcm/trunk/src/diagnostics/observable.f90	(revision 599)
@@ -6,6 +6,5 @@
   TYPE(t_field),POINTER, SAVE :: f_buf_i(:), &
        f_buf_Fel(:), f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH
-       f_buf_ulon(:), f_buf_ulat(:), &
-       f_buf_u3d(:) ! unused, remove ?
+       f_buf_ulon(:), f_buf_ulat(:)
   TYPE(t_field),POINTER, SAVE :: f_buf1_i(:), f_buf2_i(:)
   TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
@@ -24,5 +23,4 @@
     CALL allocate_field(f_buf2_i,   field_t,type_real,llm,name="buffer2_i")
     CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
-    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
     CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon")
     CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
@@ -49,4 +47,5 @@
     USE theta2theta_rhodz_mod
     USE omega_mod
+    USE diagflux_mod
     LOGICAL, INTENT(IN) :: init
     INTEGER :: l
@@ -160,4 +159,36 @@
     END IF
 
+    IF(.NOT. init) THEN
+       IF(diagflux_on) THEN
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_buf_ulon, f_buf_ulat)
+          CALL output_field("mass_t", f_masst)
+          CALL output_field("massflux_lon",f_buf_ulon)
+          CALL output_field("massflux_lat",f_buf_ulat)
+
+          CALL transfert_request(f_epotfluxt,req_e1_vect) 
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_epotfluxt, f_buf_ulon, f_buf_ulat)
+          CALL output_field("epot_t", f_epot)
+          CALL output_field("epotflux_lon",f_buf_ulon)
+          CALL output_field("epotflux_lat",f_buf_ulat)
+
+          CALL transfert_request(f_ekinfluxt,req_e1_vect) 
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_ekinfluxt, f_buf_ulon, f_buf_ulat)
+          CALL output_field("ekin_t", f_ekin)
+          CALL output_field("ekinflux_lon",f_buf_ulon)
+          CALL output_field("ekinflux_lat",f_buf_ulat)
+
+          CALL transfert_request(f_enthalpyfluxt,req_e1_vect) 
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_enthalpyfluxt, f_buf_ulon, f_buf_ulat)
+          CALL output_field("enthalpy_t", f_enthalpy)
+          CALL output_field("enthalpyflux_lon",f_buf_ulon)
+          CALL output_field("enthalpyflux_lat",f_buf_ulat)
+
+          CALL qflux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat)
+          CALL output_field("qmass_t", f_qmasst)
+          CALL output_field("qflux_lon",f_qfluxt_lon)
+          CALL output_field("qflux_lat",f_qfluxt_lat)
+          CALL zero_qfluxt  ! restart accumulating fluxes
+       END IF
+    END IF
   END SUBROUTINE write_output_fields_basic
   
Index: codes/icosagcm/trunk/src/diagnostics/wind.F90
===================================================================
--- codes/icosagcm/trunk/src/diagnostics/wind.F90	(revision 599)
+++ codes/icosagcm/trunk/src/diagnostics/wind.F90	(revision 599)
@@ -0,0 +1,355 @@
+MODULE wind_mod
+  USE omp_para
+  USE icosa
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: compute_wind_centered, compute_flux_centered, &
+       compute_wind_centered_lonlat_compound, compute_wind2d_perp_from_lonlat_compound, &
+       compute_wind_centered_from_lonlat_compound, compute_wind_perp_from_lonlat_compound, &
+       un2ulonlat, ulonlat2un
+       
+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    
+    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       ulon=f_ulon(ind)
+       ulat=f_ulat(ind)
+       CALL compute_un2ulonlat(u,ulon, ulat)
+    END DO
+
+  END SUBROUTINE un2ulonlat
+
+  SUBROUTINE ulonlat2un(f_ulon, f_ulat,f_u)
+    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
+    
+    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       ulon=f_ulon(ind)
+       ulat=f_ulat(ind)
+       CALL compute_ulonlat2un(ulon, ulat,u)
+    END DO
+
+  END SUBROUTINE ulonlat2un
+ 
+  SUBROUTINE compute_wind_centered(ue,ucenter)
+  REAL(rstd) :: ue(3*iim*jjm,llm)
+  REAL(rstd) :: ucenter(iim*jjm,llm,3)
+  INTEGER :: ij,l
+  REAL(rstd), PARAMETER :: scale=1.
+  REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz
+#include "../kernels/wind_centered.k90"
+ END SUBROUTINE compute_wind_centered
+ 
+  SUBROUTINE compute_flux_centered(scale,ue,ucenter)
+  REAL(rstd), INTENT(IN) :: scale
+  REAL(rstd) :: ue(3*iim*jjm,llm)
+  REAL(rstd) :: ucenter(iim*jjm,llm,3)
+  INTEGER :: ij,l
+  REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz
+#include "../kernels/flux_centered.k90"
+  END SUBROUTINE compute_flux_centered
+ 
+  
+ SUBROUTINE compute_wind_on_edge(ue,uedge)
+  REAL(rstd) :: ue(3*iim*jjm,llm)
+  REAL(rstd) :: uedge(3*iim*jjm,llm,3)
+
+  REAL(rstd) :: ut(3*iim*jjm,llm)
+  INTEGER :: i,j,ij,l     
+    
+    CALL compute_tangential_compound(ue,ut)
+  
+    DO l=ll_begin,ll_end
+      DO j=jj_begin,jj_end
+        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)
+        ENDDO
+      ENDDO
+    ENDDO
+ 
+ END SUBROUTINE compute_wind_on_edge
+ 
+ 
+ 
+ SUBROUTINE compute_tangential_compound(ue,ut)
+  REAL(rstd) :: ue(3*iim*jjm,llm)
+  REAL(rstd) :: ut(3*iim*jjm,llm)
+  INTEGER :: i,j,l,ij
+    
+  DO l=ll_begin,ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+   
+        ut(ij+u_right,l) = 1/de(ij+u_right) *                                            &  
+                         ( wee(ij+u_right,1,1)*ue(ij+u_rup,l)+                           &
+                           wee(ij+u_right,2,1)*ue(ij+u_lup,l)+                           &
+                           wee(ij+u_right,3,1)*ue(ij+u_left,l)+                          &
+                           wee(ij+u_right,4,1)*ue(ij+u_ldown,l)+                         &
+                           wee(ij+u_right,5,1)*ue(ij+u_rdown,l)+                         & 
+                           wee(ij+u_right,1,2)*ue(ij+t_right+u_ldown,l)+                 &
+                           wee(ij+u_right,2,2)*ue(ij+t_right+u_rdown,l)+                 &
+                           wee(ij+u_right,3,2)*ue(ij+t_right+u_right,l)+                 &
+                           wee(ij+u_right,4,2)*ue(ij+t_right+u_rup,l)+                   &
+                           wee(ij+u_right,5,2)*ue(ij+t_right+u_lup,l) )   
+     
+        ut(ij+u_lup,l) =  1/de(ij+u_lup) *                                           &  
+                         ( wee(ij+u_lup,1,1)*ue(ij+u_left,l)+                        &
+                           wee(ij+u_lup,2,1)*ue(ij+u_ldown,l)+                       &
+                           wee(ij+u_lup,3,1)*ue(ij+u_rdown,l)+                       &
+                           wee(ij+u_lup,4,1)*ue(ij+u_right,l)+                       &
+                           wee(ij+u_lup,5,1)*ue(ij+u_rup,l)+                         & 
+                           wee(ij+u_lup,1,2)*ue(ij+t_lup+u_right,l)+                 &
+                           wee(ij+u_lup,2,2)*ue(ij+t_lup+u_rup,l)+                   &
+                           wee(ij+u_lup,3,2)*ue(ij+t_lup+u_lup,l)+                   &
+                           wee(ij+u_lup,4,2)*ue(ij+t_lup+u_left,l)+                  &
+                           wee(ij+u_lup,5,2)*ue(ij+t_lup+u_ldown,l) )
+
+    
+        ut(ij+u_ldown,l) = 1/de(ij+u_ldown) *   &  
+                         ( wee(ij+u_ldown,1,1)*ue(ij+u_rdown,l)+                      &
+                           wee(ij+u_ldown,2,1)*ue(ij+u_right,l)+                      &
+                           wee(ij+u_ldown,3,1)*ue(ij+u_rup,l)+                        &
+                           wee(ij+u_ldown,4,1)*ue(ij+u_lup,l)+                        &
+                           wee(ij+u_ldown,5,1)*ue(ij+u_left,l)+                       & 
+                           wee(ij+u_ldown,1,2)*ue(ij+t_ldown+u_lup,l)+                &
+                           wee(ij+u_ldown,2,2)*ue(ij+t_ldown+u_left,l)+               &
+                           wee(ij+u_ldown,3,2)*ue(ij+t_ldown+u_ldown,l)+              &
+                           wee(ij+u_ldown,4,2)*ue(ij+t_ldown+u_rdown,l)+              &
+                           wee(ij+u_ldown,5,2)*ue(ij+t_ldown+u_right,l) )  
+        
+        ENDDO
+      ENDDO
+    ENDDO
+                        
+ END SUBROUTINE compute_tangential_compound
+ 
+! 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
+      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,:)
+        ENDDO
+      ENDDO
+    ENDDO
+ 
+  END SUBROUTINE compute_wind_from_lonlat_compound
+ 
+  SUBROUTINE compute_wind_centered_from_lonlat_compound(ulon, ulat, u)
+  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
+      DO j=jj_begin-1,jj_end+1
+        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 
+  END SUBROUTINE compute_wind_centered_from_lonlat_compound
+ 
+  SUBROUTINE compute_wind2D_from_lonlat_compound(ulon, ulat, u)
+  REAL(rstd) :: u(3*iim*jjm,3)
+  REAL(rstd) :: ulon(3*iim*jjm)
+  REAL(rstd) :: ulat(3*iim*jjm)
+  
+  INTEGER :: i,j,ij
+  
+  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,:)=ulon(ij+u_right)*elon_e(ij+u_right,:)+ ulat(ij+u_right)*elat_e(ij+u_right,:)
+        u(ij+u_lup,:)=ulon(ij+u_lup)*elon_e(ij+u_lup,:) + ulat(ij+u_lup)*elat_e(ij+u_lup,:)
+        u(ij+u_ldown,:)=ulon(ij+u_ldown)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown)*elat_e(ij+u_ldown,:)
+     ENDDO
+  ENDDO
+    
+  END SUBROUTINE compute_wind2D_from_lonlat_compound
+ 
+  SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up)
+  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,llm,3)
+
+  INTEGER :: i,j,ij,l     
+  
+   CALL compute_wind_from_lonlat_compound(ulon, ulat, u)
+
+    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
+          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
+    ENDDO
+ 
+  END SUBROUTINE compute_wind_perp_from_lonlat_compound
+   
+  SUBROUTINE compute_wind2D_perp_from_lonlat_compound(ulon, ulat, up)
+  REAL(rstd) :: up(3*iim*jjm)
+  REAL(rstd) :: ulon(3*iim*jjm)
+  REAL(rstd) :: ulat(3*iim*jjm)
+  REAL(rstd) :: u(3*iim*jjm,3)
+
+  INTEGER :: i,j,ij 
+  
+  CALL compute_wind2D_from_lonlat_compound(ulon, ulat, u)
+  DO j=jj_begin-1,jj_end+1
+     DO i=ii_begin-1,ii_end+1
+        ij=(j-1)*iim+i
+        up(ij+u_right)=sum(u(ij+u_right,:)*ep_e(ij+u_right,:))
+        up(ij+u_lup)=sum(u(ij+u_lup,:)*ep_e(ij+u_lup,:))
+        up(ij+u_ldown)=sum(u(ij+u_ldown,:)*ep_e(ij+u_ldown,:))
+     ENDDO
+  ENDDO
+   
+  END SUBROUTINE compute_wind2D_perp_from_lonlat_compound
+   
+ SUBROUTINE compute_wind_centered_lonlat_compound(uc, ulon, ulat)
+  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,:)) 
+        ENDDO
+      ENDDO
+    ENDDO
+ 
+ END SUBROUTINE compute_wind_centered_lonlat_compound
+
+ SUBROUTINE compute_wind_centered_from_wind_lonlat_centered(ulon, ulat,uc)
+  REAL(rstd) :: ulon(iim*jjm,llm)
+  REAL(rstd) :: ulat(iim*jjm,llm)
+  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
+    ENDDO
+ 
+ END SUBROUTINE compute_wind_centered_from_wind_lonlat_centered
+
+ SUBROUTINE compute_wind_perp_from_wind_centered(uc,un)
+
+  IMPLICIT NONE
+  REAL(rstd),INTENT(IN)   :: uc(iim*jjm,llm,3)
+  REAL(rstd),INTENT(OUT)  :: un(3*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
+          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
+    ENDDO
+ 
+ END SUBROUTINE compute_wind_perp_from_wind_centered
+
+
+ SUBROUTINE compute_un2ulonlat(un, ulon, ulat)
+  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,llm,3)
+    
+  CALL compute_wind_centered(un,uc)  
+  CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
+ 
+ END SUBROUTINE compute_un2ulonlat
+
+ SUBROUTINE compute_ulonlat2un(ulon, ulat,un)
+  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,llm,3)
+    
+    CALL compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc)  
+    CALL compute_wind_perp_from_wind_centered(uc, un)
+ 
+ END SUBROUTINE compute_ulonlat2un
+
+
+END MODULE wind_mod
Index: codes/icosagcm/trunk/src/diagnostics/wind.f90
===================================================================
--- codes/icosagcm/trunk/src/diagnostics/wind.f90	(revision 581)
+++ 	(revision )
@@ -1,408 +1,0 @@
-MODULE wind_mod
-
-CONTAINS
-
-  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat)
-  USE icosa
-  IMPLICIT NONE
-    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
-    
-    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
-    INTEGER :: ind
-
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       u=f_u(ind)
-       ulon=f_ulon(ind)
-       ulat=f_ulat(ind)
-       CALL compute_un2ulonlat(u,ulon, ulat)
-    END DO
-
-  END SUBROUTINE un2ulonlat
-
-  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
-    
-    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
-    INTEGER :: ind
-
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       u=f_u(ind)
-       ulon=f_ulon(ind)
-       ulat=f_ulat(ind)
-       CALL compute_ulonlat2un(ulon, ulat,u)
-    END DO
-
-  END SUBROUTINE ulonlat2un
- 
-  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)
-  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
-          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,:))          &
-                         + ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:))          &
-                         + ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:))    &
-                         + ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))&
-                         + ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:)))
-        ENDDO
-      ENDDO
-    ENDDO
- 
- END SUBROUTINE compute_wind_centered
- 
-  
- 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) :: ut(3*iim*jjm,llm)
-  INTEGER :: i,j,ij,l     
-    
-    CALL compute_tangential_compound(ue,ut)
-  
-    DO l=ll_begin,ll_end
-      DO j=jj_begin,jj_end
-        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)
-        ENDDO
-      ENDDO
-    ENDDO
- 
- END SUBROUTINE compute_wind_on_edge
- 
- 
- 
- 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)
-  INTEGER :: i,j,l,ij
-    
-  DO l=ll_begin,ll_end
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-   
-        ut(ij+u_right,l) = 1/de(ij+u_right) *                                            &  
-                         ( wee(ij+u_right,1,1)*ue(ij+u_rup,l)+                           &
-                           wee(ij+u_right,2,1)*ue(ij+u_lup,l)+                           &
-                           wee(ij+u_right,3,1)*ue(ij+u_left,l)+                          &
-                           wee(ij+u_right,4,1)*ue(ij+u_ldown,l)+                         &
-                           wee(ij+u_right,5,1)*ue(ij+u_rdown,l)+                         & 
-                           wee(ij+u_right,1,2)*ue(ij+t_right+u_ldown,l)+                 &
-                           wee(ij+u_right,2,2)*ue(ij+t_right+u_rdown,l)+                 &
-                           wee(ij+u_right,3,2)*ue(ij+t_right+u_right,l)+                 &
-                           wee(ij+u_right,4,2)*ue(ij+t_right+u_rup,l)+                   &
-                           wee(ij+u_right,5,2)*ue(ij+t_right+u_lup,l) )   
-     
-        ut(ij+u_lup,l) =  1/de(ij+u_lup) *                                           &  
-                         ( wee(ij+u_lup,1,1)*ue(ij+u_left,l)+                        &
-                           wee(ij+u_lup,2,1)*ue(ij+u_ldown,l)+                       &
-                           wee(ij+u_lup,3,1)*ue(ij+u_rdown,l)+                       &
-                           wee(ij+u_lup,4,1)*ue(ij+u_right,l)+                       &
-                           wee(ij+u_lup,5,1)*ue(ij+u_rup,l)+                         & 
-                           wee(ij+u_lup,1,2)*ue(ij+t_lup+u_right,l)+                 &
-                           wee(ij+u_lup,2,2)*ue(ij+t_lup+u_rup,l)+                   &
-                           wee(ij+u_lup,3,2)*ue(ij+t_lup+u_lup,l)+                   &
-                           wee(ij+u_lup,4,2)*ue(ij+t_lup+u_left,l)+                  &
-                           wee(ij+u_lup,5,2)*ue(ij+t_lup+u_ldown,l) )
-
-    
-        ut(ij+u_ldown,l) = 1/de(ij+u_ldown) *   &  
-                         ( wee(ij+u_ldown,1,1)*ue(ij+u_rdown,l)+                      &
-                           wee(ij+u_ldown,2,1)*ue(ij+u_right,l)+                      &
-                           wee(ij+u_ldown,3,1)*ue(ij+u_rup,l)+                        &
-                           wee(ij+u_ldown,4,1)*ue(ij+u_lup,l)+                        &
-                           wee(ij+u_ldown,5,1)*ue(ij+u_left,l)+                       & 
-                           wee(ij+u_ldown,1,2)*ue(ij+t_ldown+u_lup,l)+                &
-                           wee(ij+u_ldown,2,2)*ue(ij+t_ldown+u_left,l)+               &
-                           wee(ij+u_ldown,3,2)*ue(ij+t_ldown+u_ldown,l)+              &
-                           wee(ij+u_ldown,4,2)*ue(ij+t_ldown+u_rdown,l)+              &
-                           wee(ij+u_ldown,5,2)*ue(ij+t_ldown+u_right,l) )  
-        
-        ENDDO
-      ENDDO
-    ENDDO
-                        
- 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     
-    
-  
-    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)
-  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,:)
-        ENDDO
-      ENDDO
-    ENDDO
- 
-  END SUBROUTINE compute_wind_from_lonlat_compound
- 
-  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) :: 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-1,jj_end+1
-        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
- 
-  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)
-  REAL(rstd) :: ulat(3*iim*jjm)
-  
-  INTEGER :: i,j,ij
-  
-  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,:)=ulon(ij+u_right)*elon_e(ij+u_right,:)+ ulat(ij+u_right)*elat_e(ij+u_right,:)
-        u(ij+u_lup,:)=ulon(ij+u_lup)*elon_e(ij+u_lup,:) + ulat(ij+u_lup)*elat_e(ij+u_lup,:)
-        u(ij+u_ldown,:)=ulon(ij+u_ldown)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown)*elat_e(ij+u_ldown,:)
-     ENDDO
-  ENDDO
-    
-  END SUBROUTINE compute_wind2D_from_lonlat_compound
- 
-  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)
-
-  INTEGER :: i,j,ij,l     
-  
-   CALL compute_wind_from_lonlat_compound(ulon, ulat, u)
-
-    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
-          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
-    ENDDO
- 
-  END SUBROUTINE compute_wind_perp_from_lonlat_compound
-   
-  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)
-  REAL(rstd) :: ulat(3*iim*jjm)
-  REAL(rstd) :: u(3*iim*jjm,3)
-
-  INTEGER :: i,j,ij 
-  
-  CALL compute_wind2D_from_lonlat_compound(ulon, ulat, u)
-  DO j=jj_begin-1,jj_end+1
-     DO i=ii_begin-1,ii_end+1
-        ij=(j-1)*iim+i
-        up(ij+u_right)=sum(u(ij+u_right,:)*ep_e(ij+u_right,:))
-        up(ij+u_lup)=sum(u(ij+u_lup,:)*ep_e(ij+u_lup,:))
-        up(ij+u_ldown)=sum(u(ij+u_ldown,:)*ep_e(ij+u_ldown,:))
-     ENDDO
-  ENDDO
-   
-  END SUBROUTINE compute_wind2D_perp_from_lonlat_compound
-   
- 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) :: 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,:)) 
-        ENDDO
-      ENDDO
-    ENDDO
- 
- END SUBROUTINE compute_wind_centered_lonlat_compound
-
- 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,:)
-        ENDDO
-      ENDDO
-    ENDDO
- 
- 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(OUT)  :: un(3*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
-          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
-    ENDDO
- 
- END SUBROUTINE compute_wind_perp_from_wind_centered
-
-
- 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)
-    
-  CALL compute_wind_centered(un,uc)  
-  CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
- 
- END SUBROUTINE compute_un2ulonlat
-
- 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)
-    
-    CALL compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc)  
-    CALL compute_wind_perp_from_wind_centered(uc, un)
- 
- END SUBROUTINE compute_ulonlat2un
-
-
-END MODULE wind_mod
Index: codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90	(revision 581)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90	(revision 599)
@@ -5,4 +5,6 @@
   USE caldyn_kernels_base_mod
   USE caldyn_kernels_mod
+  USE omp_para
+  USE mpipara
   IMPLICIT NONE
   PRIVATE
@@ -13,9 +15,4 @@
   
   SUBROUTINE init_caldyn
-    USE icosa
-    USE observable_mod
-    USE mpipara
-    USE omp_para
-    IMPLICIT NONE
     CHARACTER(len=255) :: def
     INTEGER            :: ind
@@ -124,7 +121,4 @@
 
   SUBROUTINE allocate_caldyn
-  USE icosa
-  IMPLICIT NONE
-
     CALL allocate_field(f_out_u,field_u,type_real,llm) 
     CALL allocate_field(f_qu,field_u,type_real,llm) 
@@ -142,7 +136,4 @@
 
   SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux)
-    USE icosa
-    USE mpipara
-    USE omp_para
     TYPE(t_field),POINTER :: f_phis(:)
     TYPE(t_field),POINTER :: f_geopot(:)
@@ -186,17 +177,7 @@
   SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
        f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
-    USE icosa
     USE observable_mod
     USE disvert_mod, ONLY : caldyn_eta, eta_mass
-    USE vorticity_mod
-    USE kinetic_mod
-    USE theta2theta_rhodz_mod
-    USE wind_mod
-    USE mpipara
     USE trace
-    USE omp_para
-    USE output_field_mod
-    USE checksum_mod
-    IMPLICIT NONE
     LOGICAL,INTENT(IN)    :: write_out
     TYPE(t_field),POINTER :: f_phis(:)
@@ -375,7 +356,4 @@
 
   SUBROUTINE check_mass_conservation(f_ps,f_dps)
-  USE icosa
-  USE mpipara
-  IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_dps(:)
Index: codes/icosagcm/trunk/src/icosa_init.f90
===================================================================
--- codes/icosagcm/trunk/src/icosa_init.f90	(revision 581)
+++ codes/icosagcm/trunk/src/icosa_init.f90	(revision 599)
@@ -1,5 +1,3 @@
 MODULE icosa_init_mod
-
-
 
 CONTAINS
@@ -22,4 +20,5 @@
   USE restart_mod
   USE etat0_mod
+  USE diagflux_mod
   IMPLICIT NONE
   
@@ -55,5 +54,6 @@
     CALL init_timeloop
     CALL init_physics
-   
+
+    CALL init_diagflux
     CALL timeloop
     CALL switch_omp_no_distrib_level
Index: codes/icosagcm/trunk/src/kernels/advect_horiz.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/advect_horiz.k90	(revision 599)
+++ codes/icosagcm/trunk/src/kernels/advect_horiz.k90	(revision 599)
@@ -0,0 +1,74 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- advect_horiz ----------------------------------
+   ! evaluate tracer field at point cc using piecewise linear reconstruction
+   ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon
+   ! sign*hfluxt>0 iff outgoing
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         IF(ne_right*hfluxt(ij+u_right,l)>0.) THEN
+            qe = qi(ij,l)
+            qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
+            qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
+            qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
+         ELSE
+            qe = qi(ij+t_right,l)
+            qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij+t_right,1))*gradq3d(ij+t_right,l,1)
+            qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij+t_right,2))*gradq3d(ij+t_right,l,2)
+            qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij+t_right,3))*gradq3d(ij+t_right,l,3)
+         END IF
+         qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
+         IF(diagflux_on) qfluxt(ij+u_right,l) = qfluxt(ij+u_right,l)+qflux(ij+u_right,l)
+         IF(ne_lup*hfluxt(ij+u_lup,l)>0.) THEN
+            qe = qi(ij,l)
+            qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
+            qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
+            qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
+         ELSE
+            qe = qi(ij+t_lup,l)
+            qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij+t_lup,1))*gradq3d(ij+t_lup,l,1)
+            qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij+t_lup,2))*gradq3d(ij+t_lup,l,2)
+            qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij+t_lup,3))*gradq3d(ij+t_lup,l,3)
+         END IF
+         qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe
+         IF(diagflux_on) qfluxt(ij+u_lup,l) = qfluxt(ij+u_lup,l)+qflux(ij+u_lup,l)
+         IF(ne_ldown*hfluxt(ij+u_ldown,l)>0.) THEN
+            qe = qi(ij,l)
+            qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
+            qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
+            qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
+         ELSE
+            qe = qi(ij+t_ldown,l)
+            qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij+t_ldown,1))*gradq3d(ij+t_ldown,l,1)
+            qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij+t_ldown,2))*gradq3d(ij+t_ldown,l,2)
+            qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij+t_ldown,3))*gradq3d(ij+t_ldown,l,3)
+         END IF
+         qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
+         IF(diagflux_on) qfluxt(ij+u_ldown,l) = qfluxt(ij+u_ldown,l)+qflux(ij+u_ldown,l)
+      END DO
+   END DO
+   ! update q and, if update_mass, update mass
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         dmass=0.
+         dq=0.
+         dmass = dmass + ne_rup*hfluxt(ij+u_rup,l)
+         dq = dq + ne_rup*qflux(ij+u_rup,l)
+         dmass = dmass + ne_lup*hfluxt(ij+u_lup,l)
+         dq = dq + ne_lup*qflux(ij+u_lup,l)
+         dmass = dmass + ne_left*hfluxt(ij+u_left,l)
+         dq = dq + ne_left*qflux(ij+u_left,l)
+         dmass = dmass + ne_ldown*hfluxt(ij+u_ldown,l)
+         dq = dq + ne_ldown*qflux(ij+u_ldown,l)
+         dmass = dmass + ne_rdown*hfluxt(ij+u_rdown,l)
+         dq = dq + ne_rdown*qflux(ij+u_rdown,l)
+         dmass = dmass + ne_right*hfluxt(ij+u_right,l)
+         dq = dq + ne_right*qflux(ij+u_right,l)
+         newmass = mass(ij,l) - dmass/Ai(ij)
+         qi(ij,l) = (qi(ij,l)*mass(ij,l)-dq/Ai(ij)) / newmass
+         IF(update_mass) mass(ij,l)=newmass
+      END DO
+   END DO
+   !---------------------------- advect_horiz ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/energy_fluxes.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/energy_fluxes.k90	(revision 599)
+++ codes/icosagcm/trunk/src/kernels/energy_fluxes.k90	(revision 599)
@@ -0,0 +1,100 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- energy_fluxes ----------------------------------
+   ! First diagnose geopotential and temperature, column-wise
+   !$OMP BARRIER
+   DO ij=ij_omp_begin_ext,ij_omp_end_ext
+      pk(ij,llm) = ptop + .5*g*rhodz(ij,llm)
+   END DO
+   DO l = llm-1,1,-1
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l)+rhodz(ij,l+1) )
+      END DO
+   END DO
+   SELECT CASE(caldyn_thermo)
+   CASE(thermo_theta)
+      DO l = 1,llm
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            p_ik = pk(ij,l)
+            theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
+            temp_ik = theta_ik*(p_ik/preff)**kappa
+            gv = (g*Rd)*temp_ik/p_ik
+            pk(ij,l) = temp_ik
+            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
+         END DO
+      END DO
+   CASE(thermo_entropy)
+      DO l = 1,llm
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            p_ik = pk(ij,l)
+            theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
+            temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
+            pk(ij,l) = temp_ik
+            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
+         END DO
+      END DO
+   END SELECT
+   !$OMP BARRIER
+   ! Now accumulate energies and energy fluxes
+   ! enthalpy
+   ! NB : at this point temperature is stored in array pk
+   ! pk then serves as buffer to store energy
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         energy = cpp*pk(ij,l)
+         enthalpy(ij,l) = enthalpy(ij,l) + frac*rhodz(ij,l)*energy
+         pk(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         enthalpy_flux(ij+u_right,l) = enthalpy_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
+         enthalpy_flux(ij+u_lup,l) = enthalpy_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
+         enthalpy_flux(ij+u_ldown,l) = enthalpy_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
+      END DO
+   END DO
+   ! potential energy
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         energy = .5*(geopot(ij,l+1)+geopot(ij,l))
+         epot(ij,l) = epot(ij,l) + frac*rhodz(ij,l)*energy
+         pk(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         epot_flux(ij+u_right,l) = epot_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
+         epot_flux(ij+u_lup,l) = epot_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
+         epot_flux(ij+u_ldown,l) = epot_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
+      END DO
+   END DO
+   ! kinetic energy
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         energy=0.d0
+         energy = energy + le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2
+         energy = energy + le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2
+         energy = energy + le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2
+         energy = energy + le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2
+         energy = energy + le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2
+         energy = energy + le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2
+         energy = energy * (.25/Ai(ij))
+         ekin(ij,l) = ekin(ij,l) + frac*rhodz(ij,l)*energy
+         pk(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ekin_flux(ij+u_right,l) = ekin_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
+         ekin_flux(ij+u_lup,l) = ekin_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
+         ekin_flux(ij+u_ldown,l) = ekin_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
+      END DO
+   END DO
+   !---------------------------- energy_fluxes ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/flux_centered.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/flux_centered.k90	(revision 599)
+++ codes/icosagcm/trunk/src/kernels/flux_centered.k90	(revision 599)
@@ -0,0 +1,45 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- flux_centered ----------------------------------
+   ! NB : here the input data is a flux and has already the factor l_e in it
+   ! Input data is rescaled by factor scale
+   ! Perot reconstruction based on Gauss theorem
+   ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         cx=centroid(ij,1)
+         cy=centroid(ij,2)
+         cz=centroid(ij,3)
+         ux=0. ; uy=0. ; uz=0.
+         ue_le = ne_rup*ue(ij+u_rup,l)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_up,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_up,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_up,3))-cz)
+         ue_le = ne_lup*ue(ij+u_lup,l)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_up,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_up,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_up,3))-cz)
+         ue_le = ne_left*ue(ij+u_left,l)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_ldown,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_ldown,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_ldown,3))-cz)
+         ue_le = ne_ldown*ue(ij+u_ldown,l)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_ldown,1)+xyz_v(ij+z_down,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_ldown,2)+xyz_v(ij+z_down,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_ldown,3)+xyz_v(ij+z_down,3))-cz)
+         ue_le = ne_rdown*ue(ij+u_rdown,l)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rdown,1)+xyz_v(ij+z_down,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rdown,2)+xyz_v(ij+z_down,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rdown,3)+xyz_v(ij+z_down,3))-cz)
+         ue_le = ne_right*ue(ij+u_right,l)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_rdown,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_rdown,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_rdown,3))-cz)
+         fac = scale*(1./Ai(ij))
+         ucenter(ij,l,1)=ux*fac
+         ucenter(ij,l,2)=uy*fac
+         ucenter(ij,l,3)=uz*fac
+      END DO
+   END DO
+   !---------------------------- flux_centered ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/gradient.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/gradient.k90	(revision 581)
+++ 	(revision )
@@ -1,12 +1,0 @@
-   !--------------------------------------------------------------------------
-   !---------------------------- gradient ----------------------------------
-   DO l = ll_begin, ll_end
-      !DIR$ SIMD
-      DO ij=ij_begin_ext, ij_end_ext
-         grad(ij+u_right,l) = ne_right*(b(ij+t_right,l)-b(ij,l))
-         grad(ij+u_lup,l) = ne_lup*(b(ij+t_lup,l)-b(ij,l))
-         grad(ij+u_ldown,l) = ne_ldown*(b(ij+t_ldown,l)-b(ij,l))
-      END DO
-   END DO
-   !---------------------------- gradient ----------------------------------
-   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/wind_centered.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/wind_centered.k90	(revision 599)
+++ codes/icosagcm/trunk/src/kernels/wind_centered.k90	(revision 599)
@@ -0,0 +1,43 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- wind_centered ----------------------------------
+   ! Perot reconstruction based on Gauss theorem
+   ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         cx=centroid(ij,1)
+         cy=centroid(ij,2)
+         cz=centroid(ij,3)
+         ux=0. ; uy=0. ; uz=0.
+         ue_le = ne_rup*ue(ij+u_rup,l)*le(ij+u_rup)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_up,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_up,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_up,3))-cz)
+         ue_le = ne_lup*ue(ij+u_lup,l)*le(ij+u_lup)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_up,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_up,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_up,3))-cz)
+         ue_le = ne_left*ue(ij+u_left,l)*le(ij+u_left)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_ldown,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_ldown,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_ldown,3))-cz)
+         ue_le = ne_ldown*ue(ij+u_ldown,l)*le(ij+u_ldown)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_ldown,1)+xyz_v(ij+z_down,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_ldown,2)+xyz_v(ij+z_down,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_ldown,3)+xyz_v(ij+z_down,3))-cz)
+         ue_le = ne_rdown*ue(ij+u_rdown,l)*le(ij+u_rdown)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rdown,1)+xyz_v(ij+z_down,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rdown,2)+xyz_v(ij+z_down,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rdown,3)+xyz_v(ij+z_down,3))-cz)
+         ue_le = ne_right*ue(ij+u_right,l)*le(ij+u_right)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_rdown,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_rdown,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_rdown,3))-cz)
+         fac = scale*(1./Ai(ij))
+         ucenter(ij,l,1)=ux*fac
+         ucenter(ij,l,2)=uy*fac
+         ucenter(ij,l,3)=uz*fac
+      END DO
+   END DO
+   !---------------------------- wind_centered ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/physics/physics.f90
===================================================================
--- codes/icosagcm/trunk/src/physics/physics.f90	(revision 581)
+++ codes/icosagcm/trunk/src/physics/physics.f90	(revision 599)
@@ -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/trunk/src/physics/physics_lmdz_generic.F90
===================================================================
--- codes/icosagcm/trunk/src/physics/physics_lmdz_generic.F90	(revision 581)
+++ codes/icosagcm/trunk/src/physics/physics_lmdz_generic.F90	(revision 599)
@@ -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/trunk/src/time/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/time/timeloop_gcm.f90	(revision 581)
+++ codes/icosagcm/trunk/src/time/timeloop_gcm.f90	(revision 599)
@@ -168,4 +168,5 @@
     USE caldyn_mod
     USE advect_tracer_mod
+    USE diagflux_mod
     USE physics_mod
     USE mpipara
@@ -181,7 +182,7 @@
     REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:)
 
-    INTEGER :: ind
-    INTEGER :: it,i,j,l,n,  stage
-    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
+    REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out
+    INTEGER :: ind, it,i,j,l,n,  stage
+    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time
     LOGICAL, PARAMETER :: check_rhodz=.FALSE.
     INTEGER :: start_clock, stop_clock, rate_clock
@@ -210,8 +211,13 @@
 
     IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
+    IF(diagflux_on) THEN
+       adv_over_out = itau_adv*(1./itau_out)
+    ELSE
+       adv_over_out = 0.
+    END IF
 
     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
@@ -267,5 +273,5 @@
           CALL HEVI_scheme(it, fluxt_zero)
        END SELECT
-       
+
        IF (MOD(it,itau_dissip)==0) THEN
           
@@ -298,8 +304,12 @@
        
        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
-          fluxt_zero=.TRUE.
-          ! FIXME : check that rhodz is consistent with ps
-          IF (check_rhodz) THEN
+          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
+               adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt)  ! accumulate mass and fluxes if diagflux_on
+          fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
+          ! At this point advect_tracer has obtained the halos of u and rhodz,
+          ! needed for correct computation of kinetic energy
+          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt)
+
+          IF (check_rhodz) THEN ! check that rhodz is consistent with ps
              DO ind=1,ndomain
                 IF (.NOT. assigned_domain(ind)) CYCLE
Index: codes/icosagcm/trunk/src/transport/advect.F90
===================================================================
--- codes/icosagcm/trunk/src/transport/advect.F90	(revision 599)
+++ codes/icosagcm/trunk/src/transport/advect.F90	(revision 599)
@@ -0,0 +1,449 @@
+MODULE advect_mod
+  USE icosa
+  IMPLICIT NONE 
+
+  PRIVATE
+  
+  PUBLIC :: init_advect, compute_backward_traj, compute_gradq3d, compute_advect_horiz
+
+CONTAINS
+
+  !==========================================================================
+
+  SUBROUTINE init_advect(normal,tangent,sqrt_leng)
+    IMPLICIT NONE
+    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
+    REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
+    REAL(rstd),INTENT(OUT) :: sqrt_leng(iim*jjm)
+
+    INTEGER :: ij
+
+!DIR$ SIMD
+      DO ij=ij_begin,ij_end
+
+          CALL cross_product2(xyz_v(ij+z_rdown,:),xyz_v(ij+z_rup,:),normal(ij+u_right,:))
+          normal(ij+u_right,:)=normal(ij+u_right,:)/sqrt(sum(normal(ij+u_right,:)**2)+1e-50)*ne(ij,right)
+
+          CALL cross_product2(xyz_v(ij+z_up,:),xyz_v(ij+z_lup,:),normal(ij+u_lup,:))
+          normal(ij+u_lup,:)=normal(ij+u_lup,:)/sqrt(sum(normal(ij+u_lup,:)**2)+1e-50)*ne(ij,lup)
+
+          CALL cross_product2(xyz_v(ij+z_ldown,:),xyz_v(ij+z_down,:),normal(ij+u_ldown,:))	 
+          normal(ij+u_ldown,:)=normal(ij+u_ldown,:)/sqrt(sum(normal(ij+u_ldown,:)**2)+1e-50)*ne(ij,ldown)
+
+          tangent(ij+u_right,:)=xyz_v(ij+z_rup,:)-xyz_v(ij+z_rdown,:) 
+          tangent(ij+u_right,:)=tangent(ij+u_right,:)/sqrt(sum(tangent(ij+u_right,:)**2)+1e-50)
+
+          tangent(ij+u_lup,:)=xyz_v(ij+z_lup,:)-xyz_v(ij+z_up,:) 
+          tangent(ij+u_lup,:)=tangent(ij+u_lup,:)/sqrt(sum(tangent(ij+u_lup,:)**2)+1e-50)
+
+          tangent(ij+u_ldown,:)=xyz_v(ij+z_down,:)-xyz_v(ij+z_ldown,:)
+          tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50)
+
+          sqrt_leng(ij) = sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), &
+                                   sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2),  &
+                                   sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) )
+    ENDDO
+
+  END SUBROUTINE init_advect
+
+  !=======================================================================================
+
+  SUBROUTINE compute_gradq3d(qi_in,sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: qi_in(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: sqrt_leng_in(iim*jjm)
+    REAL(rstd),INTENT(IN)  :: xyz_i(iim*jjm,3)
+    REAL(rstd),INTENT(IN)  :: xyz_v(2*iim*jjm,3)
+    REAL(rstd),INTENT(OUT) :: gradq3d_out(iim*jjm,llm,3) 
+    REAL(rstd) :: maxq,minq,minq_c,maxq_c 
+    REAL(rstd) :: alphamx,alphami,alpha	,maggrd
+    REAL(rstd) :: leng1,leng2
+    REAL(rstd) :: arr(2*iim*jjm)
+    REAL(rstd) :: ar(iim*jjm)
+    REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
+    INTEGER :: ij,k,ind,l
+    REAL(rstd)  :: qi(iim*jjm,llm)
+    REAL(rstd)  :: sqrt_leng(iim*jjm)
+    REAL(rstd) :: gradq3d(iim*jjm,llm,3) 
+    REAL(rstd) :: detx,dety,detz,det
+    REAL(rstd) :: A(3,3), a11,a12,a13,a21,a22,a23,a31,a32,a33
+    REAL(rstd) :: x1,x2,x3
+    REAL(rstd) :: dq(3)
+
+    qi=qi_in
+    sqrt_leng=sqrt_leng_in
+    
+    CALL trace_start("compute_gradq3d1")
+
+    ! 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 = ll_begin,ll_end 
+!!DIR$ SIMD
+!      DO ij=ij_begin_ext,ij_end_ext
+!!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
+!!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
+!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
+!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
+!       END DO
+!    END DO
+
+     DO l = ll_begin,ll_end 
+!DIR$ SIMD
+      DO ij=ij_begin_ext,ij_end_ext
+!       CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
+
+        
+        A(1,1)=xyz_i(ij+t_rup,1)-xyz_i(ij,1);  A(1,2)=xyz_i(ij+t_rup,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_rup,3)-xyz_i(ij,3) 
+        A(2,1)=xyz_i(ij+t_lup,1)-xyz_i(ij,1);  A(2,2)=xyz_i(ij+t_lup,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_lup,3)-xyz_i(ij,3) 
+        A(3,1)=xyz_v(ij+z_up,1);               A(3,2)= xyz_v(ij+z_up,2);             A(3,3)=xyz_v(ij+z_up,3)
+    
+        dq(1) = qi(ij+t_rup,l)-qi(ij,l)
+        dq(2) = qi(ij+t_lup,l)-qi(ij,l)
+        dq(3) = 0.0
+
+
+!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         det =  x1 - x2 + x3
+                 
+!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
+
+         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         detx =  x1 - x2 + x3
+        
+!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         dety =  x1 - x2 + x3
+
+!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         detz =  x1 - x2 + x3
+
+        gradtri(ij+z_up,l,1) = detx
+        gradtri(ij+z_up,l,2) = dety
+        gradtri(ij+z_up,l,3) = detz
+        arr(ij+z_up) = det
+        
+      ENDDO
+      
+      DO ij=ij_begin_ext,ij_end_ext
+
+
+!        CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
+
+        A(1,1)=xyz_i(ij+t_ldown,1)-xyz_i(ij,1);  A(1,2)=xyz_i(ij+t_ldown,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_ldown,3)-xyz_i(ij,3) 
+        A(2,1)=xyz_i(ij+t_rdown,1)-xyz_i(ij,1);  A(2,2)=xyz_i(ij+t_rdown,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_rdown,3)-xyz_i(ij,3) 
+        A(3,1)=xyz_v(ij+z_down,1);               A(3,2)= xyz_v(ij+z_down,2);             A(3,3)=xyz_v(ij+z_down,3)
+ 
+        dq(1) = qi(ij+t_ldown,l)-qi(ij,l)
+        dq(2) = qi(ij+t_rdown,l)-qi(ij,l)
+        dq(3) = 0.0
+
+
+!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         det =  x1 - x2 + x3
+                 
+!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
+
+         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         detx =  x1 - x2 + x3
+        
+!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         dety =  x1 - x2 + x3
+
+!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         detz =  x1 - x2 + x3
+
+         gradtri(ij+z_down,l,1) = detx
+         gradtri(ij+z_down,l,2) = dety
+         gradtri(ij+z_down,l,3) = detz
+         arr(ij+z_down) = det
+
+      END DO
+    END DO
+
+!DIR$ SIMD
+    DO ij=ij_begin,ij_end
+       ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50
+    ENDDO
+    CALL trace_end("compute_gradq3d1")
+    CALL trace_start2("compute_gradq3d2")
+      
+    DO k=1,3
+      DO l =ll_begin,ll_end
+!DIR$ SIMD
+      DO ij=ij_begin,ij_end
+              gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k)   +   &
+                                 gradtri(ij+z_rup,l,k) +  gradtri(ij+z_ldown,l,k)   +   & 
+                                 gradtri(ij+z_lup,l,k)+    gradtri(ij+z_rdown,l,k) ) / ar(ij)  
+        END DO
+      END DO
+    ENDDO
+    CALL trace_end2("compute_gradq3d2")
+    CALL trace_start("compute_gradq3d3")
+
+    !============================================================================================= LIMITING 
+    DO l =ll_begin,ll_end
+!DIR$ SIMD
+      DO ij=ij_begin,ij_end
+!             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:))
+             maggrd = gradq3d(ij,l,1)*gradq3d(ij,l,1) + gradq3d(ij,l,2)*gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3) 
+             maggrd = sqrt(maggrd) 
+             maxq_c = qi(ij,l) + maggrd*sqrt_leng(ij)
+             minq_c = qi(ij,l) - maggrd*sqrt_leng(ij)
+             maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
+                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
+             minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
+                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
+             alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
+             alphamx = max(alphamx,0.0)
+             alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
+             alphami = max(alphami,0.0) 
+             alpha   = min(alphamx,alphami,1.0)
+!             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
+             gradq3d(ij,l,1) = alpha*gradq3d(ij,l,1)
+             gradq3d(ij,l,2) = alpha*gradq3d(ij,l,2)
+             gradq3d(ij,l,3) = alpha*gradq3d(ij,l,3)
+       END DO
+    END DO
+
+  CALL trace_end("compute_gradq3d3")
+  
+  
+  DO k=1,3
+    DO l = ll_begin,ll_end
+      DO ij=ij_begin,ij_end 
+        gradq3d_out(ij,l,k)=gradq3d(ij,l,k)
+      ENDDO
+    ENDDO
+  ENDDO
+  
+  CONTAINS
+
+    SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq1,dq2,dq3,det)
+    IMPLICIT NONE
+    INTEGER, INTENT(IN) :: n0,l,n1,n2,n3
+    REAL(rstd), INTENT(IN)     :: q(iim*jjm,llm)
+!    REAL(rstd), INTENT(OUT)    :: dq(3), det
+    REAL(rstd), INTENT(OUT)    :: dq1,dq2,dq3,det
+    REAL(rstd)    :: dq(3)
+    
+    REAL(rstd)  ::detx,dety,detz
+    INTEGER    :: info
+    INTEGER    :: IPIV(3)
+
+    REAL(rstd) :: A(3,3)
+    REAL(rstd) :: B(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,l)-q(n0,l)
+    dq(2) = q(n2,l)-q(n0,l)
+    dq(3) = 0.0
+
+    !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
+!    CALL determinant(A(:,1),A(:,2),A(:,3),det)
+!    CALL determinant(dq,A(:,2),A(:,3),detx)
+!    CALL determinant(A(:,1),dq,A(:,3),dety)
+!    CALL determinant(A(:,1),A(:,2),dq,detz)
+!    dq(1) = detx
+!    dq(2) = dety
+!    dq(3) = detz 
+
+    CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
+    CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),dq1)
+    CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dq2)
+    CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),dq3)
+
+  END SUBROUTINE gradq
+
+  !==========================================================================
+!  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))
+!    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
+!    det =  x1 - x2 + x3
+!  END SUBROUTINE determinant
+
+   SUBROUTINE determinant(a11,a12,a13,a21,a22,a23,a31,a32,a33,det)
+    IMPLICIT NONE
+    REAL(rstd), INTENT(IN) :: a11,a12,a13,a21,a22,a23,a31,a32,a33
+    REAL(rstd), INTENT(OUT) :: det
+    REAL(rstd) ::  x1,x2,x3
+    x1 =  a11 * (a22 * a33 - a23 * a32)
+    x2 =  a12 * (a21 * a33 - a23 * a31)
+    x3 =  a13 * (a21 * a32 - a22 * a31)
+    det =  x1 - x2 + x3
+  END SUBROUTINE determinant
+
+
+  END SUBROUTINE compute_gradq3d
+
+  ! Backward trajectories, for use with Miura approach
+  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
+    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
+    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT)   :: cc(3*iim*jjm,llm,3) ! start of backward trajectory
+    REAL(rstd),INTENT(IN)    :: tau
+
+    REAL(rstd) :: v_e(3), up_e, qe, ed(3)    
+    INTEGER :: ij,l
+
+    CALL trace_start("compute_backward_traj")
+
+    ! 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 = ll_begin,ll_end
+!DIR$ SIMD
+      DO ij=ij_begin,ij_end
+             up_e =1/de(ij+u_right)*(                                                   &
+                  wee(ij+u_right,1,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
+                  wee(ij+u_right,2,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
+                  wee(ij+u_right,3,1)*le(ij+u_left)*ue(ij+u_left,l)+                      &
+                  wee(ij+u_right,4,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                    &
+                  wee(ij+u_right,5,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    & 
+                  wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ue(ij+t_right+u_ldown,l)+    &
+                  wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ue(ij+t_right+u_rdown,l)+    &
+                  wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ue(ij+t_right+u_right,l)+    &
+                  wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ue(ij+t_right+u_rup,l)+        &
+                  wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ue(ij+t_right+u_lup,l)         &
+                  )
+             v_e = ue(ij+u_right,l)*normal(ij+u_right,:) + up_e*tangent(ij+u_right,:)
+             cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau
+
+             up_e=1/de(ij+u_lup)*(                                                      &
+                  wee(ij+u_lup,1,1)*le(ij+u_left)*ue(ij+u_left,l)+                        &
+                  wee(ij+u_lup,2,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                      &
+                  wee(ij+u_lup,3,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                      &
+                  wee(ij+u_lup,4,1)*le(ij+u_right)*ue(ij+u_right,l)+                      &
+                  wee(ij+u_lup,5,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                          & 
+                  wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ue(ij+t_lup+u_right,l)+          &
+                  wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ue(ij+t_lup+u_rup,l)+              &
+                  wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ue(ij+t_lup+u_lup,l)+              &
+                  wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ue(ij+t_lup+u_left,l)+            &
+                  wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ue(ij+t_lup+u_ldown,l)           &
+                  )
+             v_e = ue(ij+u_lup,l)*normal(ij+u_lup,:) + up_e*tangent(ij+u_lup,:)
+             cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau
+            
+
+             up_e=1/de(ij+u_ldown)*(                                                    &
+                  wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    &
+                  wee(ij+u_ldown,2,1)*le(ij+u_right)*ue(ij+u_right,l)+                    &
+                  wee(ij+u_ldown,3,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
+                  wee(ij+u_ldown,4,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
+                  wee(ij+u_ldown,5,1)*le(ij+u_left)*ue(ij+u_left,l)+                      & 
+                  wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ue(ij+t_ldown+u_lup,l)+        &
+                  wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ue(ij+t_ldown+u_left,l)+      &
+                  wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ue(ij+t_ldown+u_ldown,l)+    &
+                  wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ue(ij+t_ldown+u_rdown,l)+    &
+                  wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ue(ij+t_ldown+u_right,l)     &
+                  )
+             v_e = ue(ij+u_ldown,l)*normal(ij+u_ldown,:) + up_e*tangent(ij+u_ldown,:)
+             cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau
+       ENDDO
+    END DO
+
+    CALL trace_end("compute_backward_traj")
+
+  END SUBROUTINE compute_backward_traj
+
+  ! Horizontal transport (S. Dubey, T. Dubos)
+  ! Slope-limited van Leer approach with hexagons
+  SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass,qi,qfluxt)
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    LOGICAL, INTENT(IN)       :: update_mass, diagflux_on
+    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), INTENT(INOUT) :: qfluxt(3*iim*jjm,llm) ! time-integrated tracer flux
+
+    REAL(rstd) :: dq,dmass,qe,ed(3), newmass
+    REAL(rstd) :: qflux(3*iim*jjm,llm)
+    INTEGER :: ij,k,l
+
+    CALL trace_start("compute_advect_horiz")
+#include "../kernels/advect_horiz.k90"
+    CALL trace_end("compute_advect_horiz")
+  END SUBROUTINE compute_advect_horiz
+
+
+END MODULE advect_mod
Index: codes/icosagcm/trunk/src/transport/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/transport/advect.f90	(revision 581)
+++ 	(revision )
@@ -1,528 +1,0 @@
-MODULE advect_mod
-  USE icosa
-  IMPLICIT NONE 
-
-  PRIVATE
-  
-  PUBLIC :: init_advect, compute_backward_traj, compute_gradq3d, compute_advect_horiz
-
-CONTAINS
-
-  !==========================================================================
-
-  SUBROUTINE init_advect(normal,tangent,sqrt_leng)
-    IMPLICIT NONE
-    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
-    REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
-    REAL(rstd),INTENT(OUT) :: sqrt_leng(iim*jjm)
-
-    INTEGER :: ij
-
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-
-          CALL cross_product2(xyz_v(ij+z_rdown,:),xyz_v(ij+z_rup,:),normal(ij+u_right,:))
-          normal(ij+u_right,:)=normal(ij+u_right,:)/sqrt(sum(normal(ij+u_right,:)**2)+1e-50)*ne(ij,right)
-
-          CALL cross_product2(xyz_v(ij+z_up,:),xyz_v(ij+z_lup,:),normal(ij+u_lup,:))
-          normal(ij+u_lup,:)=normal(ij+u_lup,:)/sqrt(sum(normal(ij+u_lup,:)**2)+1e-50)*ne(ij,lup)
-
-          CALL cross_product2(xyz_v(ij+z_ldown,:),xyz_v(ij+z_down,:),normal(ij+u_ldown,:))	 
-          normal(ij+u_ldown,:)=normal(ij+u_ldown,:)/sqrt(sum(normal(ij+u_ldown,:)**2)+1e-50)*ne(ij,ldown)
-
-          tangent(ij+u_right,:)=xyz_v(ij+z_rup,:)-xyz_v(ij+z_rdown,:) 
-          tangent(ij+u_right,:)=tangent(ij+u_right,:)/sqrt(sum(tangent(ij+u_right,:)**2)+1e-50)
-
-          tangent(ij+u_lup,:)=xyz_v(ij+z_lup,:)-xyz_v(ij+z_up,:) 
-          tangent(ij+u_lup,:)=tangent(ij+u_lup,:)/sqrt(sum(tangent(ij+u_lup,:)**2)+1e-50)
-
-          tangent(ij+u_ldown,:)=xyz_v(ij+z_down,:)-xyz_v(ij+z_ldown,:)
-          tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50)
-
-          sqrt_leng(ij) = sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), &
-                                   sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2),  &
-                                   sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) )
-    ENDDO
-
-  END SUBROUTINE init_advect
-
-  !=======================================================================================
-
-  SUBROUTINE compute_gradq3d(qi_in,sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: qi_in(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: sqrt_leng_in(iim*jjm)
-    REAL(rstd),INTENT(IN)  :: xyz_i(iim*jjm,3)
-    REAL(rstd),INTENT(IN)  :: xyz_v(2*iim*jjm,3)
-    REAL(rstd),INTENT(OUT) :: gradq3d_out(iim*jjm,llm,3) 
-    REAL(rstd) :: maxq,minq,minq_c,maxq_c 
-    REAL(rstd) :: alphamx,alphami,alpha	,maggrd
-    REAL(rstd) :: leng1,leng2
-    REAL(rstd) :: arr(2*iim*jjm)
-    REAL(rstd) :: ar(iim*jjm)
-    REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
-    INTEGER :: ij,k,ind,l
-    REAL(rstd)  :: qi(iim*jjm,llm)
-    REAL(rstd)  :: sqrt_leng(iim*jjm)
-    REAL(rstd) :: gradq3d(iim*jjm,llm,3) 
-    REAL(rstd) :: detx,dety,detz,det
-    REAL(rstd) :: A(3,3), a11,a12,a13,a21,a22,a23,a31,a32,a33
-    REAL(rstd) :: x1,x2,x3
-    REAL(rstd) :: dq(3)
-
-    qi=qi_in
-    sqrt_leng=sqrt_leng_in
-    
-    CALL trace_start("compute_gradq3d1")
-
-    ! 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 = ll_begin,ll_end 
-!!DIR$ SIMD
-!      DO ij=ij_begin_ext,ij_end_ext
-!!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
-!!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
-!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
-!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
-!       END DO
-!    END DO
-
-     DO l = ll_begin,ll_end 
-!DIR$ SIMD
-      DO ij=ij_begin_ext,ij_end_ext
-!       CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
-
-        
-        A(1,1)=xyz_i(ij+t_rup,1)-xyz_i(ij,1);  A(1,2)=xyz_i(ij+t_rup,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_rup,3)-xyz_i(ij,3) 
-        A(2,1)=xyz_i(ij+t_lup,1)-xyz_i(ij,1);  A(2,2)=xyz_i(ij+t_lup,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_lup,3)-xyz_i(ij,3) 
-        A(3,1)=xyz_v(ij+z_up,1);               A(3,2)= xyz_v(ij+z_up,2);             A(3,3)=xyz_v(ij+z_up,3)
-    
-        dq(1) = qi(ij+t_rup,l)-qi(ij,l)
-        dq(2) = qi(ij+t_lup,l)-qi(ij,l)
-        dq(3) = 0.0
-
-
-!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
-
-         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
-         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
-         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
-
-         x1 =  a11 * (a22 * a33 - a23 * a32)
-         x2 =  a12 * (a21 * a33 - a23 * a31)
-         x3 =  a13 * (a21 * a32 - a22 * a31)
-         det =  x1 - x2 + x3
-                 
-!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
-
-         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3)
-         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
-         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
-
-         x1 =  a11 * (a22 * a33 - a23 * a32)
-         x2 =  a12 * (a21 * a33 - a23 * a31)
-         x3 =  a13 * (a21 * a32 - a22 * a31)
-         detx =  x1 - x2 + x3
-        
-!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
-
-         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
-         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3)
-         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
-
-         x1 =  a11 * (a22 * a33 - a23 * a32)
-         x2 =  a12 * (a21 * a33 - a23 * a31)
-         x3 =  a13 * (a21 * a32 - a22 * a31)
-         dety =  x1 - x2 + x3
-
-!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
-
-         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
-         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
-         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3)
-
-         x1 =  a11 * (a22 * a33 - a23 * a32)
-         x2 =  a12 * (a21 * a33 - a23 * a31)
-         x3 =  a13 * (a21 * a32 - a22 * a31)
-         detz =  x1 - x2 + x3
-
-        gradtri(ij+z_up,l,1) = detx
-        gradtri(ij+z_up,l,2) = dety
-        gradtri(ij+z_up,l,3) = detz
-        arr(ij+z_up) = det
-        
-      ENDDO
-      
-      DO ij=ij_begin_ext,ij_end_ext
-
-
-!        CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
-
-        A(1,1)=xyz_i(ij+t_ldown,1)-xyz_i(ij,1);  A(1,2)=xyz_i(ij+t_ldown,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_ldown,3)-xyz_i(ij,3) 
-        A(2,1)=xyz_i(ij+t_rdown,1)-xyz_i(ij,1);  A(2,2)=xyz_i(ij+t_rdown,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_rdown,3)-xyz_i(ij,3) 
-        A(3,1)=xyz_v(ij+z_down,1);               A(3,2)= xyz_v(ij+z_down,2);             A(3,3)=xyz_v(ij+z_down,3)
- 
-        dq(1) = qi(ij+t_ldown,l)-qi(ij,l)
-        dq(2) = qi(ij+t_rdown,l)-qi(ij,l)
-        dq(3) = 0.0
-
-
-!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
-
-         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
-         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
-         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
-
-         x1 =  a11 * (a22 * a33 - a23 * a32)
-         x2 =  a12 * (a21 * a33 - a23 * a31)
-         x3 =  a13 * (a21 * a32 - a22 * a31)
-         det =  x1 - x2 + x3
-                 
-!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
-
-         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3)
-         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
-         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
-
-         x1 =  a11 * (a22 * a33 - a23 * a32)
-         x2 =  a12 * (a21 * a33 - a23 * a31)
-         x3 =  a13 * (a21 * a32 - a22 * a31)
-         detx =  x1 - x2 + x3
-        
-!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
-
-         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
-         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3)
-         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
-
-         x1 =  a11 * (a22 * a33 - a23 * a32)
-         x2 =  a12 * (a21 * a33 - a23 * a31)
-         x3 =  a13 * (a21 * a32 - a22 * a31)
-         dety =  x1 - x2 + x3
-
-!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
-
-         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
-         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
-         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3)
-
-         x1 =  a11 * (a22 * a33 - a23 * a32)
-         x2 =  a12 * (a21 * a33 - a23 * a31)
-         x3 =  a13 * (a21 * a32 - a22 * a31)
-         detz =  x1 - x2 + x3
-
-         gradtri(ij+z_down,l,1) = detx
-         gradtri(ij+z_down,l,2) = dety
-         gradtri(ij+z_down,l,3) = detz
-         arr(ij+z_down) = det
-
-      END DO
-    END DO
-
-!DIR$ SIMD
-    DO ij=ij_begin,ij_end
-       ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50
-    ENDDO
-    CALL trace_end("compute_gradq3d1")
-    CALL trace_start2("compute_gradq3d2")
-      
-    DO k=1,3
-      DO l =ll_begin,ll_end
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-              gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k)   +   &
-                                 gradtri(ij+z_rup,l,k) +  gradtri(ij+z_ldown,l,k)   +   & 
-                                 gradtri(ij+z_lup,l,k)+    gradtri(ij+z_rdown,l,k) ) / ar(ij)  
-        END DO
-      END DO
-    ENDDO
-    CALL trace_end2("compute_gradq3d2")
-    CALL trace_start("compute_gradq3d3")
-
-    !============================================================================================= LIMITING 
-    DO l =ll_begin,ll_end
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-!             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:))
-             maggrd = gradq3d(ij,l,1)*gradq3d(ij,l,1) + gradq3d(ij,l,2)*gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3) 
-             maggrd = sqrt(maggrd) 
-             maxq_c = qi(ij,l) + maggrd*sqrt_leng(ij)
-             minq_c = qi(ij,l) - maggrd*sqrt_leng(ij)
-             maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
-                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
-             minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
-                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
-             alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
-             alphamx = max(alphamx,0.0)
-             alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
-             alphami = max(alphami,0.0) 
-             alpha   = min(alphamx,alphami,1.0)
-!             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
-             gradq3d(ij,l,1) = alpha*gradq3d(ij,l,1)
-             gradq3d(ij,l,2) = alpha*gradq3d(ij,l,2)
-             gradq3d(ij,l,3) = alpha*gradq3d(ij,l,3)
-       END DO
-    END DO
-
-  CALL trace_end("compute_gradq3d3")
-  
-  
-  DO k=1,3
-    DO l = ll_begin,ll_end
-      DO ij=ij_begin,ij_end 
-        gradq3d_out(ij,l,k)=gradq3d(ij,l,k)
-      ENDDO
-    ENDDO
-  ENDDO
-  
-  CONTAINS
-
-    SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq1,dq2,dq3,det)
-    IMPLICIT NONE
-    INTEGER, INTENT(IN) :: n0,l,n1,n2,n3
-    REAL(rstd), INTENT(IN)     :: q(iim*jjm,llm)
-!    REAL(rstd), INTENT(OUT)    :: dq(3), det
-    REAL(rstd), INTENT(OUT)    :: dq1,dq2,dq3,det
-    REAL(rstd)    :: dq(3)
-    
-    REAL(rstd)  ::detx,dety,detz
-    INTEGER    :: info
-    INTEGER    :: IPIV(3)
-
-    REAL(rstd) :: A(3,3)
-    REAL(rstd) :: B(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,l)-q(n0,l)
-    dq(2) = q(n2,l)-q(n0,l)
-    dq(3) = 0.0
-
-    !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
-!    CALL determinant(A(:,1),A(:,2),A(:,3),det)
-!    CALL determinant(dq,A(:,2),A(:,3),detx)
-!    CALL determinant(A(:,1),dq,A(:,3),dety)
-!    CALL determinant(A(:,1),A(:,2),dq,detz)
-!    dq(1) = detx
-!    dq(2) = dety
-!    dq(3) = detz 
-
-    CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
-    CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),dq1)
-    CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dq2)
-    CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),dq3)
-
-  END SUBROUTINE gradq
-
-  !==========================================================================
-!  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))
-!    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
-!    det =  x1 - x2 + x3
-!  END SUBROUTINE determinant
-
-   SUBROUTINE determinant(a11,a12,a13,a21,a22,a23,a31,a32,a33,det)
-    IMPLICIT NONE
-    REAL(rstd), INTENT(IN) :: a11,a12,a13,a21,a22,a23,a31,a32,a33
-    REAL(rstd), INTENT(OUT) :: det
-    REAL(rstd) ::  x1,x2,x3
-    x1 =  a11 * (a22 * a33 - a23 * a32)
-    x2 =  a12 * (a21 * a33 - a23 * a31)
-    x3 =  a13 * (a21 * a32 - a22 * a31)
-    det =  x1 - x2 + x3
-  END SUBROUTINE determinant
-
-
-  END SUBROUTINE compute_gradq3d
-
-  ! Backward trajectories, for use with Miura approach
-  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
-    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
-    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
-    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
-    REAL(rstd),INTENT(OUT)   :: cc(3*iim*jjm,llm,3) ! start of backward trajectory
-    REAL(rstd),INTENT(IN)    :: tau
-
-    REAL(rstd) :: v_e(3), up_e, qe, ed(3)    
-    INTEGER :: ij,l
-
-    CALL trace_start("compute_backward_traj")
-
-    ! 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 = ll_begin,ll_end
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-             up_e =1/de(ij+u_right)*(                                                   &
-                  wee(ij+u_right,1,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
-                  wee(ij+u_right,2,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
-                  wee(ij+u_right,3,1)*le(ij+u_left)*ue(ij+u_left,l)+                      &
-                  wee(ij+u_right,4,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                    &
-                  wee(ij+u_right,5,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    & 
-                  wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ue(ij+t_right+u_ldown,l)+    &
-                  wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ue(ij+t_right+u_rdown,l)+    &
-                  wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ue(ij+t_right+u_right,l)+    &
-                  wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ue(ij+t_right+u_rup,l)+        &
-                  wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ue(ij+t_right+u_lup,l)         &
-                  )
-             v_e = ue(ij+u_right,l)*normal(ij+u_right,:) + up_e*tangent(ij+u_right,:)
-             cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau
-
-             up_e=1/de(ij+u_lup)*(                                                      &
-                  wee(ij+u_lup,1,1)*le(ij+u_left)*ue(ij+u_left,l)+                        &
-                  wee(ij+u_lup,2,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                      &
-                  wee(ij+u_lup,3,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                      &
-                  wee(ij+u_lup,4,1)*le(ij+u_right)*ue(ij+u_right,l)+                      &
-                  wee(ij+u_lup,5,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                          & 
-                  wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ue(ij+t_lup+u_right,l)+          &
-                  wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ue(ij+t_lup+u_rup,l)+              &
-                  wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ue(ij+t_lup+u_lup,l)+              &
-                  wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ue(ij+t_lup+u_left,l)+            &
-                  wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ue(ij+t_lup+u_ldown,l)           &
-                  )
-             v_e = ue(ij+u_lup,l)*normal(ij+u_lup,:) + up_e*tangent(ij+u_lup,:)
-             cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau
-            
-
-             up_e=1/de(ij+u_ldown)*(                                                    &
-                  wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    &
-                  wee(ij+u_ldown,2,1)*le(ij+u_right)*ue(ij+u_right,l)+                    &
-                  wee(ij+u_ldown,3,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
-                  wee(ij+u_ldown,4,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
-                  wee(ij+u_ldown,5,1)*le(ij+u_left)*ue(ij+u_left,l)+                      & 
-                  wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ue(ij+t_ldown+u_lup,l)+        &
-                  wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ue(ij+t_ldown+u_left,l)+      &
-                  wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ue(ij+t_ldown+u_ldown,l)+    &
-                  wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ue(ij+t_ldown+u_rdown,l)+    &
-                  wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ue(ij+t_ldown+u_right,l)     &
-                  )
-             v_e = ue(ij+u_ldown,l)*normal(ij+u_ldown,:) + up_e*tangent(ij+u_ldown,:)
-             cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau
-       ENDDO
-    END DO
-
-    CALL trace_end("compute_backward_traj")
-
-  END SUBROUTINE compute_backward_traj
-
-  ! Horizontal transport (S. Dubey, T. Dubos)
-  ! Slope-limited van Leer approach with hexagons
-  SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 
-    USE trace
-    USE omp_para
-    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 :: ij,k,l
-
-    CALL trace_start("compute_advect_horiz")
-
-    ! 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 = ll_begin,ll_end
-
-!DIR$ SIMD
-      DO ij=ij_begin_ext,ij_end_ext
-
-             IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN 
-                ed = cc(ij+u_right,l,:) - xyz_i(ij,:)
-!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
-                 qe = qi(ij,l)+gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
-             ELSE
-                ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:)
-!                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 
-                qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)*ed(1)+gradq3d(ij+t_right,l,2)*ed(2)+gradq3d(ij+t_right,l,3)*ed(3) 
-             ENDIF
-             qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
-             
-             IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN 
-                ed = cc(ij+u_lup,l,:) - xyz_i(ij,:)
-!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
-                qe = qi(ij,l)  + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
-             ELSE
-                ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:)
-!                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 
-                qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)*ed(1)+gradq3d(ij+t_lup,l,2)*ed(2)+gradq3d(ij+t_lup,l,3)*ed(3) 
-             ENDIF
-             qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe 
-
-             IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN 
-                ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:)
-!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
-                qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
-             ELSE
-                ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:)
-!                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 
-                qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)*ed(1)+gradq3d(ij+t_ldown,l,2)*ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3) 
-             ENDIF
-             qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
-       END DO
-    END DO
-
-    ! update q and, if update_mass, update mass
-    DO l = ll_begin,ll_end
-!DIR$ SIMD
-      DO ij=ij_begin,ij_end
-             ! sign convention as Ringler et al. (2010) eq. 21
-             dmass =   hfluxt(ij+u_right,l)*ne(ij,right)   &
-                     + hfluxt(ij+u_lup,l)  *ne(ij,lup)     &
-                     + hfluxt(ij+u_ldown,l)*ne(ij,ldown)   &
-                     + hfluxt(ij+u_rup,l)  *ne(ij,rup)     &
-                     + hfluxt(ij+u_left,l) *ne(ij,left)    &
-                     + hfluxt(ij+u_rdown,l)*ne(ij,rdown)
-
-             dq    =   qflux(ij+u_right,l) *ne(ij,right)   &
-                     + qflux(ij+u_lup,l)   *ne(ij,lup)     &
-                     + qflux(ij+u_ldown,l) *ne(ij,ldown)   &
-                     + qflux(ij+u_rup,l)   *ne(ij,rup)     &
-                     + qflux(ij+u_left,l)  *ne(ij,left)    &
-                     + qflux(ij+u_rdown,l) *ne(ij,rdown)
-
-             
-             newmass = mass(ij,l) - dmass/Ai(ij)
-             qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass
-             IF(update_mass) mass(ij,l) = newmass
-
-       END DO
-    END DO
-
-    CALL trace_end("compute_advect_horiz")
-  CONTAINS
-
-    !====================================================================================
-    PURE REAL(rstd) FUNCTION sum2(a1,a2)
-      IMPLICIT NONE
-      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
-
-
-END MODULE advect_mod
Index: codes/icosagcm/trunk/src/transport/advect_tracer.f90
===================================================================
--- codes/icosagcm/trunk/src/transport/advect_tracer.f90	(revision 581)
+++ codes/icosagcm/trunk/src/transport/advect_tracer.f90	(revision 599)
@@ -1,4 +1,5 @@
 MODULE advect_tracer_mod
   USE icosa
+  USE advect_mod
   IMPLICIT NONE
   PRIVATE
@@ -25,5 +26,4 @@
 
   SUBROUTINE init_advect_tracer
-    USE advect_mod
     USE omp_para
     REAL(rstd),POINTER :: tangent(:,:)
@@ -54,12 +54,10 @@
   END SUBROUTINE init_advect_tracer
 
-  SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz)
-    USE advect_mod
-    USE mpipara
+  SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz,&
+       frac, f_masst,f_qmasst,f_massfluxt,f_qfluxt)
+    USE omp_para
     USE trace
     USE write_field_mod
     USE tracer_mod
-    IMPLICIT NONE
-    
     TYPE(t_field),POINTER :: f_hfluxt(:)   ! time-integrated horizontal mass flux
     TYPE(t_field),POINTER :: f_wfluxt(:)   ! time-integrated vertical mass flux
@@ -67,7 +65,12 @@
     TYPE(t_field),POINTER :: f_q(:)        ! tracer
     TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step
+    REAL(rstd), INTENT(in):: frac          ! ratio itau_adv/itau_out or 0. if not diagflux_on
+    TYPE(t_field),POINTER :: f_masst(:)    ! time-integrated mass
+    TYPE(t_field),POINTER :: f_qmasst(:)   ! time-integrated tracer mass
+    TYPE(t_field),POINTER :: f_massfluxt(:)! time-integrated horizontal mass flux
+    TYPE(t_field),POINTER :: f_qfluxt(:)   ! time-integrated horizontal tracer flux
 
     REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:)
-    REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:)
+    REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:), masst(:,:), qmasst(:,:,:), massfluxt(:,:), qfluxt(:,:,:)
     REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 
 ! temporary shared variable for vlz
@@ -77,5 +80,5 @@
     REAL(rstd),POINTER ::  wq(:,:)           ! time-integrated flux of q
     
-     INTEGER :: ind,k, nq_last
+    INTEGER :: ind,k, nq_last
     LOGICAL,SAVE :: first=.TRUE.
 !$OMP THREADPRIVATE(first)
@@ -141,5 +144,4 @@
     CALL send_message(f_cc,req_cc)
 
-
     ! horizontal transport - split in two to place transfer of gradq3d
     DO k = 1, nqtot
@@ -154,5 +156,4 @@
           sqrt_leng=f_sqrt_leng(ind)
           CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v)
-
         END DO
 
@@ -160,5 +161,4 @@
         CALL wait_message(req_cc)
         CALL wait_message(req_gradq3d)
-
 
         DO ind=1,ndomain
@@ -170,7 +170,21 @@
           rhodz   = f_rhodz(ind)
           hfluxt  = f_hfluxt(ind) 
+          qfluxt  = f_qfluxt(ind) 
           gradq3d = f_gradq3d(ind)
-          CALL compute_advect_horiz(k==nq_last,hfluxt,cc,gradq3d, rhodz,q(:,:,k))
+
+          IF(frac>0.) THEN ! accumulate mass, mass flux and tracer mass
+             qmasst  = f_qmasst(ind) 
+             qmasst(:,ll_begin:ll_end,k) = qmasst(:,ll_begin:ll_end,k) + &
+                  frac*rhodz(:,ll_begin:ll_end)*q(:,ll_begin:ll_end,k)
+             IF(k==nq_last) THEN
+                masst  = f_masst(ind)
+                massfluxt  = f_massfluxt(ind)
+                masst(:,ll_begin:ll_end) = masst(:,ll_begin:ll_end)+frac*rhodz(:,ll_begin:ll_end)
+                massfluxt(:,ll_begin:ll_end) = massfluxt(:,ll_begin:ll_end)+hfluxt(:,ll_begin:ll_end)
+             END IF
+          END IF
+          CALL compute_advect_horiz(k==nq_last,frac>0., hfluxt,cc,gradq3d, rhodz, q(:,:,k), qfluxt(:,:,k))
         END DO
+
       ENDIF
     END DO 
@@ -214,5 +228,4 @@
     USE trace
     USE omp_para
-    IMPLICIT NONE
     LOGICAL, INTENT(IN)       :: update_mass
     REAL(rstd), INTENT(IN)    :: fac, wfluxt(iim*jjm,llm+1) ! vertical mass flux
