     Calculate_Roe_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Lz, Q_Rz, F_z, nx_s, ny_s, nz_s, n_ghost, gama, 2, n_fields);
     #endif //ROE
     #ifdef HLLC 
-    #ifdef TIME
-    cudaEventRecord(start, 0);
-    #endif //TIME 
     Calculate_HLLC_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Lx, Q_Rx, F_x, nx_s, ny_s, nz_s, n_ghost, gama, 0, n_fields);
-    #ifdef TIME
-    cudaEventRecord(stop, 0);
-    cudaEventSynchronize(stop);
-    cudaEventElapsedTime(&elapsedTime, start, stop);
-    //printf("x2 fluxes: %5.3f ms\n", elapsedTime);
-    //fflush(stdout);
-    r2x += elapsedTime;
-    cudaEventRecord(start, 0);
-    #endif //TIME 
     Calculate_HLLC_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Ly, Q_Ry, F_y, nx_s, ny_s, nz_s, n_ghost, gama, 1, n_fields);
-    #ifdef TIME
-    cudaEventRecord(stop, 0);
-    cudaEventSynchronize(stop);
-    cudaEventElapsedTime(&elapsedTime, start, stop);
-    //printf("y2 fluxes: %5.3f ms\n", elapsedTime);
-    //fflush(stdout);
-    r2y += elapsedTime;
-    cudaEventRecord(start, 0);
-    #endif //TIME 
     Calculate_HLLC_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Lz, Q_Rz, F_z, nx_s, ny_s, nz_s, n_ghost, gama, 2, n_fields);
-    #ifdef TIME
-    cudaEventRecord(stop, 0);
-    cudaEventSynchronize(stop);
-    cudaEventElapsedTime(&elapsedTime, start, stop);
-    //printf("z2 fluxes: %5.3f ms\n", elapsedTime);
-    //fflush(stdout);
-    r2z += elapsedTime;
-    cudaEventRecord(start, 0);
-    #endif //TIME 
     #endif //HLLC
     CudaCheckError();
 
 
     // Step 6: Update the conserved variable array
-    #ifdef TIME
-    cudaEventRecord(start, 0);
-    #endif //TIME 
     Update_Conserved_Variables_3D<<<dim1dGrid,dim1dBlock>>>(dev_conserved, F_x, F_y, F_z, nx_s, ny_s, nz_s, x_off_s, y_off_s, z_off_s, n_ghost, dx, dy, dz, xbound, ybound, zbound, dt, gama, n_fields);
-    #ifdef TIME
-    cudaEventRecord(stop, 0);
-    cudaEventSynchronize(stop);
-    cudaEventElapsedTime(&elapsedTime, start, stop);
-    //printf("conserved variable update: %5.3f ms\n", elapsedTime);
-    //fflush(stdout);
-    cvu += elapsedTime;
-    #endif //TIME     
     CudaCheckError();
 
     #ifdef DE

    Calculate_Roe_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Ly, Q_Ry, F_y, nx_s, ny_s, nz_s, n_ghost, gama, 1, n_fields);
    Calculate_Roe_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Lz, Q_Rz, F_z, nx_s, ny_s, nz_s, n_ghost, gama, 2, n_fields);
    #endif //ROE
    #ifdef HLLC 
    #ifdef TIME
    cudaEventRecord(start, 0);
    #endif //TIME 
    Calculate_HLLC_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Lx, Q_Rx, F_x, nx_s, ny_s, nz_s, n_ghost, gama, 0, n_fields);
    CudaCheckError();
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    //printf("x1 fluxes: %5.3f ms\n", elapsedTime);
    //fflush(stdout);
    r1x += elapsedTime;
    cudaEventRecord(start, 0);
    #endif //TIME 
    Calculate_HLLC_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Ly, Q_Ry, F_y, nx_s, ny_s, nz_s, n_ghost, gama, 1, n_fields);
    CudaCheckError();
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    //printf("y1 fluxes: %5.3f ms\n", elapsedTime);
    //fflush(stdout);
    r1y += elapsedTime;
    cudaEventRecord(start, 0);
    #endif //TIME 
    Calculate_HLLC_Fluxes_CUDA<<<dim1dGrid,dim1dBlock>>>(Q_Lz, Q_Rz, F_z, nx_s, ny_s, nz_s, n_ghost, gama, 2, n_fields);
    CudaCheckError();
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    //printf("z1 fluxes: %5.3f ms\n", elapsedTime);
    //fflush(stdout);
    r1z += elapsedTime;
    cudaEventRecord(start, 0);
    #endif //TIME 
    #endif //HLLC
    CudaCheckError();


    // Step 3: Update the conserved variables half a timestep 
    #ifdef TIME
    cudaEventRecord(start, 0);
    #endif //TIME   
    Update_Conserved_Variables_3D_half<<<dim1dGrid,dim1dBlock>>>(dev_conserved, dev_conserved_half, F_x, F_y, F_z, nx_s, ny_s, nz_s, n_ghost, dx, dy, dz, 0.5*dt, gama, n_fields);
    CudaCheckError();
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    //printf("Half step update: %5.3f ms\n", elapsedTime);
    //fflush(stdout);
    ie += elapsedTime;
    #endif //TIME 

    // Step 4: Construct left and right interface values using updated conserved variables
    #ifdef PCM

    PPMP_cuda<<<dim1dGrid,dim1dBlock>>>(dev_conserved_half, Q_Ly, Q_Ry, nx_s, ny_s, nz_s, n_ghost, dy, dt, gama, 1, n_fields);
    PPMP_cuda<<<dim1dGrid,dim1dBlock>>>(dev_conserved_half, Q_Lz, Q_Rz, nx_s, ny_s, nz_s, n_ghost, dz, dt, gama, 2, n_fields);
    #endif //PPMP
    CudaCheckError();
    #ifdef PPMC
    PPMC_cuda<<<dim1dGrid,dim1dBlock>>>(dev_conserved_half, Q_Lx, Q_Rx, nx_s, ny_s, nz_s, n_ghost, dx, dt, gama, 0, n_fields);
    CudaCheckError();
    PPMC_cuda<<<dim1dGrid,dim1dBlock>>>(dev_conserved_half, Q_Ly, Q_Ry, nx_s, ny_s, nz_s, n_ghost, dy, dt, gama, 1, n_fields);
    CudaCheckError();
    PPMC_cuda<<<dim1dGrid,dim1dBlock>>>(dev_conserved_half, Q_Lz, Q_Rz, nx_s, ny_s, nz_s, n_ghost, dz, dt, gama, 2, n_fields);
    CudaCheckError();
    #endif //PPMC
    

    // Step 5: Calculate the fluxes again
