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