     vxl = mxl / dl;
     vyl = myl / dl;
     vzl = mzl / dl;
+    #ifdef DE //PRESSURE_DE
+    E_kin = 0.5 * dl * ( vxl*vxl + vyl*vyl + vzl*vzl );
+    pl = Get_Pressure_From_DE( El, El - E_kin, dgel, gamma ); 
+    #else
     pl  = (El - 0.5*dl*(vxl*vxl + vyl*vyl + vzl*vzl)) * (gamma - 1.0);
+    #endif //PRESSURE_DE
     pl  = fmax(pl, (Real) TINY_NUMBER);
     #ifdef SCALAR
     for (int i=0; i<NSCALARS; i++) {

    vxr = mxr / dr;
    vyr = myr / dr;
    vzr = mzr / dr;
    pr  = (Er - 0.5*dr*(vxr*vxr + vyr*vyr + vzr*vzr)) * (gamma - 1.0);
    pr  = fmax(pr, (Real) TINY_NUMBER);    
    #ifdef SCALAR
    for (int i=0; i<NSCALARS; i++) {

#include"global_cuda.h"
#include"hllc_cuda.h"

#ifdef DE //PRESSURE_DE
#include"hydro_cuda.h"
#endif


/*! \fn Calculate_HLLC_Fluxes_CUDA(Real *dev_bounds_L, Real *dev_bounds_R, Real *dev_flux, int nx, int ny, int nz, int n_ghost, Real gamma, int dir, int n_fields)
