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 ----------------------------------
+   !--------------------------------------------------------------------------
