diff --git a/src/calc_ion_temperature.cpp b/src/calc_ion_temperature.cpp
index 115b539b..4b0d02f5 100644
--- a/src/calc_ion_temperature.cpp
+++ b/src/calc_ion_temperature.cpp
@@ -78,7 +78,7 @@ void Ions::calc_ion_temperature(Neutrals neutrals, Grid grid,
 //    // First ion species only, currently is O+
 //    nSpecs = 1;
 //  else
-    nSpecs = nSpecies;
+  nSpecs = nSpecies;
 
   if (report.test_verbose(4)) {
     std::cout << "Bulk ion temp flag: " << input.get_do_calc_bulk_ion_temp() ?
@@ -86,71 +86,71 @@ void Ions::calc_ion_temperature(Neutrals neutrals, Grid grid,
     std::cout << " so 'number of ions' is " << nSpecs << "\n";
   }
 
-  // Loop over all species or assume only bulk calculation
-  for (iIon = 0; iIon < nSpecs; iIon++) {
-    //std::cout << "iIon : " << iIon << "\n";
-    Mi = species[iIon].mass / cAMU;
+  calc_lambda();
+
+  if (input.get_do_calc_bulk_ion_temp()) {
     for (iLon = nGCs; iLon < nLons - nGCs; iLon++) {
       for (iLat = nGCs; iLat < nLats - nGCs; iLat++) {
-        //std::cout << "iLon, iLat : " << iLon << " " << iLat << "\n";
-        // ---------------------------------------------------------------------
-        // Calculate heat flux (conduction) in 1D; loop over all lat,lon
-        // ---------------------------------------------------------------------
-        ratios.zeros();
-        for (jIon = 0; jIon < nSpecies; jIon++) {
-          if (jIon != iIon) {
-            Mj = species[jIon].mass / cAMU;
-            density_ratio = species[jIon].density_scgc.tube(iLon, iLat) / 
-              species[iIon].density_scgc.tube(iLon, iLat);
-            density_ratio.clamp(0.001, 1000.0);
-            ratios = ratios + density_ratio * 
-              (species[jIon].charge * species[jIon].charge / 
-               species[iIon].charge / species[iIon].charge) *
-              sqrt(Mj / (Mi + Mj)) *
-              (3 * Mi * Mi + 1.6 * Mi * Mj + 1.3 * Mj * Mj)/
-              ((Mi + Mj) * (Mi + Mj));
-          }
-        }
-        temp1d = species[iIon].temperature_scgc.tube(iLon, iLat);
-        temp1d(0) = neutrals.temperature_scgc(iLon, iLat, 0);
-        temp1d(1) = neutrals.temperature_scgc(iLon, iLat, 1);
-        lambda1d = 3.1e6 / sqrt(Mi) / pow(species[iIon].charge, 4) *
-           pow(temp1d, 2.5) % (1 + 1.75 * ratios) * cE;
+        temp1d = temperature_scgc.tube(iLon, iLat);
+        lambda1d = lambda.tube(iLon, iLat);
         lambda1d(1) = lambda1d(2);
         lambda1d(0) = lambda1d(2);
-        //lambda1d = 25.0 * cKB * pow(temp1d, 2.5) * (cKB / species[iIon].mass)
-        //           / species[iIon].nu_ion_ion[iIon] / 8.0;
-        front1d  = 3.0 / 2.0 * cKB * species[iIon].density_scgc.tube(iLon, iLat);
+        front1d  = 3.0 / 2.0 * cKB * density_scgc.tube(iLon, iLat);
         dalt1d   = grid.dalt_lower_scgc.tube(iLon, iLat);
         sources1d = (heating_neutral_friction_scgc.tube(iLon, iLat) +
                     heating_neutral_heat_transfer_scgc.tube(iLon, iLat));
         sources1d = sources1d / front1d;
-
-        //std::cout << "iIon lab : " << iIon << "\n" << lambda1d << "\n source:\n " << sources1d << "\n";
-
         conduction1d.zeros();    // reset temp variable to zero
-
         conduction1d = solver_conduction(temp1d, 
                                          lambda1d, 
                                          front1d, 
                                          sources1d, 
                                          dalt1d,
-                                         dt/100., 
+                                         dt/10., 
                                          nGCs, 
                                          false);
-
         // The conduction solver gives Tnew-Told, so divide by dt
         conduction1d.clamp(200, 5000);
-        species[iIon].temperature_scgc.tube(iLon, iLat) = conduction1d;
-
-        //std::cout << "temp : " << conduction1d << "\n";
-
-      } // Lats
-    } // Lons
-  } // Ions
+        temperature_scgc.tube(iLon, iLat) = conduction1d;
+      }
+    }
+    for (iIon = 0; iIon < nSpecies; iIon++)
+      species[iIon].temperature_scgc = temperature_scgc;
 
-  //if (!input.get_do_calc_bulk_ion_temp()) {
-    // Use the density averaged temperature to fill the bulk temperature
+  } else {
+
+    // Loop over all species or assume only bulk calculation
+    for (iIon = 0; iIon < nSpecs; iIon++) {
+      for (iLon = nGCs; iLon < nLons - nGCs; iLon++) {
+        for (iLat = nGCs; iLat < nLats - nGCs; iLat++) {
+          temp1d = species[iIon].temperature_scgc.tube(iLon, iLat);
+          temp1d(0) = neutrals.temperature_scgc(iLon, iLat, 0);
+          temp1d(1) = neutrals.temperature_scgc(iLon, iLat, 1);
+          lambda1d = species[iIon].lambda.tube(iLon, iLat);
+          lambda1d(1) = lambda1d(2);
+          lambda1d(0) = lambda1d(2);
+          front1d  = 3.0 / 2.0 * cKB * species[iIon].density_scgc.tube(iLon, iLat);
+          dalt1d   = grid.dalt_lower_scgc.tube(iLon, iLat);
+          sources1d = (species[iIon].heating_neutral_friction_scgc.tube(iLon, iLat) +
+                      species[iIon].heating_neutral_heat_transfer_scgc.tube(iLon, iLat));
+          sources1d = sources1d / front1d;
+
+          conduction1d.zeros();    // reset temp variable to zero
+          conduction1d = solver_conduction(temp1d, 
+                                           lambda1d, 
+                                           front1d, 
+                                           sources1d, 
+                                           dalt1d,
+                                           dt/10., 
+                                           nGCs, 
+                                           false);
+
+          // The conduction solver gives Tnew-Told, so divide by dt
+          conduction1d.clamp(200, 5000);
+          species[iIon].temperature_scgc.tube(iLon, iLat) = conduction1d;
+        } // Lats
+      } // Lons
+    } // Ions
     tempT.zeros();
     tempD.zeros();
 
@@ -162,19 +162,8 @@ void Ions::calc_ion_temperature(Neutrals neutrals, Grid grid,
     }
 
     temperature_scgc = tempT / tempD;
-  //}
-/*
-  if (input.get_do_calc_bulk_ion_temp()) {
-    // Add temperature terms together to advance bulk ion temperature
-    temperature_scgc = temperature_scgc + dt * (conduction_scgc);
-
-    // Use the bulk ion temperature to fill all ion specie temperatures
-    for (iIon = 0; iIon < nSpecies; iIon++)
-      species[iIon].temperature_scgc = temperature_scgc;
   }
-*/
 
-  //std::cout << "ion temp : " << temperature_scgc(2,2,5) << " " << neutrals.temperature_scgc(2,2,5) << "\n";
   report.exit(function);
   return;
 }