             / units;    
 }
 
-//Calculation of comp.
-// Compton cooling
-__device__ double comp_rate(double units)
-{
-    return 5.65e-36 / units;
-}
+// //Calculation of comp.
+// // Compton cooling
+// __device__ double comp_rate(double units)
+// {
+//    / units
+//   // Set compton cooling coefficients (and temperature)
+//   comp1 = 5.65e-36 * (1._DKIND + zr)**4
+//   comp2 = 2.73_DKIND * (1._DKIND + zr)
+//   !                  Compton cooling or heating
+// 
+//        &           - comp1      * (tgas(i) - comp2)     * myde(i)*dom_inv
+//     return 5.65e-36 / units;
+// }
+// 
+// __device__ Real Cooling_Rate_Compton_CMB_MillesOstriker01( Real n_e, Real temp, Real z ){
+//   // M. Coleman Miller and Eve C. Ostriker 2001 (https://iopscience.iop.org/article/10.1086/323321/fulltext/)
+//   Real T_3, T_cm_3;
+//   T_3 = temp / 1e3;
+//   // T_cm_3 = 0; //# Don't know this value
+//   T_cm_3 = 2.726 * ( 1 + z ) / 1e3; //# Don't know this value
+//   return 5.6e-33 * n_e * pow( 1 + z , 4) * ( T_3 - T_cm_3 ) ;
+// }
 
 
 

  if ( dt < 0 ) printf( "Chem_GPU: Negative dt: %e  dt_U: %e  dt_HI: %e  dt_e %e   n_HI: %e\n", dt, dt_min_U, dt_min_HI, dt_min_e, TS.n_HI );
  // if ( print ) printf("Chem GPU: dt_U: %e  dt_HI: %e  dt_e: %e  \n", dt_min_U, dt_min_HI, dt_min_e );
  
  // Record initial states of HI and e to compute derivatives
  init_n_HI = TS.n_HI;
  init_n_e  = TS.n_e;

  TS.U += dQ_dt * dt / rho ;
  
  // Bound lower
  TS.n_HI    = fmax( upd_n_HI,    n_min);
  TS.n_HII   = fmax( upd_n_HII,   n_min);
  TS.n_HeI   = fmax( upd_n_HeI,   n_min);
  TS.n_HeII  = fmax( upd_n_HeII,  n_min);
  TS.n_HeIII = fmax( upd_n_HeIII, 1e-5*n_min);
  
  // Use charge conservation to determine electron fraction
  TS.n_e = TS.n_HII + TS.n_HeII + 2 * TS.n_HeIII;
