     vxl = dev_bounds_L[o1*n_cells + tid]/dl;
     vyl = dev_bounds_L[o2*n_cells + tid]/dl;
     vzl = dev_bounds_L[o3*n_cells + tid]/dl;
+    #ifdef DE //PRESSURE_DE
+    E = dev_bounds_L[4*n_cells + tid];
+    E_kin = 0.5 * dl * ( vxl*vxl + vyl*vyl + vzl*vzl );
+    dge = dev_bounds_L[(n_fields-1)*n_cells + tid];
+    pl = Get_Pressure_From_DE( E, E - E_kin, dge, gamma ); 
+    #else   
     pl  = (dev_bounds_L[4*n_cells + tid] - 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++) {

#include"global_cuda.h"
#include"exact_cuda.h"


/*! \fn Calculate_Exact_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)
 *  \brief Exact Riemann solver based on the Fortran code given in Sec. 4.9 of Toro (1999). */

    }
    #endif
    #ifdef DE
    ger = dge / dr;
    #endif


