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