Index: codes/icosagcm/devel/src/diagnostics/diagflux.F90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/diagflux.F90	(revision 593)
+++ codes/icosagcm/devel/src/diagnostics/diagflux.F90	(revision 594)
@@ -1,8 +1,10 @@
 MODULE diagflux_mod
   USE icosa
+  USE omp_para
   IMPLICIT NONE
   SAVE
-  
-  TYPE(t_field),POINTER :: &
+  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
@@ -13,4 +15,6 @@
   !$OMP THREADPRIVATE(diagflux_on)
 
+  PUBLIC :: diagflux_on, init_diagflux, zero_qfluxt, qflux_centered_lonlat, flux_centered_lonlat, diagflux_energy
+
 CONTAINS
 
@@ -20,5 +24,5 @@
     diagflux_on = .FALSE.
     CALL getin("diagflux", diagflux_on)
-    ll = MERGE(ll,0,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")
@@ -36,6 +40,8 @@
   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
-    USE omp_para
     INTEGER :: ind
     REAL(rstd), POINTER :: buf2(:,:),buf3(:,:,:)
@@ -43,21 +49,23 @@
        IF (.NOT. assigned_domain(ind)) CYCLE
        CALL swap_dimensions(ind)
-       buf2=f_masst(ind)
-       buf2(:,ll_begin:ll_end)=0.
-       buf2=f_massfluxt(ind)
-       buf2(:,ll_begin:ll_end)=0.
-       buf3=f_qmasst(ind)
-       buf3(:,ll_begin:ll_end,:)=0.
-       buf3=f_qfluxt(ind)
-       buf3(:,ll_begin:ll_end,:)=0.
+       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
 
-  SUBROUTINE flux_centered_lonlat(scale, f_massflux, f_flux, f_massflux_lon, f_massflux_lat, f_flux_lon, f_flux_lat)
+!------------------------------------ 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(:), &
-         f_massflux(:), f_massflux_lon(:), f_massflux_lat(:)
-    REAL(rstd), POINTER :: flux(:,:,:), flux_lon(:,:,:), flux_lat(:,:,:), &
-         massflux(:,:), massflux_lon(:,:), massflux_lat(:,:)
+    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
@@ -71,8 +79,20 @@
           CALL compute_flux_centered_lonlat(scale, flux(:,:,itrac), flux_lon(:,:,itrac), flux_lat(:,:,itrac))
        END DO
-       massflux=f_massflux(ind)
-       massflux_lon=f_massflux_lon(ind)
-       massflux_lat=f_massflux_lat(ind)
-       CALL compute_flux_centered_lonlat(scale, massflux, massflux_lon, massflux_lat)
+    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
@@ -88,3 +108,51 @@
   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
+    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
+    ij_omp_begin_ext = ij_omp_begin_ext+ij_begin_ext-1
+    ij_omp_end_ext = ij_omp_end_ext+ij_begin_ext-1
+#include "../kernels/energy_fluxes.k90"
+  END SUBROUTINE compute_diagflux_energy
+
 END MODULE diagflux_mod
Index: codes/icosagcm/devel/src/diagnostics/observable.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/observable.f90	(revision 593)
+++ codes/icosagcm/devel/src/diagnostics/observable.f90	(revision 594)
@@ -163,10 +163,22 @@
     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("qmass_t", f_qmasst)
-          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_qfluxt, &
-               f_buf_ulon, f_buf_ulat, f_qfluxt_lon, f_qfluxt_lat)
           CALL output_field("massflux_lon",f_buf_ulon)
           CALL output_field("massflux_lat",f_buf_ulat)
+          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 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 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)
Index: codes/icosagcm/devel/src/kernels/energy_fluxes.k90
===================================================================
--- codes/icosagcm/devel/src/kernels/energy_fluxes.k90	(revision 594)
+++ codes/icosagcm/devel/src/kernels/energy_fluxes.k90	(revision 594)
@@ -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/devel/src/time/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 593)
+++ codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 594)
@@ -273,5 +273,5 @@
           CALL HEVI_scheme(it, fluxt_zero)
        END SELECT
-       
+
        IF (MOD(it,itau_dissip)==0) THEN
           
@@ -304,9 +304,9 @@
        
        IF(MOD(it,itau_adv)==0) THEN
+          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt)
           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.
-          ! FIXME : check that rhodz is consistent with ps
-          IF (check_rhodz) THEN
+          fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
+          IF (check_rhodz) THEN ! check that rhodz is consistent with ps
              DO ind=1,ndomain
                 IF (.NOT. assigned_domain(ind)) CYCLE
Index: codes/icosagcm/devel/src/transport/advect.F90
===================================================================
--- codes/icosagcm/devel/src/transport/advect.F90	(revision 593)
+++ codes/icosagcm/devel/src/transport/advect.F90	(revision 594)
@@ -444,14 +444,4 @@
 #include "../kernels/advect_horiz.k90"
     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
 
