     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
+    CudaCheckError();
     #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

    #endif

    // Apply cooling
    #ifdef TIME
    cudaEventRecord(start, 0);
    #endif //TIME 
    #ifdef COOLING_GPU
    cooling_kernel<<<dim1dGrid,dim1dBlock>>>(dev_conserved, nx_s, ny_s, nz_s, n_ghost, n_fields, dt, gama, dev_dt_array);  
    CudaCheckError();
    #endif
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    cool += elapsedTime;
    #endif //TIME     
 
    // Step 7: Calculate the next time step
    #ifdef TIME
    cudaEventRecord(start, 0);
    #endif //TIME 
    Calc_dt_3D<<<dim1dGrid,dim1dBlock>>>(dev_conserved, nx_s, ny_s, nz_s, n_ghost, dx, dy, dz, dev_dti_array, gama);
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    //printf("dti calc: %5.3f ms\n", elapsedTime);
    //fflush(stdout);
    dti += elapsedTime;
    #endif //TIME  
    CudaCheckError();

    // copy the updated conserved variable array back to the CPU
    #ifdef TIME
    cudaEventRecord(start, 0);
    #endif //TIME    
    CudaSafeCall( cudaMemcpy(tmp2, dev_conserved, n_fields*BLOCK_VOL*sizeof(Real), cudaMemcpyDeviceToHost) );
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    //printf("Copy from GPU: %5.3f ms\n", elapsedTime);
    //fflush(stdout);
    cpfr += elapsedTime;
    #endif //TIME    

    // copy the updated conserved variable array from the buffer into the host_conserved array on the CPU
    #ifdef TIME
    cudaEventRecord(start, 0);
    #endif //TIME 
    host_return_block_3D(nx, ny, nz, nx_s, ny_s, nz_s, n_ghost, block, block1_tot, block2_tot, block3_tot, remainder1, remainder2, remainder3, BLOCK_VOL, host_conserved1, buffer, n_fields);
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    //printf("CPU copying: %5.3f ms\n", elapsedTime);
    //fflush(stdout);
    buff += elapsedTime;
    #endif //TIME    

    // copy the dti array onto the CPU
    #ifdef TIME
    cudaEventRecord(start, 0);
    #endif //TIME 
    CudaSafeCall( cudaMemcpy(host_dti_array, dev_dti_array, ngrid*sizeof(Real), cudaMemcpyDeviceToHost) );
    // find maximum inverse timestep from CFL condition
    for (int i=0; i<ngrid; i++) {
      max_dti = fmax(max_dti, host_dti_array[i]);
    }
    #ifdef TIME
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    //printf("dti copying & calc: %5.3f ms\n", elapsedTime);
    //fflush(stdout);
    dti += elapsedTime;
    #endif //TIME  
    #ifdef COOLING_GPU
    // copy the dt array from cooling onto the CPU
    CudaSafeCall( cudaMemcpy(host_dt_array, dev_dt_array, ngrid*sizeof(Real), cudaMemcpyDeviceToHost) );

  Free_Memory_VL_3D();
  #endif


  // return the maximum inverse timestep
  return max_dti;
