diff --git a/2_SWSF_p4of5_Code_v51.R b/2_SWSF_p4of5_Code_v51.R
index f288ff6b..d550e49d 100644
--- a/2_SWSF_p4of5_Code_v51.R
+++ b/2_SWSF_p4of5_Code_v51.R
@@ -2146,7 +2146,7 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer
 					if(print.debug) print("Start of get SiteClimate")
 					do.C4vars <- any(create_treatments == "PotentialNaturalVegetation_CompositionShrubsC3C4_Paruelo1996") || aon$dailyC4_TempVar
 					#redo SiteClimate_Ambient
-					SiteClimate_Ambient <- sw_SiteClimate_Ambient(weatherList=i_sw_weatherList[[1]], year.start=min(simTime$useyrs), year.end=max(simTime$useyrs), do.C4vars=do.C4vars, simTime2=simTime2)
+					SiteClimate_Ambient <- sw_SiteClimate_Ambient(weatherList=i_sw_weatherList[[1]], year.start=min(simTime$useyrs), year.end=max(simTime$useyrs), do.C4vars=do.C4vars, simTime = simTime, simTime2=simTime2)
 				}
 			}
 			if(!getScenarioWeatherDataFromDatabase) {
@@ -2248,11 +2248,11 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer
 					if(do.C4vars){
 						SiteClimate_Scenario$dailyTempMin <- SiteClimate_Ambient$dailyTempMin + t_min_f[simTime2$month_ForEachUsedDay]
 						SiteClimate_Scenario$dailyTempMean <- SiteClimate_Ambient$dailyTempMean + tmean_f[simTime2$month_ForEachUsedDay]
-						SiteClimate_Scenario$dailyC4vars <- sw_dailyC4_TempVar(SiteClimate_Scenario$dailyTempMin, SiteClimate_Scenario$dailyTempMean, simTime2)
+						SiteClimate_Scenario$dailyC4vars <- sw_dailyC4_TempVar(SiteClimate_Scenario$dailyTempMin, SiteClimate_Scenario$dailyTempMean, simTime, simTime2)
 					}
 				}
 			} else {
-				SiteClimate_Scenario <- sw_SiteClimate_Ambient(weatherList=i_sw_weatherList[[sc]], year.start=min(simTime$useyrs), year.end=max(simTime$useyrs), do.C4vars=do.C4vars, simTime2=simTime2)
+				SiteClimate_Scenario <- sw_SiteClimate_Ambient(weatherList=i_sw_weatherList[[sc]], year.start=min(simTime$useyrs), year.end=max(simTime$useyrs), do.C4vars=do.C4vars, simTime = simTime, simTime2=simTime2)
 				if(sc > 1){
 					ppt_sc <- (temp <- swWeather_MonScalingParams(swRunScenariosData[[sc]]))[,1]
 					t_max <- temp[,2]
@@ -3388,14 +3388,15 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer
 				}
 #TODO(drs): progress state
 			#23
-				if(any(simulation_timescales=="daily") & aon$dailyC4_TempVar){#Variables to estimate percent C4 species in North America: Teeri JA, Stowe LG (1976) Climatic patterns and the distribution of C4 grasses in North America. Oecologia, 23, 1-12.
+				if(any(simulation_timescales=="daily") & aon$dailyC4_TempVar){
+					# Variables to estimate percent C4 species in North America: Teeri JA, Stowe LG (1976) Climatic patterns and the distribution of C4 grasses in North America. Oecologia, 23, 1-12.
 					if(print.debug) print("Aggregation of dailyC4_TempVar")
 					if(!exists("temp.dy")) temp.dy <- get_Temp_dy(sc, runData, simTime)
 
           temp <- sw_dailyC4_TempVar(dailyTempMin = temp.dy$min,
                                     dailyTempMean = temp.dy$mean,
-                                    simTime2) #accountNSHemispheres_agg
-          temp <- as.numeric(temp)
+                                    simTime, simTime2,
+                                    return_yearly = TRUE)
 
 					resMeans[nv:(nv+2)] <- temp[1:3]
 					resSDs[nv:(nv+2)] <- temp[4:6]
diff --git a/2_SWSF_p5of5_Functions_v51.R b/2_SWSF_p5of5_Functions_v51.R
index bd6efc38..ba7ef272 100644
--- a/2_SWSF_p5of5_Functions_v51.R
+++ b/2_SWSF_p5of5_Functions_v51.R
@@ -405,7 +405,7 @@ adjustLayersDepth <- compiler::cmpfun(function(layers_depth, d) round(layers_dep
 getLayersWidth <- compiler::cmpfun(function(layers_depth) diff(c(0, layers_depth)))
 setLayerSequence <- compiler::cmpfun(function(d) seq_len(d))
 
-sw_dailyC4_TempVar <- compiler::cmpfun(function(dailyTempMin, dailyTempMean, simTime2) {
+sw_dailyC4_TempVar <- compiler::cmpfun(function(dailyTempMin, dailyTempMean, simTime, simTime2, return_yearly = FALSE) {
   #Variables to estimate percent C4 species in North America: Teeri JA, Stowe LG (1976) Climatic patterns and the distribution of C4 grasses in North America. Oecologia, 23, 1-12.
 
   temp7 <- simTime2$month_ForEachUsedDay_NSadj == 7
@@ -420,19 +420,37 @@ sw_dailyC4_TempVar <- compiler::cmpfun(function(dailyTempMin, dailyTempMean, sim
   DegreeDaysAbove65F_DaysC <- tapply(temp_base65F, simTime2$year_ForEachUsedDay_NSadj, sum)
 
   nyrs <- seq_along(Month7th_MinTemp_C) #if southern Hemisphere, then 7th month of last year is not included
-  temp <- cbind(Month7th_MinTemp_C[nyrs],
+  ydat <- cbind(Month7th_MinTemp_C[nyrs],
                 LengthFreezeFreeGrowingPeriod_Days[nyrs],
                 DegreeDaysAbove65F_DaysC[nyrs])
-  res <- c(apply(temp, 2, mean), apply(temp, 2, sd))
-  temp <- c("Month7th_NSadj_MinTemp_C",
-            "LengthFreezeFreeGrowingPeriod_NSadj_Days",
-            "DegreeDaysAbove65F_NSadj_DaysC")
-  names(res) <- c(temp, paste0(temp, ".sd"))
+  dimnames(ydat) <- list(NULL,
+                        c("Month7th_NSadj_MinTemp_C",
+                          "LengthFreezeFreeGrowingPeriod_NSadj_Days",
+                          "DegreeDaysAbove65F_NSadj_DaysC")
+                        )
+
+  if (return_yearly) {
+    yrs_have <- as.integer(names(Month7th_MinTemp_C[nyrs]))
+
+    if (all(yrs_have == simTime$useyrs)) {
+      res <- cbind(Years = yrs_have, ydat)
+
+    } else {
+      res <- matrix(NA, nrow = simTime$no.useyr, ncol = 4,
+                    dimnames = list(NULL, c("Years", colnames(ydat))))
+      res[, "Years"] <- simTime$useyrs
+      res[, -1] <- ydat[match(simTime$useyrs, yrs_have), ]
+    }
+
+  } else {
+    res <- c(apply(ydat, 2, mean), apply(ydat, 2, sd))
+    names(res) <- c(colnames(ydat), paste0(colnames(ydat), ".sd"))
+  }
 
   res
 })
 
-sw_SiteClimate_Ambient <- compiler::cmpfun(function(weatherList, year.start, year.end, do.C4vars = FALSE, simTime2 = NULL) {
+sw_SiteClimate_Ambient <- compiler::cmpfun(function(weatherList, year.start, year.end, do.C4vars = FALSE, simTime = NULL, simTime2 = NULL) {
   sw.weather.suffix <- as.numeric(names(weatherList))
   itemp <- year.start <= sw.weather.suffix & year.end >= sw.weather.suffix
   years <- sw.weather.suffix[itemp]
@@ -479,7 +497,7 @@ sw_SiteClimate_Ambient <- compiler::cmpfun(function(weatherList, year.start, yea
   if (do.C4vars) {
     res[["dailyTempMin"]] <- dailyTempMin
     res[["dailyTempMean"]] <- dailyTempMean
-    res[["dailyC4vars"]] <- sw_dailyC4_TempVar(dailyTempMin, dailyTempMean, simTime2)
+    res[["dailyC4vars"]] <- sw_dailyC4_TempVar(dailyTempMin, dailyTempMean, simTime, simTime2)
 
   } else {
     res[["dailyTempMin"]] <- res[["dailyTempMean"]] <- res[["dailyC4vars"]] <- NA