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 593)
@@ -0,0 +1,90 @@
+MODULE diagflux_mod
+  USE icosa
+  IMPLICIT NONE
+  SAVE
+  
+  TYPE(t_field),POINTER :: &
+       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)
+
+CONTAINS
+
+  SUBROUTINE init_diagflux
+    USE getin_mod
+    INTEGER :: ll
+    diagflux_on = .FALSE.
+    CALL getin("diagflux", diagflux_on)
+    ll = MERGE(ll,0,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
+
+  SUBROUTINE zero_qfluxt
+    USE omp_para
+    INTEGER :: ind
+    REAL(rstd), POINTER :: buf2(:,:),buf3(:,:,:)
+    DO ind=1,ndomain
+       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.
+    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)
+    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(:,:)
+    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
+       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 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
+
+END MODULE diagflux_mod
Index: codes/icosagcm/devel/src/diagnostics/diagflux.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/diagflux.f90	(revision 592)
+++ 	(revision )
@@ -1,90 +1,0 @@
-MODULE diagflux_mod
-  USE icosa
-  IMPLICIT NONE
-  SAVE
-  
-  TYPE(t_field),POINTER :: &
-       f_masst(:), f_qmasst(:), & ! time-integrated mass, tracer mass, 
-       f_massfluxt(:), f_qfluxt(:), & ! mass flux and tracer flux
-       f_qfluxt_lon(:), f_qfluxt_lat(:) ! scalar flux reconstructed  cell centers
-  LOGICAL :: diagflux_on
-  !$OMP THREADPRIVATE(diagflux_on)
-
-CONTAINS
-
-  SUBROUTINE init_diagflux
-    USE getin_mod
-    diagflux_on = .FALSE.
-    CALL getin("diagflux", diagflux_on)
-    IF(diagflux_on) THEN
-       CALL allocate_field(f_masst,  field_t,type_real,llm, name="masst")
-       CALL allocate_field(f_qmasst,  field_t,type_real,llm,nqtot, name="qmasst")
-       CALL allocate_field(f_massfluxt, field_u,type_real,llm, name="massfluxt")
-       CALL allocate_field(f_qfluxt,  field_u,type_real,llm,nqtot, name="qfluxt")
-       CALL allocate_field(f_qfluxt_lon,  field_t,type_real,llm,nqtot, name="qfluxt_lon")
-       CALL allocate_field(f_qfluxt_lat,  field_t,type_real,llm,nqtot, name="qfluxt_lat")
-       CALL zero_qfluxt
-    ELSE
-       CALL allocate_field(f_masst,  field_t,type_real,0, name="masst")
-       CALL allocate_field(f_qmasst,  field_t,type_real,llm,0, name="qmasst")
-       CALL allocate_field(f_massfluxt, field_u,type_real,0, name="massfluxt")
-       CALL allocate_field(f_qfluxt,  field_u,type_real,llm,0, name="qfluxt")
-       CALL allocate_field(f_qfluxt_lon,  field_t,type_real,llm,0, name="qfluxt_lon")
-       CALL allocate_field(f_qfluxt_lat,  field_t,type_real,llm,0, name="qfluxt_lat")
-    END IF
-  END SUBROUTINE init_diagflux
-
-  SUBROUTINE zero_qfluxt
-    USE mpipara
-    USE omp_para
-    INTEGER :: ind
-    REAL(rstd), POINTER :: buf2(:,:),buf3(:,:,:)
-    DO ind=1,ndomain
-       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.
-    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)
-    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(:,:)
-    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
-       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 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
-
-END MODULE diagflux_mod
Index: codes/icosagcm/devel/src/kernels/energetics.k90
===================================================================
--- codes/icosagcm/devel/src/kernels/energetics.k90	(revision 593)
+++ codes/icosagcm/devel/src/kernels/energetics.k90	(revision 593)
@@ -0,0 +1,61 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- energetics ----------------------------------
+   ! potential energy
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         energy = .5*(geopot(ij,l+1)+geopot(ij,l))
+         epot(ij,l) = epot(ij,l) + frac*mass(ij,l)*energy
+         ebuf(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         epot_flux(ij+u_right,l) = epot_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(ebuf(ij,l)+ebuf(ij+t_right,l))
+         epot_flux(ij+u_lup,l) = epot_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(ebuf(ij,l)+ebuf(ij+t_lup,l))
+         epot_flux(ij+u_ldown,l) = epot_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(ebuf(ij,l)+ebuf(ij+t_ldown,l))
+      END DO
+   END DO
+   ! enthalpy
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         energy = cpp*temp(ij,l)
+         enthalpy(ij,l) = enthalpy(ij,l) + frac*mass(ij,l)*energy
+         ebuf(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         enthalpy_flux(ij+u_right,l) = enthalpy_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(ebuf(ij,l)+ebuf(ij+t_right,l))
+         enthalpy_flux(ij+u_lup,l) = enthalpy_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(ebuf(ij,l)+ebuf(ij+t_lup,l))
+         enthalpy_flux(ij+u_ldown,l) = enthalpy_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(ebuf(ij,l)+ebuf(ij+t_ldown,l))
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         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*mass(ij,l)*energy
+         ebuf(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         ekin_flux(ij+u_right,l) = ekin_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(ebuf(ij,l)+ebuf(ij+t_right,l))
+         ekin_flux(ij+u_lup,l) = ekin_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(ebuf(ij,l)+ebuf(ij+t_lup,l))
+         ekin_flux(ij+u_ldown,l) = ekin_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(ebuf(ij,l)+ebuf(ij+t_ldown,l))
+      END DO
+   END DO
+   !---------------------------- energetics ----------------------------------
+   !--------------------------------------------------------------------------
