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)
