!--------------------------------------------------------------------------
   !---------------------------- 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
   ! NB : at this point pressure is stored in array pk
   ! pk then serves as buffer to store temperature
   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)
            theta(ij,l) = theta_ik
            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)
            theta(ij,l) = Treff*exp(theta_ik/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
   ! NB : at this point temperature is stored in array pk
   ! pk then serves as buffer to store energy
   ! enthalpy
   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
   ! theta
   DO l = ll_begin, ll_end
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         energy = theta(ij,l)
         thetat(ij,l) = thetat(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
         thetat_flux(ij+u_right,l) = thetat_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
         thetat_flux(ij+u_lup,l) = thetat_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
         thetat_flux(ij+u_ldown,l) = thetat_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)*ue(ij+u_rup,l)**2
         energy = energy + le(ij+u_lup)*de(ij+u_lup)*ue(ij+u_lup,l)**2
         energy = energy + le(ij+u_left)*de(ij+u_left)*ue(ij+u_left,l)**2
         energy = energy + le(ij+u_ldown)*de(ij+u_ldown)*ue(ij+u_ldown,l)**2
         energy = energy + le(ij+u_rdown)*de(ij+u_rdown)*ue(ij+u_rdown,l)**2
         energy = energy + le(ij+u_right)*de(ij+u_right)*ue(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
   ! ulon
   DO l = ll_begin, ll_end
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         cx=centroid(ij,1)
         cy=centroid(ij,2)
         cz=centroid(ij,3)
         ux=0. ; uy=0. ; uz=0.
         ue_le = ne_rup*ue(ij+u_rup,l)*le(ij+u_rup)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_up,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_up,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_up,3))-cz)
         ue_le = ne_lup*ue(ij+u_lup,l)*le(ij+u_lup)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_up,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_up,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_up,3))-cz)
         ue_le = ne_left*ue(ij+u_left,l)*le(ij+u_left)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_ldown,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_ldown,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_ldown,3))-cz)
         ue_le = ne_ldown*ue(ij+u_ldown,l)*le(ij+u_ldown)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_ldown,1)+xyz_v(ij+z_down,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_ldown,2)+xyz_v(ij+z_down,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_ldown,3)+xyz_v(ij+z_down,3))-cz)
         ue_le = ne_rdown*ue(ij+u_rdown,l)*le(ij+u_rdown)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_rdown,1)+xyz_v(ij+z_down,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_rdown,2)+xyz_v(ij+z_down,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_rdown,3)+xyz_v(ij+z_down,3))-cz)
         ue_le = ne_right*ue(ij+u_right,l)*le(ij+u_right)
         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_rdown,1))-cx)
         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_rdown,2))-cy)
         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_rdown,3))-cz)
         ulon_i = ux*elon_i(ij,1) + uy*elon_i(ij,2) + uz*elon_i(ij,3)
         energy = ulon_i*(1./Ai(ij))
         ulon(ij,l) = ulon(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
         ulon_flux(ij+u_right,l) = ulon_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
         ulon_flux(ij+u_lup,l) = ulon_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
         ulon_flux(ij+u_ldown,l) = ulon_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 ----------------------------------
   !--------------------------------------------------------------------------