diff --git a/mtcars_nestedplotting.R b/mtcars_nestedplotting.R
new file mode 100644
index 0000000..905a969
--- /dev/null
+++ b/mtcars_nestedplotting.R
@@ -0,0 +1,53 @@
+library(dplyr)
+library(tidyr)
+library(purrr)
+library(ggplot2)
+
+
+by_cyl <- mtcars %>%  
+  group_by(cyl, gear) %>%  
+  nest() %>%  
+  rename(cyl_data = data)
+
+by_gear <- mtcars %>% 
+  group_by(gear) %>%  
+  nest() %>%  
+  rename(gear_data = data)
+
+mtcars_nest <- left_join(by_cyl, by_gear, by = 'gear')
+
+mtcars_nest <- mtcars_nest %>% 
+  mutate(
+    map(cyl_data, ~ ggplot(., aes(x = wt, y = mpg)) + 
+                             geom_point() +
+                             geom_smooth(se = TRUE, color = 'blue')
+    )
+  ) %>% 
+  rename(plot_cyl = `map(...)`)
+
+mtcars_nest <- mtcars_nest %>% 
+  mutate(
+    map(gear_data, ~ ggplot(., aes(x = wt, y = mpg)) + 
+          geom_point() +
+          geom_smooth(se = TRUE, color = 'red')
+    )
+  ) %>% 
+  rename(plot_gear = `map(...)`)
+
+mtcars_nest$plot_cyl[1]
+mtcars_nest$plot_gear[1]
+
+#linear_model <- lm(mpg ~ hp + wt, data = mtcars)
+
+#How to get these two plots on one figure?
+mtcars_nest <- mtcars_nest %>% 
+  mutate( 
+    plot_cyl_gear = map2(cyl_data, gear_data,
+                         ~ ggplot(.x, aes(x = wt, y = mpg)) +
+                           geom_smooth(se = TRUE, color = 'blue') +
+                           geom_smooth(data = .y, se = TRUE, color = 'red')
+    ) )
+
+mtcars_nest$plot_cyl_gear[1]
+  
+  
diff --git a/ncplot.r b/ncplot.r
new file mode 100644
index 0000000..10e5a88
--- /dev/null
+++ b/ncplot.r
@@ -0,0 +1,151 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+#workdir <- '~/tmp/us-umb/stand/lnd/hist'
+#setwd(workdir)
+#ncfname <- 'US-UMB_I1850CLM45CN_ad_spinup.clm2.h0.0701-01-01-00000.nc'
+#ncf <- open.nc(ncfname)
+#vars <- as.list(read.nc(nc))
+
+##### build set of datatables, once for each variable of interest
+
+##### create datatable from standard maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb-n/a_t_n-/lnd/hist', pattern = '*.h0.*.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i])
+  a <- (var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365)
+  b <- (var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365)
+  c <- (var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365)
+  d <- (var.get.nc(nc = ncf, "NPP") * 3600 * 24 * 365)
+  t <- (var.get.nc(nc = ncf, "time") / 365)
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #print(m) # for debugging
+}
+
+t <- seq(0, 820, by = 20) # create time variable, step by 20
+stand <- as.data.table(m) # convert matrix to data table
+stand[, "V5"] <- t # overwrite "time" with t, and rename columns
+stand <- rename(stand, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="YEARS"))
+stand <- stand[, RS := "stand"] # add a column for respiration algorithm
+
+##### create datatable from no acclimation maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb-n/mh_t_n-/lnd/hist', pattern = '*.h0.*.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i])
+  a <- (var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365)
+  b <- (var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365)
+  c <- (var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365)
+  d <- (var.get.nc(nc = ncf, "NPP") * 3600 * 24 * 365)
+  t <- (var.get.nc(nc = ncf, "time") / 365)
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #print(m) # for debugging
+}
+
+t <- seq(0, 820, by = 20) # create time variable, step by 20
+noacclim <- as.data.table(m) # convert matrix to data table
+noacclim[, "V5"] <- t # overwrite "time" with t, and rename columns
+noacclim <- rename(noacclim, c("V1" = "LEAF_MR", "V2"="MR","V3" = "GPP", "V4" = "NPP", "V5" = "YEARS"))
+noacclim <- noacclim[, RS := "noacclim"] # add a column for respiration algorithm
+
+##### create datatable from acclimation maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb-n/ch_t_n-/lnd/hist', pattern = '*.h0.*.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i])
+  a <- (var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365)
+  b <- (var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365)
+  c <- (var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365)
+  d <- (var.get.nc(nc = ncf, "NPP") * 3600 * 24 * 365)
+  t <- (var.get.nc(nc = ncf, "time") / 365)
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #print(m) # for debugging
+}
+
+t <- seq(0, 820, by = 20) # create time variable, step by 20
+acclim <- as.data.table(m) # convert matrix to data table
+acclim[, "V5"] <- t # overwrite "time" with t, and rename columns
+acclim <- rename(acclim, c("V1" = "LEAF_MR", "V2" = "MR","V3" = "GPP", "V4" = "NPP", "V5" = "YEARS"))
+acclim <- acclim[, RS := "acclim"] # add a column for respiration algorithm
+
+### concatinate the data tables into single data source
+l <- list(stand, noacclim, acclim)
+dt <- rbindlist(l, use.names=TRUE)
+
+### MAKE FANCY FIGURES
+
+lmr <- ggplot() +
+  geom_line(data = stand, aes(x = YEARS, y = LEAF_MR), color = 'black') +
+  geom_line(data = noacclim, aes(x = YEARS, y = LEAF_MR), color = 'blue') + 
+  geom_line(data = acclim, aes(x = YEARS, y = LEAF_MR), color = 'red') + 
+  scale_x_continuous(name = '') +                                                     # hide x-axis name
+  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"), limits = c(0, 400)) +   # set y-axis to be 0:400
+  theme(axis.text.x = element_blank(), axis.ticks = element_blank())                # hide x-axis tick marks
+  #panel.background  = element_rect(fill = "white", colour = NA) 
+
+mr <- ggplot() + 
+  geom_line(data = stand, aes(x = YEARS, y = MR), color='black') +
+  geom_line(data = noacclim, aes(x = YEARS, y = MR), color='blue') + 
+  geom_line(data = acclim, aes(x = YEARS, y = MR), color='red') + 
+  scale_x_continuous(name = '') +
+  scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")"), limits = c(0, 400)) + 
+  theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+  #panel.background  = element_rect(fill = "white", colour = NA) 
+
+gpp <- ggplot() + 
+  geom_line(data = stand, aes(x = YEARS, y = GPP), color = 'black') +
+  geom_line(data = noacclim, aes(x = YEARS, y = GPP), color = 'blue') + 
+  geom_line(data = acclim, aes(x = YEARS, y = GPP), color = 'red') + 
+  scale_x_continuous(name = 'Years', breaks = c(100, 200, 300, 400, 500, 600, 700, 800), 
+                     labels = c("100", "200", "300", "400", "500", "600", "700", "800")) +  # set x-axis labels to every 100 years
+  scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+  #panel.background  = element_rect(fill = "white", colour = NA)
+  
+npp <- ggplot() + 
+  geom_line(data = stand, aes(x = YEARS, y = NPP), color = 'black') +
+  geom_line(data = noacclim, aes(x = YEARS, y = NPP), color = 'blue') + 
+  geom_line(data = acclim, aes(x = YEARS, y = NPP), color = 'red') + 
+  scale_x_continuous(name = 'Years', breaks = c(100, 200, 300, 400, 500, 600, 700, 800), 
+                     labels = c("100", "200", "300", "400", "500", "600", "700", "800")) +  # set x-axis labels to every 100 years
+  scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+  #panel.background = element_rect(fill = "white", colour = NA)
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(lmr)
+
+
+grid.arrange(lmr, mr, gpp, npp, ncol = 2, nrow = 2, # display all 4 plots as 2 x 2 grid
+             top="Three Respiration Equations") 
+             #bottom = textGrob("original=black, acclim=blue, acclim at 25°C=red", 
+                               #gp = gpar(fontface=3, fontsize=9), 
+                               #hjust=1, x=1))
+
diff --git a/ncplot3.r b/ncplot3.r
new file mode 100644
index 0000000..598df5f
--- /dev/null
+++ b/ncplot3.r
@@ -0,0 +1,170 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/us-umb/acme_t_n-/lnd/hist'
+setwd(workdir)
+ncfname <- 'acme_t_n-_US-UMB_I1850CLM45CN_ad_spinup.clm2.h0.0601-01-01-00000.nc'
+ncf <- open.nc(ncfname)
+vars <- as.list(read.nc(ncf))
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from standard ACME maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/acme_t_n-/lnd/hist', pattern = '*.h0.*.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i])
+  a <- (var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365)
+  b <- (var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365)
+  c <- (var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365)
+  d <- (var.get.nc(nc = ncf, "NPP") * 3600 * 24 * 365)
+  t <- (var.get.nc(nc = ncf, "time") / 365)
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #print(m) # for debugging
+  close.nc(ncf)  
+}
+
+t <- seq(0, 800, by = 20) # create time variable, step by 20
+acme <- as.data.table(m) # convert matrix to data table
+acme[, "V5"] <- t # overwrite "time" with t, and rename columns
+acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="YEARS"))
+acme <- acme[, rs_algorithm := "ACME"] # add a column for respiration algorithm
+
+##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/fb_q_t_n-/lnd/hist', pattern = '*.h0.*.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i])
+  a <- (var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365)
+  b <- (var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365)
+  c <- (var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365)
+  d <- (var.get.nc(nc = ncf, "NPP") * 3600 * 24 * 365)
+  t <- (var.get.nc(nc = ncf, "time") / 365)
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #print(m) # for debugging
+  close.nc(ncf)  
+}
+
+t <- seq(0, 800, by = 20) # create time variable, step by 20
+fb_q <- as.data.table(m) # convert matrix to data table
+fb_q[, "V5"] <- t # overwrite "time" with t, and rename columns
+fb_q <- rename(fb_q, c("V1" = "LEAF_MR", "V2"="MR","V3" = "GPP", "V4" = "NPP", "V5" = "YEARS"))
+fb_q <- fb_q[, rs_algorithm := "FB_Q"] # add a column for respiration algorithm
+
+##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/vb_q_t_n-/lnd/hist', pattern = '*.h0.*.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i])
+  a <- (var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365)
+  b <- (var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365)
+  c <- (var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365)
+  d <- (var.get.nc(nc = ncf, "NPP") * 3600 * 24 * 365)
+  t <- (var.get.nc(nc = ncf, "time") / 365)
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #print(m) # for debugging
+  close.nc(ncf)  
+}
+
+t <- seq(0, 800, by = 20) # create time variable, step by 20
+vb_q <- as.data.table(m) # convert matrix to data table
+vb_q[, "V5"] <- t # overwrite "time" with t, and rename columns
+vb_q <- rename(vb_q, c("V1" = "LEAF_MR", "V2"="MR","V3" = "GPP", "V4" = "NPP", "V5" = "YEARS"))
+vb_q <- vb_q[, rs_algorithm := "VB_Q"] # add a column for respiration algorithm
+
+##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/vb_e_t_n-/lnd/hist', pattern = '*.h0.*.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i])
+  a <- (var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365)
+  b <- (var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365)
+  c <- (var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365)
+  d <- (var.get.nc(nc = ncf, "NPP") * 3600 * 24 * 365)
+  t <- (var.get.nc(nc = ncf, "time") / 365)
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #print(m) # for debugging
+  close.nc(ncf)  
+}
+
+t <- seq(0, 800, by = 20) # create time variable, step by 20
+vb_e <- as.data.table(m) # convert matrix to data table
+vb_e[, "V5"] <- t # overwrite "time" with t, and rename columns
+vb_e <- rename(vb_e, c("V1" = "LEAF_MR", "V2" = "MR","V3" = "GPP", "V4" = "NPP", "V5" = "YEARS"))
+vb_e <- vb_e[, rs_algorithm := "VB_E"] # add a column for respiration algorithm
+
+### concatinate the data tables into single data source
+l <- list(acme, fb_q, vb_q, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+### MAKE FANCY FIGURES
+
+p <- ggplot(dt, aes(x = YEARS, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "red", "green", "blue"))
+
+lmr <- ggplot(dt, aes(x = YEARS, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr +  scale_x_continuous(name = '')                                                     # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"), limits = c(0, 400))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "red", "green", "blue"))
+                               
+mr <- ggplot(dt, aes(x = YEARS, y = MR, color = rs_algorithm)) + geom_line() 
+mr <- mr + scale_x_continuous(name = '')
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")"), limits = c(0, 400)) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+gpp <- ggplot(dt, aes(x = YEARS, y = GPP, color = rs_algorithm)) + geom_line()  
+gpp <- gpp + scale_x_continuous(name = 'Years', breaks = c(100, 200, 300, 400, 500, 600, 700, 800), 
+                     labels = c("100", "200", "300", "400", "500", "600", "700", "800"))        # set x-axis labels to every 100 years
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+npp <- ggplot(dt, aes(x = YEARS, y = NPP, color = rs_algorithm)) + geom_line()  
+npp <- npp + scale_x_continuous(name = 'Years', breaks = c(100, 200, 300, 400, 500, 600, 700, 800), 
+                      labels = c("100", "200", "300", "400", "500", "600", "700", "800"))       # set x-axis labels to every 100 years
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), 
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Four Respiration Equations")
diff --git a/ncplot_annual.R b/ncplot_annual.R
new file mode 100644
index 0000000..7227102
--- /dev/null
+++ b/ncplot_annual.R
@@ -0,0 +1,187 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/output/us-umb/acme_t_n-/lnd/hist'
+setwd(workdir)
+ncfname <- 'acme_t_n-_US-UMB_I20TRCLM45CN.clm2.h0.1982-01-01-00000.nc'
+ncf <- open.nc(ncfname)
+vars <- as.list(read.nc(ncf))
+ncf_data <- read.nc(ncfile = ncf)
+#ncf_data$LEAF_MR
+#ncf_data$XSMRPOOL
+#ncf_data$GPP
+#ncf_data$NPP
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from standard ACME maintenence respiration files
+# open all files and initialize matrix
+
+files <- list.files(path = '~/tmp/output/us-umb/acme_t_n-/lnd/hist', pattern = '*h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1850, 2007, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+acme <- as.data.table(ar2dt) # convert array to data table
+acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+acme <- acme[, rs_algorithm := "acme"] # add a column for respiration algorithm
+
+##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/output/us-umb/fb_q_t_n-/lnd/hist', pattern = '*h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1850, 2007, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+fb_q <- as.data.table(ar2dt) # convert array to data table
+fb_q <- rename(fb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+fb_q <- fb_q[, rs_algorithm := "fb_q"] # add a column for respiration algorithm
+
+##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/output/us-umb/vb_q_t_n-/lnd/hist', pattern = '*h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1850, 2007, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_q <- as.data.table(ar2dt) # convert array to data table
+vb_q <- rename(vb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_q <- vb_q[, rs_algorithm := "vb_q"] # add a column for respiration algorithm
+
+##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/output/us-umb/vb_e_t_n-/lnd/hist', pattern = '*h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1850, 2007, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_e <- as.data.table(ar2dt) # convert array to data table
+vb_e <- rename(vb_e, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+
+### MAKE FANCY FIGURES
+
+### concatinate the data tables into single data source then draw figures
+l <- list(acme, fb_q, vb_q, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+p <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "red", "green", "blue"))
+
+lmr <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr +  scale_x_continuous(name = '', limits = c(1990, 2005))                                                     # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "red", "green", "blue"))
+                               
+mr <- ggplot(dt, aes(x = YEAR, y = MR, color = rs_algorithm)) + geom_line() 
+mr <- mr + scale_x_continuous(name = '', limits = c(1990, 2005))
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")")) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+gpp <- ggplot(dt, aes(x = YEAR, y = GPP, color = rs_algorithm)) + geom_line()  
+gpp <- gpp + scale_x_continuous(name = 'Years', limits = c(1990, 2005))        # set x-axis labels to every 100 years
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+npp <- ggplot(dt, aes(x = YEAR, y = NPP, color = rs_algorithm)) + geom_line()  
+npp <- npp + scale_x_continuous(name = 'Years', limits = c(1990, 2005))       # set x-axis labels to every 100 years
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), 
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Four Respiration Equations")
+
diff --git a/ncplot_br_sa1_annual.R b/ncplot_br_sa1_annual.R
new file mode 100644
index 0000000..50006a5
--- /dev/null
+++ b/ncplot_br_sa1_annual.R
@@ -0,0 +1,324 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+#####SPINUP##############################################################################
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/run/acme_t_n-_BR-Sa1_I1850CLM45CN_ad_spinup/run'
+setwd(workdir)
+ncfname <- 'acme_t_n-_BR-Sa1_I1850CLM45CN_ad_spinup.clm2.h0.0261-01-01-00000.nc'
+ncf <- open.nc(ncfname)
+vars <- as.list(read.nc(ncf))
+ncf_data <- read.nc(ncfile = ncf)
+#ncf_data$LEAF_MR
+#ncf_data$XSMRPOOL
+#ncf_data$GPP
+#ncf_data$NPP
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from standard ACME maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/acme_t_n-_BR-Sa1_I1850CLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365     #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365          #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365         #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365          #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 310, by = 10) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+acme <- as.data.table(ar2dt) # convert array to data table
+acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+acme <- acme[, rs_algorithm := "acme"] # add a column for respiration algorithm
+
+##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# open all files and initialize matrix
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/fb_q_t_n-_BR-Sa1_I1850CLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365    #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365         #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365        #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365         #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 310, by = 10) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+fb_q <- as.data.table(ar2dt) # convert array to data table
+fb_q <- rename(fb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+fb_q <- fb_q[, rs_algorithm := "fb_q"] # add a column for respiration algorithm
+
+##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# open all files and initialize matrix
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/vb_q_t_n-_BR-Sa1_I1850CLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365     #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365         #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365        #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365          #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 310, by = 10) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_q <- as.data.table(ar2dt) # convert array to data table
+vb_q <- rename(vb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_q <- vb_q[, rs_algorithm := "vb_q"] # add a column for respiration algorithm
+
+##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# open all files and initialize matrix
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/vb_e_t_n-_BR-Sa1_I1850CLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365     #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365          #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365         #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365          #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 310, by = 10) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_e <- as.data.table(ar2dt) # convert array to data table
+vb_e <- rename(vb_e, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+
+###############################################################################
+
+#####TRANSIENT##############################################################################
+# # set working directory, read in single file,and place variable names into "vars"
+# # use the "vars" variable to find what you are looking for
+# workdir <- '~/tmp/run/acme_t_n-_BR-Sa1_I20TRCLM45CN/run'
+# setwd(workdir)
+# ncfname <- 'acme_t_n-_BR-Sa1_I20TRCLM45CN.clm2.h0.1933-01-01-00000.nc'
+# ncf <- open.nc(ncfname)
+# vars <- as.list(read.nc(ncf))
+# ncf_data <- read.nc(ncfile = ncf)
+# #ncf_data$LEAF_MR
+# #ncf_data$XSMRPOOL
+# #ncf_data$GPP
+# #ncf_data$NPP
+# 
+# ##### build set of datatables, once for each variable of interest
+# 
+# ##### equation 1, create datatable from standard ACME maintenence respiration files
+# # open all files and initialize matrix
+# files <- list.files(path = '~/tmp/run/acme_t_n-_BR-Sa1_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+#                     full.names = TRUE)
+# 
+# ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+# 
+# # loop over all files
+# for(i in 1:length(files)) {
+#   ncf <- open.nc(files[i]) # open netcdf 
+#   a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+#   b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+#   c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+#   d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+#   ar[i, , ] <- rbind(a, b, c, d)
+#   #if(i%%50==0){print(i)}            
+#   #print(m) # for debugging
+#   close.nc(ncf) # close netcdf after reading
+# }
+# 
+# ar2d <- apply(ar, c(1, 2), sum)
+# t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+# ar2dt <- cbind(ar2d, t)
+# 
+# acme <- as.data.table(ar2dt) # convert array to data table
+# acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+# acme <- acme[, rs_algorithm := "acme"] # add a column for respiration algorithm
+# 
+# ##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# # open all files and initialize matrix
+# files <- list.files(path = '~/tmp/run/fb_q_t_n-_BR-Sa1_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+#                     full.names = TRUE)
+# 
+# ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+# 
+# # loop over all files
+# for(i in 1:length(files)) {
+#   ncf <- open.nc(files[i]) # open netcdf 
+#   a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+#   b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+#   c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+#   d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+#   ar[i, , ] <- rbind(a, b, c, d)
+#   #if(i%%50==0){print(i)}            
+#   #print(m) # for debugging
+#   close.nc(ncf) # close netcdf after reading
+# }
+# 
+# ar2d <- apply(ar, c(1, 2), sum)
+# t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+# ar2dt <- cbind(ar2d, t)
+# 
+# fb_q <- as.data.table(ar2dt) # convert array to data table
+# fb_q <- rename(fb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+# fb_q <- fb_q[, rs_algorithm := "fb_q"] # add a column for respiration algorithm
+# 
+# ##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# # open all files and initialize matrix
+# files <- list.files(path = '~/tmp/run/vb_q_t_n-_BR-Sa1_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+#                     full.names = TRUE)
+# 
+# ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+# 
+# # loop over all files
+# for(i in 1:length(files)) {
+#   ncf <- open.nc(files[i]) # open netcdf 
+#   a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+#   b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+#   c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+#   d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+#   ar[i, , ] <- rbind(a, b, c, d)
+#   #if(i%%50==0){print(i)}            
+#   #print(m) # for debugging
+#   close.nc(ncf) # close netcdf after reading
+# }
+# 
+# ar2d <- apply(ar, c(1, 2), sum)
+# t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+# ar2dt <- cbind(ar2d, t)
+# 
+# vb_q <- as.data.table(ar2dt) # convert array to data table
+# vb_q <- rename(vb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+# vb_q <- vb_q[, rs_algorithm := "vb_q"] # add a column for respiration algorithm
+# 
+# ##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# # open all files and initialize matrix
+# files <- list.files(path = '~/tmp/run/vb_e_t_n-_BR-Sa1_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+#                     full.names = TRUE)
+# 
+# ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+# 
+# # loop over all files
+# for(i in 1:length(files)) {
+#   ncf <- open.nc(files[i]) # open netcdf 
+#   a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+#   b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+#   c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+#   d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+#   ar[i, , ] <- rbind(a, b, c, d)
+#   #if(i%%50==0){print(i)}            
+#   #print(m) # for debugging
+#   close.nc(ncf) # close netcdf after reading
+# }
+# 
+# ar2d <- apply(ar, c(1, 2), sum)
+# t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+# ar2dt <- cbind(ar2d, t)
+# 
+# vb_e <- as.data.table(ar2dt) # convert array to data table
+# vb_e <- rename(vb_e, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+# vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+###############################################################################
+
+
+### MAKE FANCY FIGURES
+
+### concatinate the data tables into single data source then draw figures
+l <- list(acme, fb_q, vb_q, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+p <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "red", "green", "blue"))
+
+#, limits = c(1990, 2005)
+
+lmr <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr +  scale_x_continuous(name = '')                                                     # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+mr <- ggplot(dt, aes(x = YEAR, y = MR, color = rs_algorithm)) + geom_line() 
+mr <- mr + scale_x_continuous(name = '')
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")")) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+gpp <- ggplot(dt, aes(x = YEAR, y = GPP, color = rs_algorithm)) + geom_line()  
+gpp <- gpp + scale_x_continuous(name = 'Years')        # set x-axis labels to every 100 years
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+npp <- ggplot(dt, aes(x = YEAR, y = NPP, color = rs_algorithm)) + geom_line()  
+npp <- npp + scale_x_continuous(name = 'Years')       # set x-axis labels to every 100 years
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), 
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Four Respiration Equations")
+
diff --git a/ncplot_br_sa1_annula_test_t-.R b/ncplot_br_sa1_annula_test_t-.R
new file mode 100644
index 0000000..ab2713e
--- /dev/null
+++ b/ncplot_br_sa1_annula_test_t-.R
@@ -0,0 +1,137 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+########################################################################################
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/run/vb_e_t_n-_BR-Sa1_I20TRCLM45CN/run'
+setwd(workdir)
+ncfname <- 'vb_e_t_n-_BR-Sa1_I20TRCLM45CN.clm2.h0.2010-01-01-00000.nc'
+ncf <- open.nc(ncfname)
+vars <- as.list(read.nc(ncf))
+ncf_data <- read.nc(ncfile = ncf)
+#ncf_data$LEAF_MR
+#ncf_data$XSMRPOOL
+#ncf_data$GPP
+#ncf_data$NPP
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from ve_e maintenence respiration files and 
+##### standard temperature forcing 
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/vb_e_t_n-_BR-Sa1_I20TRCLM45CN/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365     #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365          #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365         #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365          #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 652, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_e_t <- as.data.table(ar2dt) # convert array to data table
+vb_e_t <- rename(vb_e_t, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_e_t <- vb_e_t[, rs_algorithm := "vb_e_t"] # add a column for respiration algorithm
+
+##### equation 2, create datatable from ve_e maintenence respiration files and 
+##### -10°C temperature forcing
+# open all files and initialize matrix
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/vb_e_t-_n-_BR-Sa1_I20TRCLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365    #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365         #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365        #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365         #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 310, by = 10) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_e_tminus10 <- as.data.table(ar2dt) # convert array to data table
+vb_e_tminus10 <- rename(vb_e_tminus10, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_e_tminus10 <- vb_e_tminus10[, rs_algorithm := "vb_e_tminus10"] # add a column for respiration algorithm
+
+
+### MAKE FANCY FIGURES
+
+### concatinate the data tables into single data source then draw figures
+l <- list(acme, fb_q, vb_q, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+p <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "red", "green", "blue"))
+
+#, limits = c(1990, 2005)
+
+lmr <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr +  scale_x_continuous(name = '')                                                     # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+mr <- ggplot(dt, aes(x = YEAR, y = MR, color = rs_algorithm)) + geom_line() 
+mr <- mr + scale_x_continuous(name = '')
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")")) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+gpp <- ggplot(dt, aes(x = YEAR, y = GPP, color = rs_algorithm)) + geom_line()  
+gpp <- gpp + scale_x_continuous(name = 'Years')        # set x-axis labels to every 100 years
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+npp <- ggplot(dt, aes(x = YEAR, y = NPP, color = rs_algorithm)) + geom_line()  
+npp <- npp + scale_x_continuous(name = 'Years')       # set x-axis labels to every 100 years
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), 
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Four Respiration Equations")
+
diff --git a/ncplot_ca-man_annual.R b/ncplot_ca-man_annual.R
new file mode 100644
index 0000000..b3ea073
--- /dev/null
+++ b/ncplot_ca-man_annual.R
@@ -0,0 +1,322 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+#####SPINUP##############################################################################
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/run/acme_t_n-_CA-MAN_I1850CLM45CN_ad_spinup/run'
+setwd(workdir)
+ncfname <- 'acme_t_n-_CA-MAN_I1850CLM45CN_ad_spinup.clm2.h0.0241-01-01-00000.nc'
+ncf <- open.nc(ncfname)
+vars <- as.list(read.nc(ncf))
+ncf_data <- read.nc(ncfile = ncf)
+#ncf_data$LEAF_MR
+#ncf_data$XSMRPOOL
+#ncf_data$GPP
+#ncf_data$NPP
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from standard ACME maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/acme_t_n-_CA-MAN_I1850CLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365     #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365          #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365         #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365          #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 301, by = 15) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+acme <- as.data.table(ar2dt) # convert array to data table
+acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+acme <- acme[, rs_algorithm := "acme"] # add a column for respiration algorithm
+
+##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# open all files and initialize matrix
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/fb_q_t_n-_CA-MAN_I1850CLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365    #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365         #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365        #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365         #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 301, by = 15) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+fb_q <- as.data.table(ar2dt) # convert array to data table
+fb_q <- rename(fb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+fb_q <- fb_q[, rs_algorithm := "fb_q"] # add a column for respiration algorithm
+
+##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# open all files and initialize matrix
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/vb_q_t_n-_CA-MAN_I1850CLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365     #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365         #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365        #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365          #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 301, by = 15) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_q <- as.data.table(ar2dt) # convert array to data table
+vb_q <- rename(vb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_q <- vb_q[, rs_algorithm := "vb_q"] # add a column for respiration algorithm
+
+##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# open all files and initialize matrix
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/vb_e_t_n-_CA-MAN_I1850CLM45CN_ad_spinup/run', pattern = '*clm2.h0.*.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24 * 365     #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24 * 365          #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24 * 365         #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24 * 365          #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1, 301, by = 15) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_e <- as.data.table(ar2dt) # convert array to data table
+vb_e <- rename(vb_e, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+
+###############################################################################
+
+# #####TRANSIENT##############################################################################
+# 
+# # set working directory, read in single file,and place variable names into "vars"
+# # use the "vars" variable to find what you are looking for
+# workdir <- '~/tmp/run/acme_t_n-_CA-MAN_I20TRCLM45CN/run'
+# setwd(workdir)
+# ncfname <- 'acme_t_n-_CA-MAN_I20TRCLM45CN.clm2.h0.1933-01-01-00000.nc'
+# ncf <- open.nc(ncfname)
+# vars <- as.list(read.nc(ncf))
+# ncf_data <- read.nc(ncfile = ncf)
+# #ncf_data$LEAF_MR
+# #ncf_data$XSMRPOOL
+# #ncf_data$GPP
+# #ncf_data$NPP
+# 
+# ##### build set of datatables, once for each variable of interest
+# 
+# ##### equation 1, create datatable from standard ACME maintenence respiration files
+# # open all files and initialize matrix
+# files <- list.files(path = '~/tmp/run/acme_t_n-_CA-MAN_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+#                     full.names = TRUE)
+# 
+# ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+# 
+# # loop over all files
+# for(i in 1:length(files)) {
+#   ncf <- open.nc(files[i]) # open netcdf 
+#   a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+#   b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+#   c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+#   d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+#   ar[i, , ] <- rbind(a, b, c, d)
+#   #if(i%%50==0){print(i)}            
+#   #print(m) # for debugging
+#   close.nc(ncf) # close netcdf after reading
+# }
+# 
+# ar2d <- apply(ar, c(1, 2), sum)
+# t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+# ar2dt <- cbind(ar2d, t)
+# 
+# acme <- as.data.table(ar2dt) # convert array to data table
+# acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+# acme <- acme[, rs_algorithm := "acme"] # add a column for respiration algorithm
+# 
+# ##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# # open all files and initialize matrix
+# files <- list.files(path = '~/tmp/run/fb_q_t_n-_CA-MAN_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+#                     full.names = TRUE)
+# 
+# ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+# 
+# # loop over all files
+# for(i in 1:length(files)) {
+#   ncf <- open.nc(files[i]) # open netcdf 
+#   a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+#   b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+#   c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+#   d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+#   ar[i, , ] <- rbind(a, b, c, d)
+#   #if(i%%50==0){print(i)}            
+#   #print(m) # for debugging
+#   close.nc(ncf) # close netcdf after reading
+# }
+# 
+# ar2d <- apply(ar, c(1, 2), sum)
+# t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+# ar2dt <- cbind(ar2d, t)
+# 
+# fb_q <- as.data.table(ar2dt) # convert array to data table
+# fb_q <- rename(fb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+# fb_q <- fb_q[, rs_algorithm := "fb_q"] # add a column for respiration algorithm
+# 
+# ##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# # open all files and initialize matrix
+# files <- list.files(path = '~/tmp/run/vb_q_t_n-_CA-MAN_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+#                     full.names = TRUE)
+# 
+# ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+# 
+# # loop over all files
+# for(i in 1:length(files)) {
+#   ncf <- open.nc(files[i]) # open netcdf 
+#   a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+#   b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+#   c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+#   d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+#   ar[i, , ] <- rbind(a, b, c, d)
+#   #if(i%%50==0){print(i)}            
+#   #print(m) # for debugging
+#   close.nc(ncf) # close netcdf after reading
+# }
+# 
+# ar2d <- apply(ar, c(1, 2), sum)
+# t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+# ar2dt <- cbind(ar2d, t)
+# 
+# vb_q <- as.data.table(ar2dt) # convert array to data table
+# vb_q <- rename(vb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+# vb_q <- vb_q[, rs_algorithm := "vb_q"] # add a column for respiration algorithm
+# 
+# ##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# # open all files and initialize matrix
+# files <- list.files(path = '~/tmp/run/vb_e_t_n-_CA-MAN_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+#                     full.names = TRUE)
+# 
+# ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+# 
+# # loop over all files
+# for(i in 1:length(files)) {
+#   ncf <- open.nc(files[i]) # open netcdf 
+#   a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+#   b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+#   c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+#   d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+#   ar[i, , ] <- rbind(a, b, c, d)
+#   #if(i%%50==0){print(i)}            
+#   #print(m) # for debugging
+#   close.nc(ncf) # close netcdf after reading
+# }
+# 
+# ar2d <- apply(ar, c(1, 2), sum)
+# t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+# ar2dt <- cbind(ar2d, t)
+# 
+# vb_e <- as.data.table(ar2dt) # convert array to data table
+# vb_e <- rename(vb_e, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+# vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+###############################################################################
+
+### MAKE FIGURES
+
+### concatinate the data tables into single data source then draw figures
+l <- list(acme, fb_q, vb_q, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+p <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "red", "green", "blue"))
+
+lmr <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr +  scale_x_continuous(name = '')                                                     # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+mr <- ggplot(dt, aes(x = YEAR, y = MR, color = rs_algorithm)) + geom_line() 
+mr <- mr + scale_x_continuous(name = '')
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")")) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+gpp <- ggplot(dt, aes(x = YEAR, y = GPP, color = rs_algorithm)) + geom_line()  
+gpp <- gpp + scale_x_continuous(name = 'Years')        # set x-axis labels to every 100 years
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+npp <- ggplot(dt, aes(x = YEAR, y = NPP, color = rs_algorithm)) + geom_line()  
+npp <- npp + scale_x_continuous(name = 'Years')       # set x-axis labels to every 100 years
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), 
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Four Respiration Equations")
+
diff --git a/ncplot_daily.R b/ncplot_daily.R
new file mode 100644
index 0000000..29a65e4
--- /dev/null
+++ b/ncplot_daily.R
@@ -0,0 +1,207 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/run/acme_t_n-_US-UMB_I20TRCLM45CN/run'
+setwd(workdir)
+ncfname <- 'acme_t_n-_US-UMB_I20TRCLM45CN.clm2.h0.1982-01-01-00000.nc'
+ncf <- open.nc(ncfname)
+vars <- as.list(read.nc(ncf))
+ncf_data <- read.nc(ncfile = ncf)
+#ncf_data$LEAF_MR
+#ncf_data$XSMRPOOL
+#ncf_data$GPP
+#ncf_data$NPP
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from standard ACME maintenence respiration files
+# open all files and initialize matrix
+
+#[1-2]{1}[0-9]{3}(?<!2007)
+
+files <- list.files(path = '~/tmp/run/acme_t_n-_US-UMB_I20TRCLM45CN/run', pattern = '*.h0.[1-2][0-9][0-9][0-9]-01-01-00000.nc', 
+                    full.names = TRUE)
+
+#files <- list.files(path = '~/tmp/us-umb/acme_t_n-/lnd/hist', pattern = '*US-UMB_I20TRCLM45CN.clm2.h0.*.nc', full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 5, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+#m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR")                  #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR")                       #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP")                      #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP")                      #gC/m^2/s
+  t <- var.get.nc(nc = ncf, "mcdate")                   #yyyymmdd
+  ar[i, , ] <- rbind(a, b, c, d, t)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- Reduce('rbind', lapply(1:365, function(k) ar[, , k])) 
+
+#t <- seq(0, 800, by = 1) # create time variable, step by 20
+acme <- as.data.table(ar2d) # convert array to data table
+#acme[, "V5"] <- t # overwrite "time" with t, and rename columns
+acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="DATE"))
+acme <- acme[, rs_algorithm := "ACME"] # add a column for respiration algorithm
+acme <- acme[order(DATE), ]
+
+##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/fb_q_t_n-/lnd/hist', pattern = '*[1-2][0-9][0-9][0-9]-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 5, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+#m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR")                       #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  t <- var.get.nc(nc = ncf, "mcdate")                   #yyyymmdd
+  ar[i, , ] <- rbind(a, b, c, d, t)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- Reduce('rbind', lapply(1:365, function(k) ar[, , k])) 
+
+#t <- seq(0, 800, by = 1) # create time variable, step by 20
+fb_q <- as.data.table(ar2d) # convert array to data table
+#acme[, "V5"] <- t # overwrite "time" with t, and rename columns
+fb_q <- rename(fb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="DATE"))
+fb_q <- fb_q[, rs_algorithm := "fb_q"] # add a column for respiration algorithm
+fb_q <- fb_q[order(DATE), ]
+
+##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/vb_q_t_n-/lnd/hist', pattern = '*[1-2][0-9][0-9][0-9]-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 5, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+#m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR")                       #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  t <- var.get.nc(nc = ncf, "mcdate")                   #yyyymmdd
+  ar[i, , ] <- rbind(a, b, c, d, t)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- Reduce('rbind', lapply(1:365, function(k) ar[, , k])) 
+
+#t <- seq(0, 800, by = 1) # create time variable, step by 20
+vb_q <- as.data.table(ar2d) # convert array to data table
+#acme[, "V5"] <- t # overwrite "time" with t, and rename columns
+vb_q <- rename(vb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="DATE"))
+vb_q <- vb_q[, rs_algorithm := "vb_q"] # add a column for respiration algorithm
+vb_q <- vb_q[order(DATE), ]
+
+##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/vb_e_t_n-/lnd/hist', pattern = '*[1-2][0-9][0-9][0-9]-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 5, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+#m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR")                       #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  t <- var.get.nc(nc = ncf, "mcdate")                   #yyyymmdd
+  ar[i, , ] <- rbind(a, b, c, d, t)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- Reduce('rbind', lapply(1:365, function(k) ar[, , k])) 
+
+#t <- seq(0, 800, by = 1) # create time variable, step by 20
+vb_e <- as.data.table(ar2d) # convert array to data table
+#acme[, "V5"] <- t # overwrite "time" with t, and rename columns
+vb_e <- rename(vb_e, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="DATE"))
+vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+vb_e <- vb_e[order(DATE), ]
+### concatinate the data tables into single data source
+l <- list(acme, fb_q, vb_q, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+### MAKE FANCY FIGURES
+
+p <- ggplot(dt, aes(x = DATE, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "red", "green", "blue"))
+
+lmr <- ggplot(dt, aes(x = DATE, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr +  scale_x_continuous(name = '')                                                     # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"), limits = c(0, 400))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "red", "green", "blue"))
+                               
+mr <- ggplot(dt, aes(x = DATE, y = MR, color = rs_algorithm)) + geom_line() 
+mr <- mr + scale_x_continuous(name = '')
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")"), limits = c(0, 400)) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+gpp <- ggplot(dt, aes(x = DATE, y = GPP, color = rs_algorithm)) + geom_line()  
+gpp <- gpp + scale_x_continuous(name = 'Years', breaks = c(100, 200, 300, 400, 500, 600, 700, 800), 
+                     labels = c("100", "200", "300", "400", "500", "600", "700", "800"))        # set x-axis labels to every 100 years
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+npp <- ggplot(dt, aes(x = DATE, y = NPP, color = rs_algorithm)) + geom_line()  
+npp <- npp + scale_x_continuous(name = 'Years', breaks = c(100, 200, 300, 400, 500, 600, 700, 800), 
+                      labels = c("100", "200", "300", "400", "500", "600", "700", "800"))       # set x-axis labels to every 100 years
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), 
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Four Respiration Equations")
diff --git a/ncplot_gpp_comp.R b/ncplot_gpp_comp.R
new file mode 100644
index 0000000..02270e6
--- /dev/null
+++ b/ncplot_gpp_comp.R
@@ -0,0 +1,174 @@
+# Load some libraries
+#library(rPython)
+#library(RNetCDF)
+library(ncdf4)
+library(data.table)
+library(plyr)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+library(zoo)
+library(scales)
+
+##### Read in csv data from fluxnet directory on google drive
+
+csvfiles <- list.files(path = '~/Google Drive/globland/data/fluxnet/AMF_US-UMB_BASE-BADM_10-1', pattern = 'AMF_USUMB_[2][0][1][0-4]_L2_GF_V010.csv', 
+                        full.names = TRUE)
+temp <- lapply(csvfiles, fread, skip = 20)
+data <- rbindlist(temp)
+cols = paste("V", c(1, 4, 39), sep="")
+tgpp <- data[, cols, with = F]
+tgpp[tgpp == -9999] <- NA
+tgpp <- rename(tgpp, c("V1"="YEAR","V4"="DOY", "V39"="GPP.umols"))
+tgpp[,`:=`(GPP.grams = GPP.umols * 3600 * 12 * 1E-6)] # convert umols CO2m^-2s^-1 to g C^-2s^-1
+key <- c("YEAR", "DOY")
+daily.tgpp <- tgpp[ , lapply(.SD, sum), by = key] # aggrate 24 hourly gpp values to a daily sum
+
+
+
+
+##### Read in netcdf data from ACME
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/run/acme_t_n-_US-UMB_I20TRCLM45CN/run'
+setwd(workdir)
+ncfname <- 'acme_t_n-_US-UMB_I20TRCLM45CN.clm2.h0.2000-01-01-00000.nc'
+#vars <- as.list(read.nc(ncf))
+#ncf_data <- read.nc(ncfile = ncf)
+
+# read individual netcdf variables into arrays
+nc <- nc_open(ncfname)
+
+lmr <- ncvar_get(nc, "LEAF_MR")
+mr <- ncvar_get(nc, "MR")
+gpp <- ncvar_get(nc, "GPP")
+npp <- ncvar_get(nc, "NPP")
+mcdate <- ncvar_get(nc, "mcdate") #ncdf4
+
+nc_close(ncfname)
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from standard ACME maintenence respiration files
+# open all files and initialize matrix
+
+ncfiles <- list.files(path = '~/tmp/run/acme_t_n-_US-UMB_I20TRCLM45CN/run', pattern = '*.clm2.h0.[2][0][1][0-4]-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(ncfiles), 5, length(gpp))) # create an array of the proper dimensions and fill it with NAs
+
+# loop over all files
+for(i in 1:length(ncfiles)) {
+  nc <- nc_open(ncfiles[i])
+  a <- ncvar_get(nc, "mcdate")
+  b <- ncvar_get(nc, "LEAF_MR") * 3600 * 24
+  c <- ncvar_get(nc, "MR") * 3600 * 24
+  d <- ncvar_get(nc, "GPP") * 3600 * 24
+  e <- ncvar_get(nc, "NPP") * 3600 * 24
+  ar[i, , ] <- rbind(a, b, c, d, e)
+  # print(ar) # for debugging
+  # print(i)
+  nc_close(nc)  
+}
+
+ar2d <- (Reduce('rbind', lapply(1:365, function(k) ar[, , k]))) 
+
+#t <- seq(0, 800, by = 20) # create time variable, step by 20
+acme <- as.data.table(ar2d) # convert matrix to data table
+#acme[, "V5"] <- t # overwrite "time" with t, and rename columns
+acme <- rename(acme, c("V1"="DATE", "V2"="LEAF_MR","V3"="MR", "V4"="GPP", "V5"="NPP"))
+acme <- acme[, rs_algorithm := "ACME"] # add a column for respiration algorithm
+acme <- acme[order(DATE)] # order data by date
+acme$DATE <- as.Date(as.character(acme$DATE), "%Y%m%d")
+
+
+# sum(acme$DATE == 20150101) # looks like 365 obs for 2015, removed odd 20150101 netcdf from file list
+# sum(acme$DATE == 20140101) # only 1 obs for everything else
+
+
+##### equation 2, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# open all files and initialize matrix
+
+ncfiles <- list.ncfiles(path = '~/tmp/run/vb_e_t_n-_US-UMB_I20TRCLM45CN/run', pattern = '*.clm2.h0.[2][0][1][0-4]-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(ncfiles), 5, length(gpp))) # create an array of the proper dimensions and fill it with NAs
+
+# loop over all files
+for(i in 1:length(ncfiles)) {
+  nc <- nc_open(ncfiles[i])
+  a <- ncvar_get(nc, "mcdate")
+  b <- ncvar_get(nc, "LEAF_MR") * 3600 * 24
+  c <- ncvar_get(nc, "MR") * 3600 * 24
+  d <- ncvar_get(nc, "GPP") * 3600 * 24
+  e <- ncvar_get(nc, "NPP") * 3600 * 24 
+  ar[i, , ] <- rbind(a, b, c, d, e)
+  # print(ar) # for debugging
+  # print(i)
+  nc_close(nc)  
+}
+
+ar2d <- (Reduce('rbind', lapply(1:365, function(k) ar[, , k]))) 
+# ar2dtest <- (Reduce('rbind', lapply(1:165, function(k) ar[k, , ])))
+
+#t <- seq(0, 800, by = 20) # create time variable, step by 20
+vb_e <- as.data.table(ar2d) # convert matrix to data table
+#vb_e[, "V5"] <- t # overwrite "time" with t, and rename columns
+vb_e <- rename(vb_e, c("V1"="DATE", "V2"="LEAF_MR","V3"="MR", "V4"="GPP", "V5"="NPP"))
+vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+vb_e <- vb_e[order(DATE)] # order data by date
+vb_e$DATE <- as.Date(as.character(vb_e$DATE), "%Y%m%d")
+
+### concatinate the data tables into single data source
+
+l <- list(acme, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+### MAKE FANCY FIGURES
+
+p <- ggplot(dt, aes(x = DATE, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "green"))
+
+lmr <- ggplot(dt, aes(x = DATE, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr + scale_x_date(labels = date_format("%Y"), name = '', limits = as.Date(c("2010-01-01","2014-12-31")))#, limits = as.Date(c("2000-01-01","2012-01-01"))) 
+#lmr <- lmr +  scale_x_continuous(name = '')    #limits = c(2010, 2012),                                                  # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~d^-1*")"), limits = c(0, 5))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "green"))
+                               
+mr <- ggplot(dt, aes(x = DATE, y = MR, color = rs_algorithm)) + geom_line()
+mr <- mr  + scale_x_date(labels = date_format("%Y"), name = '', limits = as.Date(c("2010-01-01","2014-12-31")))#, limits = as.Date(c("2000-01-01","2012-01-01"))) 
+#mr <- mr + scale_x_continuous(name = '') 
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~d^-1*")"), limits = c(0, 5)) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "green"))
+
+gpp <- ggplot(dt, aes(x = DATE, y = GPP, color = rs_algorithm)) + geom_line()
+gpp <- gpp  + scale_x_date(labels = date_format("%Y"), name = 'Year', limits = as.Date(c("2010-01-01","2014-12-31")))#, limits = as.Date(c("2000-01-01","2012-01-01"))) 
+#gpp <- gpp + scale_x_continuous(name = 'Years') #limits = c(2010, 2012), 
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~d^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "green"))
+
+npp <- ggplot(dt, aes(x = DATE, y = NPP, color = rs_algorithm)) + geom_line()
+npp <- npp  + scale_x_date(labels = date_format("%Y"), name = 'Year', limits = as.Date(c("2010-01-01","2014-12-31")))#, limits = as.Date(c("2000-01-01","2012-01-01"))) 
+#npp <- npp + scale_x_continuous(name = 'Years') #limits = c("2010", "2012"), 
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~d^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "green"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), #, origin = "1970-01-01"
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Two Respiration Equations")
diff --git a/ncplot_transient.R b/ncplot_transient.R
new file mode 100644
index 0000000..3d00625
--- /dev/null
+++ b/ncplot_transient.R
@@ -0,0 +1,179 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/us-umb/acme_t_n-/lnd/hist'
+setwd(workdir)
+ncfname <- 'acme_t_n-_US-UMB_I20TRCLM45CN.clm2.h0.1982-01-01-00000.nc'
+ncf <- open.nc(ncfname)
+vars <- as.list(read.nc(ncf))
+ncf_data <- read.nc(ncfile = ncf)
+#ncf_data$LEAF_MR
+#ncf_data$XSMRPOOL
+#ncf_data$GPP
+#ncf_data$NPP
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from standard ACME maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/acme_t_n-/lnd/hist', pattern = '*[1-2][0-9][0-9][0-9]-01-01-00000.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- sum(var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24)      #gC/m^2/s
+  b <- sum(var.get.nc(nc = ncf, "MR") * 3600 * 24)                 #gC/m^2
+  c <- sum(var.get.nc(nc = ncf, "GPP") * 3600 * 24)          #gC/m^2/s
+  d <- sum(var.get.nc(nc = ncf, "NPP") * 3600 *24)           #gC/m^2/s
+  t <- sum(var.get.nc(nc = ncf, "mcdate"))                   #yyyymmdd
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+t <- seq(1850, 2007, by = 1) # create time variable, step by 20
+acme <- as.data.table(m) # convert matrix to data table
+acme[, "V5"] <- t # overwrite "time" with t, and rename columns
+acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="YEAR"))
+acme <- acme[, rs_algorithm := "acme"] # add a column for respiration algorithm
+
+##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/fb_q_t_n-/lnd/hist', pattern = '*[1-2][0-9][0-9][0-9]-01-01-00000.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- sum(var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24)      #gC/m^2/s
+  b <- sum(var.get.nc(nc = ncf, "MR") * 3600 *24)                 #gC/m^2
+  c <- sum(var.get.nc(nc = ncf, "GPP") * 3600 * 24)          #gC/m^2/s
+  d <- sum(var.get.nc(nc = ncf, "NPP") * 3600 *24)           #gC/m^2/s
+  t <- sum(var.get.nc(nc = ncf, "mcdate"))                   #yyyymmdd
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+t <- seq(1850, 2007, by = 1) # create time variable, step by 20
+fb_q <- as.data.table(m) # convert matrix to data table
+fb_q[, "V5"] <- t # overwrite "time" with t, and rename columns
+fb_q <- rename(fb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="YEAR"))
+fb_q <- fb_q[, rs_algorithm := "fb_q"] # add a column for respiration algorithm
+
+##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/vb_q_t_n-/lnd/hist', pattern = '*[1-2][0-9][0-9][0-9]-01-01-00000.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- sum(var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24)      #gC/m^2/s
+  b <- sum(var.get.nc(nc = ncf, "MR") * 3600 * 24)                 #gC/m^2
+  c <- sum(var.get.nc(nc = ncf, "GPP") * 3600 * 24)          #gC/m^2/s
+  d <- sum(var.get.nc(nc = ncf, "NPP") * 3600 *24)           #gC/m^2/s
+  t <- sum(var.get.nc(nc = ncf, "mcdate"))                   #yyyymmdd
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+t <- seq(1850, 2007, by = 1) # create time variable, step by 20
+vb_q <- as.data.table(m) # convert matrix to data table
+vb_q[, "V5"] <- t # overwrite "time" with t, and rename columns
+vb_q <- rename(vb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="YEAR"))
+vb_q <- vb_q[, rs_algorithm := "vb_q"] # add a column for respiration algorithm
+
+##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/us-umb/vb_e_t_n-/lnd/hist', pattern = '*[1-2][0-9][0-9][0-9]-01-01-00000.nc', 
+                    full.names = TRUE)
+m <- matrix(NA, length(files), 5)
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- sum(var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24)      #gC/m^2/s
+  b <- sum(var.get.nc(nc = ncf, "MR") * 3600 * 24)                 #gC/m^2
+  c <- sum(var.get.nc(nc = ncf, "GPP") * 3600 * 24)          #gC/m^2/s
+  d <- sum(var.get.nc(nc = ncf, "NPP") * 3600 *24)           #gC/m^2/s
+  t <- sum(var.get.nc(nc = ncf, "mcdate"))                   #yyyymmdd
+  m[i, ] <- t(rbind(a, b, c, d, t))
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+t <- seq(1850, 2007, by = 1) # create time variable, step by 20
+vb_e <- as.data.table(m) # convert matrix to data table
+vb_e[, "V5"] <- t # overwrite "time" with t, and rename columns
+vb_e <- rename(vb_e, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "V5"="YEAR"))
+vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+
+### concatinate the data tables into single data source
+l <- list(acme, fb_q, vb_q, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+### MAKE FANCY FIGURES
+
+p <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "red", "green", "blue"))
+
+lmr <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr +  scale_x_continuous(name = '')                                                     # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"), limits = c(0, 400))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "red", "green", "blue"))
+                               
+mr <- ggplot(dt, aes(x = YEAR, y = MR, color = rs_algorithm)) + geom_line() 
+mr <- mr + scale_x_continuous(name = '')
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")"), limits = c(0, 400)) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+gpp <- ggplot(dt, aes(x = YEAR, y = GPP, color = rs_algorithm)) + geom_line()  
+gpp <- gpp + scale_x_continuous(name = 'YEAR', breaks = c(1850, 1900, 1950, 2000), 
+                     labels = c("1850", "1900", "1950", "2000"))        # set x-axis labels to every 50 years
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+npp <- ggplot(dt, aes(x = YEAR, y = NPP, color = rs_algorithm)) + geom_line()  
+npp <- npp + scale_x_continuous(name = 'YEAR', breaks = c(1850, 1900, 1950, 2000), 
+                      labels = c("1850", "1900", "1950", "2000"))       # set x-axis labels to every 50 years
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), 
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Four Respiration Equations")
diff --git a/ncplot_us-umb_annual.R b/ncplot_us-umb_annual.R
new file mode 100644
index 0000000..5e6d69b
--- /dev/null
+++ b/ncplot_us-umb_annual.R
@@ -0,0 +1,186 @@
+# Load some libraries
+#library(rPython)
+#library(ncdf4)
+library(data.table)
+library(plyr)
+library(RNetCDF)
+library(ggplot2)
+library(grid)
+library(gridExtra)
+
+# set working directory, read in single file,and place variable names into "vars"
+# use the "vars" variable to find what you are looking for
+workdir <- '~/tmp/run/acme_t_n-_US-UMB_I20TRCLM45CN/run'
+setwd(workdir)
+ncfname <- 'acme_t_n-_US-UMB_I20TRCLM45CN.clm2.h0.1933-01-01-00000.nc'
+ncf <- open.nc(ncfname)
+vars <- as.list(read.nc(ncf))
+ncf_data <- read.nc(ncfile = ncf)
+#ncf_data$LEAF_MR
+#ncf_data$XSMRPOOL
+#ncf_data$GPP
+#ncf_data$NPP
+
+##### build set of datatables, once for each variable of interest
+
+##### equation 1, create datatable from standard ACME maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/acme_t_n-_US-UMB_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+acme <- as.data.table(ar2dt) # convert array to data table
+acme <- rename(acme, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+acme <- acme[, rs_algorithm := "acme"] # add a column for respiration algorithm
+
+##### equation 2, create datatable from fixed base and Q10=2 respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/fb_q_t_n-_US-UMB_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+fb_q <- as.data.table(ar2dt) # convert array to data table
+fb_q <- rename(fb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+fb_q <- fb_q[, rs_algorithm := "fb_q"] # add a column for respiration algorithm
+
+##### equation 3, create datatable from variable base and Q10=2 maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/vb_q_t_n-_US-UMB_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_q <- as.data.table(ar2dt) # convert array to data table
+vb_q <- rename(vb_q, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_q <- vb_q[, rs_algorithm := "vb_q"] # add a column for respiration algorithm
+
+##### equation 4, create datatable from variable base and exponential (variable Q10) maintenence respiration files
+# open all files and initialize matrix
+files <- list.files(path = '~/tmp/run/vb_e_t_n-_US-UMB_I20TRCLM45CN/run', pattern = '*clm2.h0.[1-2]{1}[0-9]{3}-01-01-00000.nc', 
+                    full.names = TRUE)
+
+ar <- array(data = NA, c(length(files), 4, length(ncf_data$NPP))) # create an array of the proper length and fill it with NAs
+
+# loop over all files
+for(i in 1:length(files)) {
+  ncf <- open.nc(files[i]) # open netcdf 
+  a <- var.get.nc(nc = ncf, "LEAF_MR") * 3600 * 24      #gC/m^2/s
+  b <- var.get.nc(nc = ncf, "MR") * 3600 * 24           #gC/m^2
+  c <- var.get.nc(nc = ncf, "GPP") * 3600 * 24          #gC/m^2/s
+  d <- var.get.nc(nc = ncf, "NPP") * 3600 *24           #gC/m^2/s
+  ar[i, , ] <- rbind(a, b, c, d)
+  #if(i%%50==0){print(i)}            
+  #print(m) # for debugging
+  close.nc(ncf) # close netcdf after reading
+}
+
+ar2d <- apply(ar, c(1, 2), sum)
+t <- seq(1850, 2015, by = 1) # create yearly time variable, step by 1
+ar2dt <- cbind(ar2d, t)
+
+vb_e <- as.data.table(ar2dt) # convert array to data table
+vb_e <- rename(vb_e, c("V1"="LEAF_MR", "V2"="MR","V3"="GPP", "V4"="NPP", "t"="YEAR"))
+vb_e <- vb_e[, rs_algorithm := "vb_e"] # add a column for respiration algorithm
+
+### MAKE FANCY FIGURES
+
+### concatinate the data tables into single data source then draw figures
+l <- list(acme, fb_q, vb_q, vb_e)
+dt <- rbindlist(l, use.names=TRUE)
+
+p <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + 
+  geom_line() + theme(legend.position = "bottom") + 
+  scale_color_manual(values=c("black", "red", "green", "blue"))
+
+lmr <- ggplot(dt, aes(x = YEAR, y = LEAF_MR, color = rs_algorithm)) + geom_line()
+lmr <- lmr +  scale_x_continuous(name = '', limits = c(1990, 2005))                                                     # hide x-axis name
+lmr <- lmr +  scale_y_continuous(expression("LMR (g C m"^-2~yr^-1*")"))     # set y-axis to be 0:400 and label
+lmr <- lmr + theme(axis.text.x = element_blank(), axis.ticks = element_blank())                 # hide ledgend, x-axis tick marks
+lmr <- lmr + scale_color_manual(values=c("black", "red", "green", "blue"))
+                               
+mr <- ggplot(dt, aes(x = YEAR, y = MR, color = rs_algorithm)) + geom_line() 
+mr <- mr + scale_x_continuous(name = '', limits = c(1990, 2005))
+mr <- mr + scale_y_continuous(expression("MR (g C m"^-2~yr^-1*")")) 
+mr <- mr + theme(axis.text.x = element_blank(), axis.ticks = element_blank()) 
+mr <- mr + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+gpp <- ggplot(dt, aes(x = YEAR, y = GPP, color = rs_algorithm)) + geom_line()  
+gpp <- gpp + scale_x_continuous(name = 'Years', limits = c(1990, 2005))        # set x-axis labels to every 100 years
+gpp <- gpp + scale_y_continuous(expression("GPP (g C m"^-2~yr^-1*")")) 
+gpp <- gpp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+npp <- ggplot(dt, aes(x = YEAR, y = NPP, color = rs_algorithm)) + geom_line()  
+npp <- npp + scale_x_continuous(name = 'Years', limits = c(1990, 2005))       # set x-axis labels to every 100 years
+npp <- npp + scale_y_continuous(expression("NPP (g C m"^-2~yr^-1*")")) 
+npp <- npp + scale_color_manual(values=c("black", "red", "green", "blue"))
+
+# create custom function to plot single ledgend for all pannels 
+grid_arrange_legend<-function(a.gplot){
+  tmp <- ggplot_gtable(ggplot_build(a.gplot))
+  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
+  legend <- tmp$grobs[[leg]]
+  return(legend)
+}
+
+legend <- grid_arrange_legend(p)
+
+grid.arrange(arrangeGrob(lmr + theme_bw() + theme(legend.position = "none"), 
+                         mr + theme_bw() + theme(legend.position = "none"),
+                         gpp + theme_bw() + theme(legend.position = "none"), 
+                         npp + theme_bw() + theme(legend.position = "none"),
+                         legend, heights=c(0.45, 0.45, 0.1), nrow=3), 
+             top = "Four Respiration Equations")
+