diff --git a/config/output_spec.csv b/config/output_spec.csv
index a8478b3c..4a7c862d 100644
--- a/config/output_spec.csv
+++ b/config/output_spec.csv
@@ -85,6 +85,7 @@ SNOWLAYERDZ,Snow pack thickness by layer,m,,,invalid,invalid,invalid,forced,doub
 SNOWLAYERTEMP,Snow temperature by layer,degree_C,,,invalid,invalid,invalid,forced,double,
 SNOWSTART,DOY of first snow fall,DOY,,invalid,invalid,invalid,invalid,invalid,int,
 SNOWTHICK,Snow pack thickness,m,,,,invalid,invalid,invalid,double,
+SOC,Soil organic C total,g/m2,,,invalid,invalid,invalid,,double,
 SOMA,Soil organic C active,g/m2,,,invalid,invalid,invalid,,double,
 SOMCR,Soil organic C chemically resistant,g/m2,,,invalid,invalid,invalid,,double,
 SOMPR,Soil organic C physically resistant,g/m2,,,invalid,invalid,invalid,,double,
diff --git a/include/OutputHolder.h b/include/OutputHolder.h
index 46abe38b..823ec9ab 100644
--- a/include/OutputHolder.h
+++ b/include/OutputHolder.h
@@ -68,6 +68,9 @@ class OutputHolder{
   std::vector<double> somrawc_tot_for_output;
   std::vector<std::array<double, MAX_SOI_LAY>> somrawc_for_output;
 
+  std::vector<double> soc_tot_for_output;
+  std::vector<std::array<double, MAX_SOI_LAY>> soc_for_output;
+
   std::vector<double> avln_tot_for_output;
   std::vector<std::array<double, MAX_SOI_LAY>> avln_for_output;
 
diff --git a/src/Runner.cpp b/src/Runner.cpp
index a33cabf1..f066fbe3 100644
--- a/src/Runner.cpp
+++ b/src/Runner.cpp
@@ -4755,6 +4755,93 @@ void Runner::output_netCDF(std::map<std::string, OutputSpec> &netcdf_outputs, in
   map_itr = netcdf_outputs.end();
 
 
+  // SOC
+  map_itr = netcdf_outputs.find("SOC");
+  if (map_itr != netcdf_outputs.end()) {
+    BOOST_LOG_SEV(glg, debug) << "NetCDF output: SOC";
+    curr_spec = map_itr->second;
+
+    #pragma omp critical(outputSOC)
+    {
+      // By layer
+      if (curr_spec.layer) {
+        // Monthly
+        if (curr_spec.monthly) {
+          std::array<double, MAX_SOI_LAY> m_soc;
+          for (int il=0; il<MAX_SOI_LAY; il++) {
+            m_soc[il] = cohort.bdall->m_sois.rawc[il]
+                      + cohort.bdall->m_sois.soma[il]
+                      + cohort.bdall->m_sois.sompr[il]
+                      + cohort.bdall->m_sois.somcr[il];
+          }
+          outhold.soc_for_output.push_back(m_soc);
+
+          if (output_this_timestep) {
+            output_nc_4dim(&curr_spec, file_stage_suffix, &outhold.soc_for_output[0][0], MAX_SOI_LAY, month_start_idx, months_to_output);
+            outhold.soc_for_output.clear();
+          }
+        }
+        // Yearly
+        else if (curr_spec.yearly) {
+          std::array<double, MAX_SOI_LAY> y_soc;
+          for (int il=0; il<MAX_SOI_LAY; il++) {
+            y_soc[il] = cohort.bdall->y_sois.rawc[il]
+                      + cohort.bdall->y_sois.soma[il]
+                      + cohort.bdall->y_sois.sompr[il]
+                      + cohort.bdall->y_sois.somcr[il];
+          }
+          outhold.soc_for_output.push_back(y_soc);
+
+          if (output_this_timestep) {
+            output_nc_4dim(&curr_spec, file_stage_suffix, &outhold.soc_for_output[0][0], MAX_SOI_LAY, year_start_idx, years_to_output);
+            outhold.soc_for_output.clear();
+          }
+        }
+      }
+      // Total, instead of by layer
+      else if (!curr_spec.layer) {
+        // Monthly
+        if (curr_spec.monthly) {
+
+          double m_soc_tot = 0.0;
+          for (int il=0; il<MAX_SOI_LAY; il++) {
+            m_soc_tot += cohort.bdall->m_sois.rawc[il]
+                       + cohort.bdall->m_sois.soma[il]
+                       + cohort.bdall->m_sois.sompr[il]
+                       + cohort.bdall->m_sois.somcr[il];
+          }
+
+          outhold.soc_tot_for_output.push_back(m_soc_tot);
+
+          if (output_this_timestep) {
+            output_nc_3dim(&curr_spec, file_stage_suffix, &outhold.soc_tot_for_output[0], 1, month_start_idx, months_to_output);
+            outhold.soc_tot_for_output.clear();
+          }
+        }
+        // Yearly
+        else if (curr_spec.yearly) {
+
+          double y_soc_tot = 0.0;
+          for (int il=0; il<MAX_SOI_LAY; il++) {
+            y_soc_tot += cohort.bdall->y_sois.rawc[il]
+                       + cohort.bdall->y_sois.soma[il]
+                       + cohort.bdall->y_sois.sompr[il]
+                       + cohort.bdall->y_sois.somcr[il];
+          }
+
+          outhold.soc_tot_for_output.push_back(y_soc_tot);
+
+          if (output_this_timestep) {
+            output_nc_3dim(&curr_spec, file_stage_suffix, &outhold.soc_tot_for_output[0], 1, year_start_idx, years_to_output);
+            outhold.soc_tot_for_output.clear();
+          }
+        }
+      }
+    } //end critical(outputSOC)
+  } //end SOC
+  map_itr = netcdf_outputs.end();
+
+
   //SOMA - soil organic matter, active
   map_itr = netcdf_outputs.find("SOMA");
   if(map_itr != netcdf_outputs.end()){