---
title: "Supplementary Material to Eckes-Shephard et al 2020"
author: "Annemarie Eckes-Shephard"
date: "16/03/2020"
output: pdf_document
---

This document contains all analysis code used in the publication of Eckes-Shephard et al 2020.
The first section (Functions..) contains all homemade functions used in this publication.
The rest of the document is organised by the topics covered in the publication.
Prior to running this script, make sure all data necessary for the analysis is in a folder called data. 
See README for a description of the datasets used in this study.

# Functions used in the below script
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,cache=TRUE)
library(xts)
library(dplR)
library(zoo)
library(lubridate)
library(dplyr)
library(tseries)

library(Hmisc)
library(plsgenomics)
library(fields)
library(Metrics)
library(RColorBrewer)
library(readxl)
library(qmap) # for empirical quantile mapping for rainfall bias reduction

# library(CaviaR) only necessary, if you want to reproduce the entire procedure. The output from this package is saved and loaded into the script
library(mapproj)
library(gridExtra)
library(rworldmap)

# colour for all CRUNCEP forced simulations:
CRUNCEP_col <- "orange"
obs_colour <- "light green"
ds_col<-"black"
Spruce_col ="black"
Larch_col= "dark green"
  
#preparing for column-plot, to plot water content in middle of each simulated soil layer
layer<-c(-45,-91,-166,-289,-493,-829,-1383,-2296,-5000)
layer_upper<- c(0,layer)
layermiddle<-layer_upper[-10] - (layer_upper[-10]-layer)/2
ymax=0.4
```

```{r fct_get_soilMoisture}

 get_soilMoisture<-function(filename,start_year_want,info=TRUE,file_trunc=FALSE){
  
  # to speed up the read-in time for large files, determine where the file should be read from:
  d=365
  hhrs=48
 
  if(file_trunc==TRUE){# account for the possibility of the soilmoisture file being shortened, and therefore having no head
  first_year_file <- read.table(filename, nrows = 1, header = FALSE)
    first_year_file <- as.integer(first_year_file$V1[1])
  }else{
     header <- read.table(filename, nrows = 1, header = FALSE, stringsAsFactors = FALSE)
     first_year_file <- read.table(filename, nrows = 1, header = TRUE)
       first_year_file <- as.integer(first_year_file$jyear[1])
  }
    
   #check years wanted are in file
   if(start_year_want < first_year_file){
    print('no such date in soil moisture data. First year available is:', first_year_file)
   }
  
  year_read_in <- start_year_want - first_year_file
  line_read_in <-  year_read_in*d*hhrs 
  
  # read in from correct year:
  sm_sel    <- read.table(filename, skip = line_read_in+1, header = FALSE)
  
if(file_trunc==TRUE){# account for the possibility of the soilmoisture file being shortened, and therefore having no head
   last_year_file<-max(unique(sm_sel$V1))
}else{
  colnames( sm_sel ) <- unlist(header)
  last_year_file<-max(unique(sm_sel$jyear))
}

  if(info==TRUE){
  print(paste('soil moisture from was available from ',first_year_file,'to',last_year_file,sep=" "))
  print(paste('read in started from ',start_year_want,sep=" "))
  }
  
  #create timestamp for wanted duration:
  # avoid leapyear issues, so go year by year

  if(file_trunc==TRUE){# account for the possibility of the soilmoisture file being shortened, and therefore having no head
   ys <-unique(sm_sel$V1)
}else{
  ys <-unique(sm_sel$jyear)

}
  for (i in 1:length(ys)){
 
    if(file_trunc==TRUE){# account for the possibility of the soilmoisture file being shortened, and therefore having no head
subs <- sm_sel[sm_sel$V1==ys[i],]
}else{
   subs <- sm_sel[sm_sel$jyear==ys[i],]

}
  subs$TIMESTAMP <- seq(as.POSIXct(paste(ys[i],"-01-01 00:00:00", sep=""),tz="CET"),by= "30 min",length.out=dim(subs)[1])
  if(i==1){collect<-subs}else{collect<-rbind(collect,subs)}
  }
 
  
  # select relevant columns only:
  sm_sel_lrs<- as.data.frame(collect[5:13])
  return(sm_sel_lrs)
}

```

```{r fct_polylims}
# from https://lukemiller.org/index.php/2012/12/generating-polygon-boundaries-for-plotting-simple-time-series-data-with-missing-data/
polylims <- function(xdata, ydata, ydata2) {
	# A function to find the beginings and endings of each run of real numbers
	# in a set of data, in order to create polygons for plotting. The assumption
	# is that ydata, ydata2, and xdata vectors are the same length, and that 
	# the ydata and ydata2 vectors contain NA values for missing data. The 
	# same values in ydata and ydata2 must be missing. The output will be 
	# a list of data frames of x and y values suitable for plotting polygons.

	# Use rle function to find contiguous real numbers
	rl = rle(is.na(ydata))
	starts = vector()
	ends = vector()
	indx = 1
	for (i in 1:length(rl$lengths)){
		if (rl$values[i]){
			# Value was NA, advance index without saving the numeric values
			indx = indx + rl$lengths[i]
		} else {
			# Value was a real number, extract and save those values
			starts = c(starts,indx)
			ends = c(ends, (indx + rl$lengths[i] - 1))
			indx = indx + rl$lengths[i]
		}	
	}

	# At this point the lengths of the vectors 'starts' and 'ends' should be
	# equal, and each pair of values represents the starting and ending indices
	# of a continuous set of data in the ydata vector.

	# Next separate out each set of continuous ydata, and the associated xdata,
	# and format them for plotting as a polygon.
	polylist = list()
	for (i in 1:length(starts)){
		temp = data.frame(x = c(xdata[starts[i]],xdata[starts[i]:ends[i]],
						rev(xdata[starts[i]:ends[i]])),
				y = c(ydata[starts[i]],ydata2[starts[i]:ends[i]],
						rev(ydata[starts[i]:ends[i]])))
		polylist[[i]] = temp	
	}
	polylist
	# You can iterate through the items in polylist and plot them as 
	# polygons on your plot. Use code similar to the following:
	#	for (i in 1:length(polylist)){
	#			polygon(polylist[[i]]$x, polylist[[i]]$y, 
	#				col = rgb(0.5,0.5,0.5,0.5), border = NA)
	#	}
}


```

```{r fct_make_transparent}

#note: always pass alpha on the 0-255 scale
makeTransparent<-function(someColor, alpha=100)
{
  newColor<-col2rgb(someColor)
  apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
    blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}

```

```{r fct_aggregate_vars, echo=FALSE}
library(xts)
# create datetime with appropriate start day and with correct year. Does not create sub-daily timestamps.
# Regardless if the input is hourly or sub-hourly, it will create the appropriate daily timestamp for subsequent met-variable aggregation.
#' @param start_date, string of start date in this format: "2007-01-01""
#' @param HYBRID_output df with diagnostic output produced form HYBRID, needs columns 'DOY' and 'jyear'.

aggregate_to_daily_met_variables<-function(HYBRID_output,start_date){
dat<-HYBRID_output
dat$DOY<-dat$DOY-1
# subsetting needed to avoid leapyear issues:
ys<-unique(dat$jyear)
for ( i in 1:length(ys)){
  tmp<-dat[which(dat$jyear==ys[i]),]
  start_date<-as.Date(paste(ys[i],"-01-01",sep=""))
  tmp$TIME<-as.Date(tmp$DOY, origin = start_date)
  if(i==1){subs=tmp}else{subs<-rbind(subs,tmp)}
}

#subs$TIME<-as.Date(subs$DOY, origin = start_date)
#day <- as.numeric(strftime(subs$TIME, format = "%d"))
#month <- as.numeric(strftime(subs$TIME, format = "%m"))
#year <- subs$jyear
#subs$TIME<-as.Date(paste(year,month,day,sep="-"))
subs$jyear<-NULL
subs$DOY<-NULL

ns<-names(subs)
# remove hour and TIME variables
ns_short<-ns[2:(length(ns)-1)]

out<-matrix(0,ncol=length(ns),nrow=length(unique(subs$TIME)))
  subs.xts<-xts(subs[,c("rain"  , "tair"  , "lwdown", "swdown", "qair" ,  "press")],subs$TIME)
  out<-apply.daily(subs.xts,mean)

out$rain<-apply.daily(subs.xts$rain,sum)
return(out)
}

# operation can be either sum or mean
aggregate_to_daily_variables<-function(HYBRID_output,start_date,operation="sum"){
dat     <- HYBRID_output
dat$DOY <- dat$DOY-1

# subsetting needed to avoid leapyear issues:
ys<-unique(dat$year)
for ( i in 1:length(ys)){
  tmp<-dat[which(dat$year==ys[i]),]
  start_date<-paste(ys[i],"-01-01",sep="")
  tmp$date<-as.Date(tmp$DOY, origin = start_date)
  #tmp$date<-date(tmp$TIMESTAMP)
  tmp$date <- as.Date(paste(ys[i],tmp$month,tmp$DOY,sep="-"))
  if(i==1){subs=tmp}else{subs<-rbind(subs,tmp)}
}

subs$year  <- NULL
subs$DOY   <- NULL

ns  <- names(subs)
# remove hour and TIME variables
#ns_short      <- ns[1:8]#ns[2:(length(ns)-1)]

out                <- data.frame(matrix(0,ncol=length(ns),nrow=length(unique(subs$date))))
colnames(out)          <- c(ns)
subs.xts               <- xts(subs[,1:(length(ns)-1)],subs$date)


for(n in 1:(length(ns)-1)){
varname       <- paste(ns[n])
if(operation=="sum"){
out[,varname] <- apply.daily(subs.xts[,varname],sum)
}
if(operation=="mean"){
out[,varname] <- apply.daily(subs.xts[,varname],mean)
}
}
out           <- data.frame(out)
out$date      <- unique(subs$date)

return(out)
}

#var.monthly tibble with var.monthly variable in monthly resolution
#var.name variable name which should be used for the correlation
#dataframe with TRW column which should be added to.
# the number of years for which monthly data is avialble must be the same as the number of years for which TRW are available
correlate_TRW <- function(var.monthly,var.name,df.TRW,plot=FALSE){
  months<- c("ja","f","m","ap","my","j","jl","a","s","o","n","d","Ja","F","M","Ap","My","J","Jl","A","S","O","N","D")#"MJJA","JJA")
  months_plotting<- c("j","f","m","a","m","j","j","a","s","o","n","d","J","F","M","A","M","J","J","A","S","O","N","D")#,"MJJA","JJA")
  y.rm<-length(t(var.monthly[var.monthly$month==2,][var.name]))
    
  #sometimes a month doesn't contain a value, becuase all values were filtered out. Therefore, put the appropriate length into a pre-made matrix
  tmp           <- data.frame( matrix(NA,ncol = length(months),nrow = y.rm) )
  colnames(tmp) <- months

#move all months in appropriate colum for TRW correlation comparison, previous year
tmp$ja[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==1,][var.name])[,-y.rm])
tmp$f[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==2,][var.name])[,-y.rm])
tmp$m[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==3,][var.name])[,-y.rm])
tmp$ap[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==4,][var.name])[,-y.rm])
tmp$my[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==5,][var.name])[,-y.rm])
tmp$j[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==6,][var.name])[,-y.rm])
tmp$jl[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==7,][var.name])[,-y.rm])
tmp$a[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==8,][var.name])[,-y.rm])
tmp$s[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==9,][var.name])[,-y.rm])
tmp$o[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==10,][var.name])[,-y.rm])
tmp$n[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==11,][var.name])[,-y.rm])
tmp$d[1:length(t(var.monthly[var.monthly$month==1,][var.name]))] <- c(NA,t(var.monthly[var.monthly$month==11,][var.name])[,-y.rm])

#move all months in appropriate colum for TRW correlation comparison, current year
tmp$Ja[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]  <- c(t(var.monthly[var.monthly$month==1,][var.name]))
tmp$`F`[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))] <- c(t(var.monthly[var.monthly$month==2,][var.name]))
tmp$M[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]   <- c(t(var.monthly[var.monthly$month==3,][var.name]))
tmp$Ap[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]  <- c(t(var.monthly[var.monthly$month==4,][var.name]))
tmp$My[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]  <- c(t(var.monthly[var.monthly$month==5,][var.name]))
tmp$J[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]   <- c(t(var.monthly[var.monthly$month==6,][var.name]))
tmp$Jl[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]  <- c(t(var.monthly[var.monthly$month==7,][var.name]))
tmp$A[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]   <- c(t(var.monthly[var.monthly$month==8,][var.name]))
tmp$S[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]   <- c(t(var.monthly[var.monthly$month==9,][var.name]))
tmp$O[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]   <- c(t(var.monthly[var.monthly$month==10,][var.name]))
tmp$N[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]   <- c(t(var.monthly[var.monthly$month==11,][var.name]))
tmp$D[1:length(c(t(var.monthly[var.monthly$month==1,][var.name])))]   <- c(t(var.monthly[var.monthly$month==12,][var.name]))

df.TRW <- cbind(df.TRW,tmp)

cor_df.TRW <- cor(df.TRW,use="pairwise.complete.obs")
# how does Yizhao manage to obtain mean correlations across all years?

# barplot for 1 year
months<- c("ja","f","m","ap","my","j","jl","a","s","o","n","d","Ja","F","M","Ap","My","J","Jl","A","S","O","N","D")#"MJJA","JJA")

#par(mfrow=c(2,1),mar=c(2,2,0.5,0.5))
if(plot==TRUE){
barplot(cor_df.TRW[months,"TRWi"],names.arg = months_plotting,cex.names=0.5,ylim=c(-0.7,0.7))
}
mtext(paste(var.name), side=3,adj=0.9,line=-2)
#barplot(cor_comp_decade_Larch_swp[months,"TRWi"],names.arg = months_plotting,cex.names=0.5,main= "median TRW")
return(cor_df.TRW)

}
  
```

```{r fct_Timestamp}
# leapyears are not considered by HYBRID, so create Date timestamp using this function.
# function creates subdaily timestamps from HYBRID output. And if needed automatically subsets 
#'@param t an hour of the day
#'@param df data frame containing the variable within a subdaily time series with 30 minute time stepping, year column must be provided and called 'jyear' , the default output name for year from HYBRID, or 'year'.
#'@param sub with sub=TRUE lets one chose a particular hour for which the data should be subsetted ( e.g. t= 12, only midday values or predawn (e.g. t=6) values)
#'@param start first year in data frame
#'@param end last year in data frame. 

create_Timestamp <- function(df,t, sub=FALSE,start=1901,end=2016){
  #not elegant but works: make sure the function treats df regardles of whether header has a jyear or year column and where it is located.
  if(match(c('jyear'),colnames(df))>0 | match(c('year'),colnames(df))>0){
    idx<-match(c('jyear'),colnames(df))[1]
    colnames(df)[idx] <- "year"}
days_in_month <- c(31,28,31,30,31,30,31,31,30,31,30,31)
hours_in_day  <- seq(0,23)

get_ndays     <- vector()
days_in_year  <- vector()

days_in_year <-append(days_in_year, seq(1:days_in_month[1]))
for (m in 2:12){
 days_in_year <-append(days_in_year, seq(1:days_in_month[m]))
}

days_in_year <- rep( days_in_year,each=48) # get day number of the month on each 30-minute time step across the entire year
hours_in_day <- rep(rep( hours_in_day,each=2),365) # get hours for each day on each month for one year
mins30_in_day <- rep(rep(seq(0,30,30),times=24),365) # get 30-min time-stepping for each day on each month for one year

for(y in start:end){
  subs       <- df[df$year ==y ,]
  subs$minute<- mins30_in_day
  subs$hour  <- hours_in_day
  subs$day   <- days_in_year
  subs$month <- rep(rep(1:12,times=days_in_month),each=48)
  get_ndays[y-1900]  <-  max(df[df$year ==y ,]$DOY )
  
  if(y==start){
    out<-subs
  }else{
 out<- rbind(out,subs)
  }

}
  
out$Date <- as.Date(paste(out$year,out$month,out$day,sep = "-"))

#sw_vvol.predawn <- out[out$hour == 6 & out$minute==0,]
#sw_vvol.midday  <- out[out$hour == 12 & out$minute==0,]
if(sub==TRUE){
df_sub <- out[out$hour == t & out$minute==0,]
return(df_sub)
}else{
return(out)
}

}

```

```{r fct_flex_cor}

#' @title Correlate variables with TRWi over varying time periods 
#' @description flex_cor correlates the mean climate or other variables over a combination of time period lengths, from 1 day to 365 day length over the entire year.
#' Package dplyr must be loaded manually for this function to work. Must be added to Namespace in the package
#' @param chrono data.frame containing a tree-ring chronologies, e.g. as obtained by chron of package dplR., with added column 'Year'.
#' @param climate data.frame or matrix containing climatic data in daily resolution, all days must contain a value, with variables arranged by column. Compulsory column names are year. Each year value should be repeated in 365 rows.
#' @param clim_vars column number of variable(s) in climate data.frane to be used for correlation analysis.
#' @param sp string that will be used in naming the output file, can be for example the pecies name
#' @param verbose default is TRUE,which outputs files with the correlation and significance matrices. If set to FALSE the function only returns the period of maximum correlation and maximum correlation value. 
#' @return cr List of two data frames. coef contains the Pearson correlation coefficients for each period length and sig contains the significance level for the correlation.
#' @example
#' #years   <- c(1901:2015)
#' #climate <- diag_daily
#' #chrono  <- comp_cent_Larch
#' #clim_vars <- c("temp" ,"rain")
#' #colnames(diag_daily)[sel[v]]
#' #sp<-"Larch"
flex_cor <- function(years,climate,chrono,clim_vars,sp,verbose = TRUE){
  diag_daily = climate
for(v in clim_vars){
 
  predictor_data = matrix(NA,length(years),365)
  for(y in 1:length(years)){
   
    predictor_data[y,] <- diag_daily[diag_daily$year==years[y],][[v]]
  }
  cr <- matrix(NA,365,365)  # matrix to capture correlation coefficients
  sig <- matrix(NA,365,365) # matrix to capture whether values are significant
  TRW_i <- chrono[chrono$Year >= min(years) & chrono$Year <= max(years) ,]
  counter <- 0

  #to collect the maximum correlation value and its period:
  max_val <- 0
  max_int <- c(NaN,NaN)
  
  # compute single day correlations.
  for(start_ind in 1:365){
    end_ind <- start_ind
    stats <- predictor_data[,start_ind]
      cr[366-end_ind,start_ind] <- cor.test(TRW_i$TRWi,stats,method="pearson")$estimate
      sig[366-end_ind,start_ind] <- cor.test(TRW_i$TRWi,stats,method="pearson")$p.value
      
      # collect maximum correlation value, if larger than previous maximum value:

      if((cr[366-end_ind,start_ind] > max_val) & (start_ind > 50)){

        max_val <- cr[366-end_ind,start_ind]
        max_int <- c(start_ind,end_ind)
      }
      counter <- counter + 1
  }

  for(s in 1:364){
    start_ind <- 1
    end_ind <- start_ind + s
    
    running_sums  <- rowSums(predictor_data[,start_ind:end_ind])

    stats <- running_sums # for sums
    #stats <- running_sums/(s+1) # for average
    cr[366-end_ind,start_ind] <- cor.test(TRW_i$TRWi,stats,method="pearson")$estimate #  (Step 3. in Kaczka et al)
    sig[366-end_ind,start_ind] <- cor.test(TRW_i$TRWi,stats,method="pearson")$p.value
    
    # collect maximum correlation value, if larger than previous maximum value:

    if((cr[366-end_ind,start_ind] > max_val) & (start_ind > 50)){

      max_val <- cr[366-end_ind,start_ind]
      max_int <- c(start_ind,end_ind)
    }

    if(s < 364){
      for(step in 1:(364-s)){
        end_ind <- end_ind + 1
        running_sums <- running_sums - predictor_data[,start_ind] + predictor_data[,end_ind]
        start_ind <- start_ind+1

        stats <- running_sums # for sums
        #stats <- running_sums/(s+1) # for average
        cr[366-end_ind,start_ind] <- cor.test(TRW_i$TRWi,stats,method="pearson")$estimate #  (Step 3. in Kaczka et al)
        sig[366-end_ind,start_ind] <- cor.test(TRW_i$TRWi,stats,method="pearson")$p.value
        
        # collect maximum correlation value, if larger than previous maximum value:

        if((cr[366-end_ind,start_ind] > max_val) & (start_ind > 50)){

          max_val <- cr[366-end_ind,start_ind]
          max_int <- c(start_ind,end_ind)
        }
        
        counter <- counter + 1
        # print(counter/(365*365/2))
      }
    }
  }


  if(verbose==TRUE){
    fcorr  = list(cr,sig)
    names(fcorr) <- c("coef","sign")
    save(fcorr,file=paste(v,sp,"flex_cor.RData",sep="_"))
  }else{
    # to be extended here
    max_vals<-list(max_val,max_int)
    names(max_vals) <- c("max_correlation","max_interval")
    #max_vals <- 1 # to return something for now
    return(max_vals)
  }
  
}
}
```

```{r fct_beta_response_curve}
# beta calculator function linear only:
# parameter to be varied: slope
# x = soil moisutre
# swp_crit = critical threshold below which growth happens. ( bet=0)
btc <- function(x,swp_crit){
  # bet <- (crit - x)  /crit 
  bet <- (swp_crit - x)  /swp_crit 
  return(bet)
}

btc_plot <- function(x,swp_crit){
  # bet <- (crit - x)  /crit 
  bet <- (swp_crit - x)  /swp_crit 
  bet[which(bet<=0)] <- 0
  return(bet)
}

# beta-function to plot CLM4.5 soil water critical thresholds:
btcCLM45 <- function(x,swp_crit,swp_o){
  # bet <- (crit - x)  /crit 
  bet<- (swp_crit - x)  /(swp_crit -swp_o)
  bet[which(bet<=0)] <- 0
  bet[which(bet>=1)] <- 1
  return(bet)
}

```

```{r SOX_model}

# Function to compute the plant stomatal conductance from enviromental variables and plant traits using the Stomatal Optimisation based on Xylem hydraulics (SOX)
# See full description of the numerical implementation of SOX used in this code in Eller et al (2018) 

# Function to calculate "a" from "P50" based on Christoffersen et al (2016) equations
# "P50" shapes the vulnerability curve a
# The argument is a single numeric "P50" value 
# The function returns a single numeric "a" value  
P50toa <- function(P50){
  stem_Slope <- 65.15*(-P50)^(-1.25)
  a <- -4*stem_Slope/100*P50
  return(a)
}

# Collatz et al (1991) C3 photosynthesis model following Clark et al (2011)
# This function returns a numeric "A" vector with length equal to the length of the argument "Ci" 
# All arguments are single numeric values, except "Ci" which can be a numeric vector
# See the "Example sox run" section below for the meaning and units of each argument
Photofun <- function(IPAR, Tc, Oa, Ci, vcmax25, Tup, Tlw, alfa, wPAR, beta1, beta2) {
  
  #Vcmax temperature response
  Q10leaf = 2
  Vcmax = vcmax25 * (Q10leaf ^ (0.1 * (Tc - 25))) / ((1 + exp(0.3 * (Tc - Tup))) * (1 + exp(0.3 * (Tlw - Tc))))
  
  #Photocompensation
  Q10rs = 0.57
  photocomp = Oa / (2 * 2600 * (Q10rs ^ (0.1 * (Tc - 25))))
  
  #Kc & Ko 
  Q10Kc = 2.1
  Kc = 30 * (Q10Kc ^ (0.1 * (Tc - 25)))
  Q10Ko = 1.2
  Ko = 30000 * (Q10Ko ^ (0.1 * (Tc - 25)))
  
  #A  
  Ac = Vcmax * ((Ci - photocomp) / (Ci + Kc * (1 + Oa / Ko)))
  Al = alfa * (1 - wPAR) * IPAR * ((Ci - photocomp) / (Ci + 2 * photocomp))
  Ae = rep(0.5 * Vcmax, length(Ac))
  
  #Co-limitation
  x = Re(mapply(function (Al, Ac, beta1) polyroot(c(Al*Ac,-(Al+Ac),beta1)), Al,Ac,beta1))
  Acol1 = do.call(pmin,lapply(1:nrow(x), function(i)x[i,]))
  x2 = Re(mapply(function (Acol1, Ae, beta2) polyroot(c(Acol1*Ae,-(Acol1+Ae),beta2)), Acol1,Ae,beta2))
  Acol2 = do.call(pmin,lapply(1:nrow(x2), function(i)x2[i,]))
  
  #output
  return(Acol2)
  
} 

# SOX arguments:
# SOX plant hydraulic parameters: krc_max, P50, He 
# If "a" is "NULL" the sox function uses the "P50toa" function to compute "a" from the "P50"
# soil: 1 = Assumes Pr = Ps , 2 = Assumes Pr = f(Ps). See figure 1 in Eller et al (2018) for details 
# Additional soil parameters if soil = 2: ksr_max, Ps_max, b which are parameters from Clapp & Hornberger (1978) equation. See Eller et al (2018) for details
# Photo: Dataframe containing the parameters of the photosynthesis model (vcmax25,Tup,Tlw,alfa,wPAR,beta1,beta2) 
# Env: Dataframe containing the environmental driving data (IPAR, Ta, D, Ps, Oa, Ca, Pa)
# Cilen: Resolution of the numerical optimization, higher is slower but more precise
SOX <- function(krc_max, P50, a = NULL, He, soil = 1, ksr_max, Ps_max, b, Photo, Env, Cilen = 1000) {
  
  # Creating output dataframe
  out <- data.frame(timestep = 1:length(Env$Ta), # same units of the enviromental input dataframe
                    krc = rep(NA, length(Env$Ta)), # Xylem hydraulic conductance from root to canopy (mol m-2 s-1 MPa-1)
                    Akcost = rep(NA, length(Env$Ta)), # Photosynthesis multiplied by kcosts (mol m-2 s-1)
                    Aopt = rep(NA, length(Env$Ta)), # Photosynthesis (mol m-2 s-1)
                    Ci = rep(NA, length(Env$Ta)), # Leaf internal CO2 partial pressure (Pa)
                    kcost = rep(NA, length(Env$Ta)), # Normalised xylem hydraulic conductance
                    Gc = rep(NA, length(Env$Ta)), # Stomatal conductance to CO2 (mol m-2 s-1)
                    E = rep(NA, length(Env$Ta)), # Transpiration (mol m-2 s-1)
                    Pc = rep(NA, length(Env$Ta)), # Canopy water potential (MPa)
                    Pr = rep(NA, length(Env$Ta))) # Root water potential (MPa)
  

    # Start loop
    for(i in 1:length(Env$Ta)) {
    
    # calculating a from P50
    if (is.null(a) == TRUE) {
      a = P50toa(P50)
      }

    # Creating Ci vector
    Ci = seq(0, 0.99 * Env$Ca[i], length.out = Cilen)
    
    # Run photo model
    A.vec = Photofun(IPAR = Env$IPAR[i], Tc = Env$Ta[i], Oa = Env$Oa[i], Ci = Ci,
                     vcmax25 = Photo$vcmax25, Tup = Photo$Tup, Tlw = Photo$Tlw, 
                     alfa = Photo$alfa, wPAR = Photo$wPAR, beta1 = Photo$beta1, beta2 = Photo$beta2)

    # Getting initial krc from previous timestep
    if (i == 1) {
      krc_prev = krc_max} else {
        krc_prev = out$krc[i-1] 
      }
      
    # Cost function
    E.vec = A.vec * 1.6 * Env$D[i]/((Env$Ca[i] -  Ci)/Env$Pa[i])
    
    if (soil == 1) {
      
      Pr = Env$Ps[i] 
      Pcpd = Pr - He * 0.01
      Pc.vec = Pcpd - E.vec / krc_prev
      Pc.midvec = (Pcpd + Pc.vec) / 2
      kcost = (1 + (Pc.midvec / P50) ^ a) ^ -1
      
      } else if (soil == 2) {
        
        W = ( Env$Ps[i] / Ps_max ) ^ (1/-b)
        ksr = ksr_max * W ^ (2 * b + 3)
        Pr = Env$Ps[i] - E.vec / ksr
        Pcpd = Pr - He * 0.01
        Pc.vec = Pcpd - E.vec / krc_prev
        Pc.midvec = (Pcpd + Pc.vec) / 2
        kcost = (1 + (Pc.midvec / P50) ^ a) ^ -1
        
        }
    
    # Optimum A
    out$Aopt[i] = A.vec[which.max(A.vec * kcost)]
    
    # Optimum Akcost
    out$Akcost[i] = max(A.vec * kcost, na.rm = T)
    
    # Optimum Ci
    out$Ci[i] = Ci[which.max(A.vec * kcost)]
    
    # Optimum kcost
    out$kcost[i] = kcost[which.max(A.vec * kcost)]
    out$krc[i] = out$kcost[i] * krc_max
    
    # Optimum Gc
    out$Gc[i] = out$Aopt[i]/((Env$Ca[i] -  out$Ci[i])/Env$Pa[i]) 
    
    # Optimum E
    out$E[i] = out$Gc[i] * 1.6 * Env$D[i]
    
    # Optimum Pc
    if (soil == 1) {
      
      out$Pr[i] = Env$Ps[i] 
      Pcpd = out$Pr[i] - He * 0.01
      out$Pc[i] = Pcpd - (out$E[i]/out$krc[i])
      
      } else if (soil == 2)  {
        
        W = (Env$Ps[i] / Ps_max ) ^ (1/-b)
        ksr = ksr_max * W ^ (2 * b + 3)
        out$Pr[i] = Env$Ps[i] - out$E[i] / ksr
        Pcpd = out$Pr[i] - He * 0.01
        out$Pc[i] = Pcpd - (out$E[i]/out$krc[i])
      
        }
    
  }
  
  # output
  return(out)
  
}

# Example sox run
# Enviromental data dataframes
Env1 <- data.frame(IPAR = rep(0.002, 10e1), # Photosynthetic active radiation incident on the leaf (mol m-2 s-1)
                   Ta = rep(20, 10e1), # Air temperature (C)
                   D = seq(0.5/100,5/100, length.out = 10e1), # Leaf to air vapour pressure deficit (mol mol-1)
                   Ps = rep(-0.1, 10e1), # Water potential on the soil close to the roots (MPa)
                   Oa = rep(21000, 10e1), # O2 partial pressure (Pa)
                   Ca = rep(40, 10e1), # CO2 partial pressure (Pa)
                   Pa = rep(101325, 10e1)) # Atmospheric pressure (Pa)

Env2 <- data.frame(IPAR = rep(0.002, 10e1), 
                   Ta = rep(20, 10e1), 
                   D = rep(0.5/100, 10e1), 
                   Ps = seq(-0.1,-5, length.out = 10e1), 
                   Oa = rep(21000, 10e1), 
                   Ca = rep(40, 10e1), 
                   Pa = rep(101325, 10e1)) 

Env3 <- data.frame(IPAR = seq(0.00001,0.002, length.out =10e1),
                   Ta = rep(20, 10e1), 
                   D = rep(0.5/100,10e1), 
                   Ps = rep(-0.1, 10e1), 
                   Oa = rep(21000, 10e1),
                   Ca = rep(40, 10e1), 
                   Pa = rep(101325, 10e1))

Env4 <- data.frame(IPAR = rep(0.002, 10e1),
                   Ta = rep(20, 10e1), 
                   D = rep(0.5/100,10e1),
                   Ps = rep(-0.1, 10e1), 
                   Oa = rep(21000, 10e1),
                   Ca = seq(20,120, length.out = 10e1), 
                   Pa = rep(101325, 10e1)) 

# Plant photosynthetic parameters dataframe
Photo <- data.frame(vcmax25 = 0.0001, # Rubisco maximum carboxylation rate (mol m-2 s-1)
                    alfa = 0.1, # quantum efficicency (mol mol-1)
                    wPAR = 0.15, # leaf scattering coefficient (unitless)
                    Tup = 50, # high temperature photosynthesis range (C)
                    Tlw = 10, # low temperature photosynthesis range (C)
                    beta1=0.83, # co-limitation factor (unitless)
                    beta2=0.93) # co-limitation factor (unitless)

# SOX xylem hydraulic parameters
P50 = -1.0 #-2.5 # Water potential when plant loses 50% of krc_max (MPa)
krc_max = 0.01 # Maximum soil to canopy hydraulic conductance (mol m-2 s-1 MPa-1)
He= 20 # Canopy height (m)



# run sox
#out1 <- SOX(P50 = P50, krc_max = krc_max, He = He, Photo = Photo, Env = Env1)
#out2 <- SOX(P50 = P50, krc_max = krc_max, He = He, Photo = Photo, Env = Env2)
#out3 <- SOX(P50 = P50, krc_max = krc_max, He = He, Photo = Photo, Env = Env3)
#out4 <- SOX(P50 = P50, krc_max = krc_max, He = He, Photo = Photo, Env = Env4)

# plotting
#par(mfrow = c(2,2))
#plot(out1$Gc ~ Env1$D, type = "l", 
#     ylab = expression(g[c]~(mol~m^-2~s^-1)),
#     xlab = expression(D~(mol~mol^-1)))
#plot(out2$Gc ~ Env2$Ps, type = "l",
#     ylab = expression(g[c]~(mol~m^-2~s^-1)),
#     xlab = expression(Psi[s]~(mol~mol^-1)))
#plot(out3$Gc ~ Env3$IPAR,type = "l",
#     ylab = expression(g[c]~(mol~m^-2~s^-1)),
#     xlab = expression(I[par]~(mol~m^-2~s^-1)))
#plot(out4$Gc ~ Env4$Ca,type = "l",
#     ylab = expression(g[c]~(mol~m^-2~s^-1)),
#     xlab = expression(c[a]~(Pa)))
#par(mfrow = c(1,1))

# References:
# Christoffersen et al (2016) Geoscientific Model Development 9: 4227-4255
# Clapp & Hornberger (1978) Water resources research 14: 601-604
# Clark et al (2011) Geoscientific Model Development 4: 701-722
# Collatz et al (1991) Agricultural and Forest meteorology 54: 107-136
# Eller et al (2018) Proceedings of the Royal Society B 


# from Eller et al 2020, to reproduce fig 1.c):
Env2 <- data.frame(IPAR = rep(0.002, 10e1), 
                   Ta = rep(20, 10e1), 
                   D = rep(0.5/100, 10e1), 
                   Ps = seq(-0.1,-5, length.out = 10e1), 
                   Oa = rep(21000, 10e1), 
                   Ca = rep(40, 10e1), 
                   Pa = rep(101325, 10e1)) 
```

```{r CLM4.5_photosynthesis_and_stomatalC}

            CLM45<-function(forc_pbot,Ta,cair,oair,btran){
              # Elements taken from Subroutine CanopyFluxesMod.f90 from CLM4.5
              # git tag:clm4_5_00
              # removed 
              #-- daylength factor treatments
              #-- sun shade differentiatial treatment
              #-- high temperature inhibition + scaling
              #-- any canopy layer-related code
              #-- any nitrogen related code
              #-- removed 10 day temperature acclimation- related code
              # assumption that rubiso-limited Photosynthesis is the most important one and therefore 
              # is the only one applied for first order responses.
            ##########
              # forc_pbot; atmospheric pressure (Pa) 
              # Ta; air temperature (°C) at input, converted to K.
              # cair;  # atmospheric carbon partial pressure (Pa)
              # oair; O2 partial pressure (Pa)
              # btran;  transpiration wetness (or soil moisture beta) factor (0 to 1)
              tfrz = 273.15     # K
              rgas = 8.31432e03 # [J/Kmol]
              Ta_K = Ta + tfrz   # convert to Kelvin
              forc_pbot 
              
              
              # temperature response function
              ft <- function(tl,ha){  exp( ha / (rgas*1.e-3*(tfrz+25.)) * (1. - (tfrz+25.)/tl) )}
              
            #==============================================================================#
            # Photosynthesis and stomatal conductance parameters, from:
            # Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593
            #==============================================================================#
            
              # maximum rate of carboxylation (umol co2/m**2/s)  
              vcmax_z  = 62.5 #Table 8.1 
              
            # Activation energy, from:
            # Bernacchi et al (2001) Plant, Cell and Environment 24:253-259
            # Bernacchi et al (2003) Plant, Cell and Environment 26:1419-1430
            # except TPU from: Harley et al (1992) Plant, Cell and Environment 15:271-282
            
            kcha    = 79430. 
            koha    = 36380. 
            cpha    = 37830. 
            vcmaxha = 72000. 

            lmrha   = 46390. 
            
            kc = (404.9  / 1.e06 ) *  forc_pbot
            ko = (278.4  / 1.e03 ) *  forc_pbot
            sco  = 0.5  * 0.209  / (42.75  / 1.e06 )
            cp = 0.5  * oair  / sco
            
            # removed all temperature responses on kc, ko, sco, co

            #    if (c3flag ) then
            lmr = vcmax_z * 0.015
            
            # removed all canopy layer loops and leaf nitrogen scalings.
            
            # Maintenance respiration
            
            # if (c3flag ) then
            # removed high temperature inhibition
            lmr_z  = lmr * ft( Ta_K, lmrha) #* fth( Ta_K, lmrhd, lmrse, lmrc)
            
            
            # Adjust for soil water
            vcmax_z  = vcmax_z  * btran 
            lmr_z  = lmr_z  * btran 
            
            ###### photosynthesis
            #if (c3flag ) then
            ci_z  = 0.7  * cair 
            
            # C3: Rubisco-limited photosynthesis => Agross ( = 'ag')
            # assuming that most limitation is Rubisco limited, and therefore disregarding
            # other limitations
            ag  = vcmax_z  * max(ci_z-cp , 0. ) / (ci_z+kc *(1. +oair/ko ))
            
            # Net photosynthesis. 
            #Rd= 0.5E-6 # respiration?
            an  = ag - lmr_z # - Rd?
            
            bbbopt  = 10000. #umol/m**2/s
          
            # Soil water stress applied to Ball-Berry parameters
            bbb  = max (bbbopt *btran , 1. )
            
            mpar = 9 #table 8.1 Oleson et al2014
            hs   = 0.7
            gs   = mpar * an/(cair/forc_pbot)* hs + bbb # umol/m**2/s ! eq. 8.1 Oleson et al 2014
            gs   = gs / 1000000   # convert # umol/m**2/s to # mol/m**2/s
            out <- data.frame(Gc=gs, Anet=an, beta_soil = btran)
            return(out)
            }
```

```{r HYBRID9_photosynthesis_and_stomatalC}
########################################################### HYBRID9 Farquhar-type model photosynthesis implementation #################
#keep all stomata-related mechanisms the same 
#@param Wfac non-stomatal soil-water limitation on photosynthesis
#@param beta stomatal soil-water limitation on photosynthesis

#@param LAI_ind leaf area index (plant type specific)
#@param Nfol  amount of nitrogen in foliage (plant type specific)
#@param sla specific leaf area ( plant type specific)
#@param Vcmax_max  (plant type specific)
#@param Jmax_max (plant type specific)
#@param CO2 CO2 concentration of that year # Unit?
#@param t_air daily air temperature #K? T?
#@param VDD vapour pressure deficit

HYBRID9_photosynthesis <- function(Wfac_phot,beta,Nfol,sla,Vcmax_max, Jmax_max,LAI_ind,CO2,tair,Iapar,VDD,crown_area,Pa,Oi){
  beta = beta # stomatal and
  #nfp  renamed to Wfac # non-stomatal soil water response
  tf = 273.16 # Freezing point of water at 1 atm (273.16)                          (K)
  R  = 8.314  # universal gas constant [J/K/kmole] = 8.314
  k = 0.5
  tair = tair + tf
  
  #! from /MODELS/HX/SOURCE/H8_canopy.f90
  n1 <- exp(-((tair - tf - 31.0) / 18.0) ** 2)
  n2 <-exp (26.35 - 65.33e3 / (8.3144 * tair))
  
  #! Canopy nitrogen (kg[N]/m[leaf]^2)
  Ncan = Nfol  / sla  #! Top leaf N.

  #! kg to mmol N.
  Ncan = Ncan * 1.0E6 / 14.0
  
  rt = 8.31451 * tair
  Vcmax = n2  * Vcmax_max  * Wfac_phot  #! μmol[CO2]/m[leaf]^2/s # non-stomatal soil water reponse on photosynthesis directly
  Jmax  = n1  * Jmax_max   * Wfac_phot  #! μmol[CO2]/m[leaf]^2/s
  
 

  Kc = (404.9  / 1.e06 ) *  Pa #! from CLM45 tagCLM4_5_19  [Pa]
  Ko = (278.4  / 1.e03 ) *  Pa #! from CLM45 tagCLM4_5_19 [Pa]
  sco  = 0.5 * 0.209 / (42.75 / 1.e06) #! from CLM45 tagCLM4_5_19
  pcp = 0.5 * Oi / sco  #! from CLM45 tagCLM4_5_19 [Pa], renamed crom cp25 to pcp
  
  #! Based on hybrid8.
  #pcp = exp(19.02 - 37.83e3 / rt ) #! μmol/mol
  #Kc  = exp(38.05 - 79.43e3 / rt ) #! μmol/mol
  #Ko  = exp(20.30 - 36.38e3 / rt ) #! μmol/mol
  #Oi=208
  Km = Kc * (1.0 + Oi / Ko) #! μmol/mol
  
  Ci = beta * 0.7 * CO2  #!  ! μmol/mol
  #print(Ci)
  #! Ds is the humidity deficit at the leaf surface (kPa)
  Ds = VDD * rt / 18.0
  
  #DO I = 1, 10
  
  #! Carboxylation-limited gross photosynthesis (μmol/m[leaf]^2/s)
  Acarb = Vcmax * (Ci - pcp) / (Ci + Km)
  
  Iapar = Iapar * 1000000  # mol/m[leaf]^2/s to µmol/m[leaf]^2/s
  #1 W/m2 ≈ 4.6 μmole.m2/s
  Je = Jmax * Iapar / (Iapar + 2.1 * Jmax)
  
  #! RuBP regeneration-limited gross photosynthesis (μmol/m[leaf]^2/s)
  ARuBP =  (Ci - pcp) / (4.5 * Ci + 10.5 * pcp) *Je
  
  #! Gross photosynthesis (μmol/m[leaf]^2/s)
  tmp <- data.frame(ARuBP=ARuBP,Acarb=Acarb)
  Ag = apply(tmp,1,min)
  
  #! C3 canopy mitochondrial respiration               (μmol/m[leaf]^2/s)
  # not condisered in the calculation of gpp
  rpcrn = Wfac_phot * 0.2E3 * 14.0E-6 * Ncan * exp(18.72 - 46.39E3 / rt)
  
  #! Net photosynthesis (μmol/m[leaf]^2/s)
  An = Ag - rpcrn
  An[which(An<0)]<-0
  # ! Scale to canopy using Eqns. (7a;8) of Sellers et al. (1996).
  #LAI_ind = plant_mass (KP,KI,x,y) * sla / crown_area , done in HYBRID, but kept constant at 2 duringt his simulation
  
  FPAR = 1 -  exp(-k  * LAI_ind)
  
  #! GPP for diagnostics (kg[C]/plant/s)
  # GPP derived from photosynthesic activity only. Maintenance respiration is not considered here.
  gpp_p = crown_area * 0.012E-6 * Ag * FPAR / k
  
  ### Stomatal conductance
  #! g0 is the conductance as A -> 0 (mol/m]leaf]^2/s)
  g0 = 0.01 #! (qs Table 2)
  #! D0 is an empirical coefficient (Figure 8)
  D0 = 1.5
  #! a1 is an empirical coefficient (Figure 8)
  a1 = 4.4
  gsc = g0 + a1 * An / ((CO2 - pcp) * (1. + Ds / D0))
  
  
  out <- data.frame(gpp=gpp_p, Anet=An, Agross=Ag,FPAR=FPAR,Gc=gsc)
  return(out)
  
}

```

```{r fct_maximum_density}
#from alistaire on https://stackoverflow.com/questions/58785930/r-find-maximum-of-density-plot
# here used to find maximum density of growth start using baseR
modes <- function(d){
    i <- which(diff(sign(diff(d$y))) < 0) + 1
    data.frame(x = d$x[i], y = d$y[i])
}
```

```{r fcts_optim_related}

#############
# Select bins based on the 'clustered bin approach'. Rather than 
# breaks are set between clusters of data points (TRWis) which are assumed to be more similar to each other because they have grown under more similar soil water regimes, which better reflects categorisation of biological responses to tree ring increment.
get_bin_bounds <-function(seq,n_bins = 10,min_bin_size = 6){
    n <- length(seq)
    sorted_seq <- sort(seq)
    seq_diffs <- diff(sorted_seq)
    diff_order <- order(seq_diffs,decreasing = TRUE)
    seq_diffs <- seq_diffs[diff_order]
    sorted_seq <- sorted_seq[1:(n-1)][diff_order]
    bin_lower_bounds <- c(max(seq),min(seq)-0.001)
   # plot(seq)
    i <- 0
    while((length(bin_lower_bounds) < n_bins+1) & (i < n)){
      i <- i + 1
      bin_lower_bound <- sorted_seq[i] + seq_diffs[i]/2
      neirest_bot <- max(bin_lower_bounds[bin_lower_bounds < bin_lower_bound])
      if(length(which((seq <=  bin_lower_bound)&(seq > neirest_bot))) < min_bin_size){next}
      
      neirest_top <- min(bin_lower_bounds[bin_lower_bounds >= bin_lower_bound])
      if(length(which((seq > bin_lower_bound)&(seq <= neirest_top))) < min_bin_size){next}
      
      bin_lower_bounds <- c(bin_lower_bounds,bin_lower_bound)
      m <- length(bin_lower_bounds)
    }
    bin_lower_bounds <- sort(bin_lower_bounds)
    return(bin_lower_bounds)
}

###########
# Function creates different versions of Fig.5
# Type = "growth_model" will provide the additive growth model with fixed origin through 0 and evaluate the proportion of variance explained (R2) by the linear model.
  # type = "corr"      will provide simple correlation on TRWi vs TRWs >> r : sample correlation coefficient
  # type = "exp_var"   visualised the variance explained between TRWi and TRWs by the model.

plot_all_update <- function(start, end, w_crit,n_bins,TRWi_sel,swp_sel,type='corr',species,incr){
  
  # for a given w_crit, simulate tree growth for each year in each bin
  output <- simulate_growth_in_bins(start,end,w_crit,n_bins=n_bins,TRWi_sel=TRWi_sel,swp_sel=swp.sel,species=species)
  
  # create mean simulated growth per bin
  x_y   <- simulate_mean_bin_growth(output,sp=species)
  
  sim  <- x_y[[1]]
  obs  <- x_y[[2]]
  w    <- x_y[[3]]
  
  if(species=="Larch"){col=Larch_col}
  if(species=="Spruce"){col=Spruce_col}
  
  if(type=='exp_var'){
    sim = sim*incr
    increasing_order = order(obs)
    plot(obs[increasing_order],ylab="TRW values",xlab="bin_index",
         ylim=c(min(obs,sim)*0.9,max(obs,sim)*1.1),col=col,pch=21,xaxt ="n",yaxt="n")
    axis(2, mgp=c(3, .6, 0))
    axis(1, mgp=c(3, .6, 0))
    points(sim[increasing_order],color="blue")
    segments(x0=1:length(obs),y0=mean(obs),y1=obs,col = "red")
    segments(x0=(1:length(obs))+0.01,y0=obs,y1=sim,col="black")
    lines(c(1,length(obs)),c(mean(obs),mean(obs)),col="red")
    mtext("Bin index", side=1,line=2)
    mtext('Deviation from TRWi mean', side=2, line=2)
  }
  if(type=='corr'){
    sim = sim*incr
    plot( sim ,obs,ylab="",xlab="",ylim=c(min(obs)*0.8,max(obs)*1.1),xlim=c(min(obs)*0.8,max(obs)*1.1),col=col,pch=21,bg=col,xaxt ="n",yaxt="n")
    axis(2, mgp=c(3, .6, 0), tck=-.015, labels=NA)
    axis(side = 2, lwd = 0, line = -0.8, las = 0.5)
    
    #lines(seq(0,1.9,0.01),seq(0,1.9,0.01))
    abline(lm(obs~ sim),lty=2)
    fit <- summary(lm(obs ~ sim))
    r2   <- fit$r.squared
    rho <- cor(obs, sim)
    coeffs <- fit$coefficients ##new
    slope <- coeffs[2,1]
    print(c(rho^2,r2)) # consistency check
    mtext(side= 3, line= -1.2, paste("r =",round(rho,3),sep=" "),adj=0.04)# r='rho' or sample correlation coefficient
  }
  if(type=='growth_model'){
    #sim=sim*incr
    plot( sim ,obs,ylab="",xlab="",xlim=c(0,(max(sim)+15)),ylim=c(0,2),col=col,pch=21,bg=col,xaxt ="n",yaxt="n")
    axis(2, mgp=c(3, .6, 0), tck=-.015, labels=NA)
    axis(side = 2, lwd = 0, line = -0.8, las = 0.5)
 
    abline(lm(obs~ sim-1,weights = w),lty=2)
    fit <- summary(lm(obs ~ sim-1,weights = w)) # whether lm() fit or our  a_optim funtion, the slope is the same ( pretty much) 
    rho   <- fit$r.squared
    obs_mean_weighted <- sum(w*obs)/sum(w)
    r2 <- 1-sum(w*(sim*incr-obs)^2)/sum(w*(obs-obs_mean_weighted)^2) # R2 (distance from additive growth model) must be calculated taking into account weighted samples and centering. 
    coeffs <- fit$coefficients ##new
    slope <- coeffs[1,1]
    #print(c(rho^2,r2)) # consistency check
    #https://stats.stackexchange.com/questions/134167/is-there-any-difference-between-r2-and-r2
    numb <- round(max(r2),3)
    mtext(side= 3, line= -1.6,  bquote(R^2~' = '~.(numb)),adj=0.03) # R= R2 , proportion of variance explained
 
    mtext(side= 3, line= -2.6,  paste("   y = ",(round(slope-0.0001,3)),'x'),adj=0.03 ) #,round(slope[1,1],3))
    # PLEASE NOTE: somehow round(slope_optim_L,3) != round(slope,3), but when printing out, they are both 0.0245.
    # since this is the case, I will be removing 0.0001 for plotting, but that makes the code less flexible unfortunately and needs to be considered when running for other conditions.

    
  }
  
  if(species=="Spruce"){
    axis(1, mgp=c(3, .6, 0), tck=-.015, labels = NA)
    axis(side = 1, lwd=0,line=-1,las=1)
  }
  samples_per_bin <- output[[2]]
  return(samples_per_bin$sample_size)
}

###########
  # function works over all number of bins and finds the years and the soil moisture content within each of these years.
  # Within each bin, for each year it then generates separate 'beta response factor' time series based on the respective soil water time series.
  # This is then passed out of the function, together with the binned TRWis.
simulate_growth_in_bins<-function(start,end, w_crit,n_bins,TRWi_sel,swp_sel,species){
  
  # Binning of the TRWs and calculating of the median swp (layer numbers are flexible)
  TRW_bin_bounds <- get_bin_bounds(TRWi_sel$TRWi,n_bins=n_bins) #quantile(TRWi_sel$TRWi,seq(0,1,binsize))
  TRW_bins       <- n_bins # number of bins
  swp_bins       <- list()
  
  #TRWi_sim will contain binned observed TRWi data and modelled (simulated) TRWi data. 
  TRWi_sim           <- data.frame(TRW_i=rep(NA,TRW_bins), bin=rep(NA,TRW_bins))
  swp.pred_TRWbins <- list()
  swp.mdd_TRWbins  <- list()
  interval_data    <- list()
  output <- list()
  record_TRW_in_bins <- vector()
  #important for later, must be initialised with the above TRWi_sim
  significant_period_betas <- matrix(0,ncol=2,nrow=dim(TRWi_sim)[1]) # col 1 record beta sum, col2 record percentage of days set to 0 because of response function
  
  sorted_TRWi <- TRWi_sel[order(TRWi_sel$TRWi),]
  
  for (i in 1:(length(TRW_bin_bounds)-1)){
    lb <- TRW_bin_bounds[i]
    ub <- TRW_bin_bounds[i+1]
    bin_indices <- which((sorted_TRWi > lb) & (sorted_TRWi <= ub))
    # obtain the TRWs in each bin
    tmp <- sorted_TRWi[bin_indices,]
    # record the number of trw in each bin, to see whether the binning skewes the samples too much.
    record_TRW_in_bins[i] <- dim(tmp)[1]
    
    # select years with TRWs
    ys                        <- sort(tmp$Year)
    # collect binned mean observed TRWi data.
    TRWi_sim$TRW_i[i]       <- mean(tmp$TRWi)
    TRWi_sim$bin[i]         <- i
    TRWi_sim$sample_size[i] <- length(tmp$TRWi)
    # obtain the swps from these TRWs ( via years)
    swp.mdd_fromTRWbins   <- swp_sel[year(swp_sel$Date) %in% ys ,]
    
    predictor_data <- matrix(rep(NaN,(end-start+1)*length(ys)),nrow = length(ys))
   
     for(y in 1:length(ys)){
      subs <- swp.mdd_fromTRWbins[year(swp.mdd_fromTRWbins$Date) == ys[y],]
      subs <- subs[subs$DOY>= start & subs$DOY<=end,]
      
      if(species=="Larch"){
      predictor_data[y,] <- subs[["swp4"]]
      }
      if(species=="Spruce"){
      predictor_data[y,] <- subs[["swp3"]]  # spruce has more shallow roots compared to Larch.
      }
    }
    swp.mdd_TRWbins[[i]] <- predictor_data
  }
  
  # collect betas over the significant period 
  # collects betas for selected years in each bin separately
  for( i in 1:(length(TRW_bin_bounds)-1)){
    subs     <- swp.mdd_TRWbins[[i]]
    beta_rw  <- btc(subs,w_crit)
    
    significant_period_betas[i,2]      <- length(beta_rw[beta_rw < 0 ])/dim(beta_rw)[1] 
    beta_rw[beta_rw < 0] <- 0 # filter out values the function cannot handle. as done in HYBRID9
    beta_rw[beta_rw > 1] <- 1
    # sum up for entire period, for selected years in each bin:
    output[[i]] <- rowSums(beta_rw)
  }
  return(list(output,TRWi_sim))
}

########### 
  # function calculates growth for each year
simulate_growth <-function(start,end, w_crit,n_bins,TRWi_sel,swp_sel,incr,species){
 
 
  TRW_samples     <- dim(TRWi_sel)[1] 
  swp_bins       <- list()
  
  #TRWi_sim will contain binned observed TRWi data and modelled (simulated) TRWi data. 
  TRWi_sim           <- data.frame(TRW_i=rep(NA,TRW_samples ), bin=rep(NA,TRW_samples ))
  swp.pred_TRWbins <- list()
  swp.mdd_TRWbins  <- list()
  interval_data    <- list()
  output <- list()
  record_TRW_in_bins <- vector()
  #important for later, must be initialised with the above TRWi_sim
  significant_period_betas <- matrix(0,ncol=2,nrow=dim(TRWi_sim)[1]) # col 1 record beta sum, col2 record percentage of days set to 0 because of response function
  

  ys <- TRWi_sel$Year
  predictor_data <- matrix(rep(NaN,(end-start+1)*length(ys)),nrow = length(ys))
  
  
  swp.mdd_fromTRWbins   <- swp_sel[year(swp_sel$Date) %in% ys ,]
  
    for(y in 1:length(ys)){
      subs <- swp.mdd_fromTRWbins[year(swp.mdd_fromTRWbins$Date) == ys[y],]
      subs <- subs[subs$DOY>= start & subs$DOY<=end,]
      if(species=="Larch"){
      predictor_data[y,] <- subs[["swp4"]]
      }
      if(species=="Spruce"){
        predictor_data[y,] <- subs[["swp3"]]
      }
    }
    swp.mdd_TRWbins <- predictor_data
    #swp.mdd_fromTRWbins[which(swp.mdd_fromTRWbins$swp4==-0.013926),]# test whether correct days are picked
    
    TRWi_sel$TRW_sim <-NA
  # collect betas over the significant period for each year and sum them up as simulated TRW
  for( i in 1:length(ys)){
    subs     <- swp.mdd_TRWbins[i,]
    beta_rw  <- btc(subs,w_crit)
    
    significant_period_betas[i,2]      <- length(beta_rw[beta_rw < 0 ])/length(beta_rw)
    beta_rw[beta_rw < 0] <- 0 # filter out values the function cannot handle. as done in HYBRID9
    beta_rw[beta_rw > 1] <- 1
    # sum up for entire period, is now one value per year:
    # store for each year
    #output[[i]] <- rowSums(beta_rw)
    TRWi_sel$TRW_sim[i]     <- sum(beta_rw)
  }

# Scale using the provided value of characteristic daily increment
    TRWi_sel$TRW_sim_norm <-TRWi_sel$TRW_sim *incr # + intercept? But for us, we want intercept to be 0
return(TRWi_sel)
}

###########
# simple step-wise optimisation, which always rejects obj, when it is larger than before. We only have one global maximum, so this is an appropriate approach.
scaling_factor_optim <- function(sim,obs,w,a_min=0.001,a_max=0.1,step = 0.0001){
 
  optim_a <- NaN
  optim_obj <- Inf
  n_steps <- (a_max-a_min)/step
  for(i in 0:n_steps){
    a <-a_min + i*step
    obj <- sum(w*(a*sim-obs)^2) # sum of squared residual errors (SSRE or SRSE)
    if(obj <= optim_obj){
      optim_a <- a
      optim_obj <- obj
    }else if(obj > 1.1*optim_obj){
      break }
  }
  return(list(optim_a,optim_obj))
}

###########
# function receives output from simulate_growth_in_bins 
  #    sim_output[[1]] = simulated growth for selected years within each bin
  #    sim_output[[2]] = binned TRWis from observations. Bin identifier and the number of TRWi - values in each bin.
  
  # computes the mean simulated growth (per bin)  from simulated growth for each year within a bin
  # compares mean simulated growth (per bin) to mean observed TRWis (per bin) 
simulate_mean_bin_growth <- function(sim_output,sp){
  
  simulated_bin_growth <- sim_output[[1]]
  TRWi_sim <- sim_output[[2]]
  
  significant_period_betas <- matrix(0,ncol=2,nrow=dim(TRWi_sim)[1])
  
  for(i in 1:length(simulated_bin_growth)){
    significant_period_betas[i,1]     <- mean(simulated_bin_growth[[i]])
  }

  TRWi_sim$TRW_sim  <- significant_period_betas[,1]
  
  y <- TRWi_sim$TRW_i     # observed
  x <- TRWi_sim$TRW_sim   # simulated
  w <- TRWi_sim$sample_size # bin weights proportional to bin sample size
  
  return(list(x,y,w))
}


###########
  # compares simulated and observed TRWis using a linear regression through the origin.
  # output are R2 and regression slope.
  # slope is important for computing MEF later, as it is used in the scaling.
  # w is used to weigh the bins according to their number of samples ( as clustered binning results in uneven number of samples per bin)

compare_growth <- function(x_y){
 
  x <- x_y[[1]]
  y <- x_y[[2]]
  w <- x_y[[3]]
  # useful: https://stats.stackexchange.com/questions/26176/removal-of-statistically-significant-intercept-term-increases-r2-in-linear-mo
  # https://stats.stackexchange.com/questions/134167/is-there-any-difference-between-r2-and-r2
  optim_solution <- scaling_factor_optim(x,y,w)
  a_optim <- optim_solution[[1]]
  srse <- optim_solution[[2]]
  weighted_mean <- sum(w*y)/sum(w)
  r2 <- 1-srse/sum(w*(y-weighted_mean)^2)
  return(list(r2,a_optim))
}

compute_r_squared<-function(y_pred,y,intercept=TRUE){
  # formula from ?summary.lm
  if(intercept){
    mean_y <- mean(y)
  }else{
    mean_y <- 0.
  }
  res <- y_pred - y
  r.squared <- 1 - (sum(abs(res)^2)  / sum( abs(y-mean(y))^2))
  return(r.squared)
}

MEF_optim <- function(start, end, w_crit,n_bins,TRWi_sel,swp_sel,incr,sp){
  
  # for a given w_crit, simulate tree growth for each year in each bin
  output <- simulate_growth_in_bins(start,end,w_crit,n_bins=n_bins,TRWi_sel=TRWi_sel,swp_sel=swp.sel,sp=sp)
  
  # create mean simulated growth per bin
  x_y   <- simulate_mean_bin_growth(output,sp=sp)
  
  sim <- x_y[[1]]
  obs  <- x_y[[2]]
  w <- x_y[[3]]
  
  MEF = 1 - (sum(w*(( incr * sim )- obs )^2)  / sum(w*(obs - mean(obs))^2 ))
  
  return(c(MEF))
  
}

```


#Downscaling

```{r plot_downscaling_decade, eval=FALSE}
jpeg(filename='figures/APPDX_Downscaling_decade.jpeg')
par(mfrow=c(2,1),mar=c(1,4,1,0.5),oma=c(3,0,0,0))

# CRUNCEP gridded data product for grid containing our study site
HYBRID_CRUNCEP <- read.table('data/CRUNCEP_MET.dat', head=TRUE)
HYBRID_CRUNCEP.xts   <- aggregate_to_daily_met_variables(HYBRID_CRUNCEP,"1901-01-01")
met_loc_CRUNCEP      <- data.frame(date=index(HYBRID_CRUNCEP.xts), coredata(HYBRID_CRUNCEP.xts))  

# locally observed temperature and rainfall close to our study site
met_loc_obs         <- read.table('data/MIX_MET_for_HYBRID9_N08.txt',head=TRUE) # rain and tair are from local observations

##################################  Bias Correction or 'spatial downscaling'  ##################################

################# TEMPERATURE
t_2011_obs     <- met_loc_obs[met_loc_obs$year>=2011,]$tair
t_2011_CRUNCEP <- met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2011-01-01",]$tair
t_ref_obs      <- met_loc_obs[met_loc_obs$year<2011,]$tair                                                       # reference period for observations
t_ref_CRUNCEP  <- met_loc_CRUNCEP[met_loc_CRUNCEP$date<"2011-01-01" & met_loc_CRUNCEP$date>="2006-01-01",]$tair  # reference period for CRUNCEP
##


# perform temperature bias correction/downscaling period 
# divide into two periods, above and = or below 0. obtain sd for both these periods
# use only the above 0 for the bias correction:

sigma_t_obs_a0 <- sd(t_ref_obs[t_ref_obs > 273.15])
sigma_t_obs_b0 <- sd(t_ref_obs[t_ref_obs <= 273.15])

sigma_t_CRUNCEP_a0  <- sd(t_ref_CRUNCEP[t_ref_obs>273.15])
sigma_t_CRUNCEP_b0  <- sd(t_ref_CRUNCEP[t_ref_obs<=273.15])

tair_CRUNCEP_update                         <- rep(0,length(t_2011_CRUNCEP))
tair_CRUNCEP_update[t_2011_CRUNCEP>273.15]  <- mean(t_ref_obs) + (sigma_t_obs_a0/sigma_t_CRUNCEP_a0) * (t_2011_CRUNCEP[t_2011_CRUNCEP>273.15] - mean(t_ref_CRUNCEP))
tair_CRUNCEP_update[t_2011_CRUNCEP<=273.15] <- mean(t_ref_obs) + (sigma_t_obs_b0/sigma_t_CRUNCEP_b0) * (t_2011_CRUNCEP[t_2011_CRUNCEP<=273.15] - mean(t_ref_CRUNCEP))

met_loc_CRUNCEP$tair_obs                                          <- NA
met_loc_CRUNCEP$tair_update                                       <- NA
met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2006-01-01",]$tair_obs     <- met_loc_obs$tair
met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2011-01-01",]$tair_update  <- tair_CRUNCEP_update  # corrected temperature, for all years above 2011


# plot temp dosnwcaled, CRUNCEP and observed; for comparison
plotprep <- met_loc_CRUNCEP[met_loc_CRUNCEP$date>= '2006-01-01',]
plotprep$tair_update <- NA
plotprep$tair_obs    <- met_loc_obs$tair
plotprep[plotprep$date>="2011-01-01",]$tair_update<-tair_CRUNCEP_update
plot_met <- plotprep

plot(plot_met$date,plot_met$tair-273.15,type="l",ylim=c(-25,30),col=CRUNCEP_col,ylab='Air Temp [°C]',xlab='Year')
lines(plot_met$date,plot_met$tair_obs-273.15,col=obs_colour)
lines(plot_met$date,plot_met$tair_update-273.15,col="black")


########## PRECIPITATION 
met_loc_CRUNCEP$month <- lubridate::month(met_loc_CRUNCEP$date)
met_loc_CRUNCEP$year  <- lubridate::year(met_loc_CRUNCEP$date)
met_loc_CRUNCEP$rain_update <-NA

met_loc_obs$date <- as.Date(met_loc_obs$day-1, origin = paste(met_loc_obs$year,"-01-01",sep=""))
met_loc_obs$month <- month(met_loc_obs$date)

# perform precip downscaling on a monthly basis
 for( m in 1:12){
  modprecip       <- met_loc_CRUNCEP[which(met_loc_CRUNCEP$date >="2011-01-01" & met_loc_CRUNCEP$month == m),]$rain
   modprecip.ref  <- met_loc_CRUNCEP[which(met_loc_CRUNCEP$date >="2006-01-01"  & met_loc_CRUNCEP$month == m),]$rain
    obsprecip     <- met_loc_obs[which(met_loc_obs$date >="2006-01-01"  & met_loc_obs$month ==m), ]$rain
    fit.mon       <- fitQmap(obsprecip,modprecip.ref,
               method="QUANT",qstep=0.001)
      qm1         <- doQmap(modprecip.ref,fit.mon,type="tricub")
    met_loc_CRUNCEP[which(met_loc_CRUNCEP$date >="2006-01-01" & met_loc_CRUNCEP$month == m),]$rain_update <- qm1
 }
 
# plot precip downscaling result for comparison:
plot_biased          <- met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2006-01-01",]
plot_biased$rain_obs <- met_loc_obs$rain

plot(plot_biased$date,plot_biased$rain ,col="orange",ylab= "precipitation [mm/day]",type="h")
mtext("Year", side=1,line=2)
lines(plot_biased$date,plot_biased$rain_obs ,col = obs_colour,type="h")
lines(plot_biased[plot_biased$date>="2011-01-01",]$date,plot_biased[plot_biased$date>="2011-01-01",]$rain_update,type="h")

legend("bottomright",horiz=TRUE,legend=c("CRUNCEP_d","CRUNCEP","Obs"), fill=c(ds_col,CRUNCEP_col, obs_colour),cex=0.7,inset=c(0.2,-0.25),xpd=NA,bty="n")
dev.off()
```

```{r plot_downscaling_century}
jpeg(filename='figures/APPDX_Downscaling_century.jpeg')
par(mfrow=c(2,1),mar=c(1,4,1,0.5),oma=c(4,0,0,0.5))

# locally observed temperature and rainfall close to our study site
met_loc_obs         <- read.table('data/MIX_MET_for_HYBRID9_N08_2006_2016.txt',head=TRUE) # rain and tair are from local observations

# CRUNCEP gridded data product for grid containing our study site
HYBRID_CRUNCEP <- read.table('data/CRUNCEP_MET.dat', head=TRUE)
HYBRID_CRUNCEP.xts   <- aggregate_to_daily_met_variables(HYBRID_CRUNCEP,"1901-01-01")
met_loc_CRUNCEP      <- data.frame(date=index(HYBRID_CRUNCEP.xts), coredata(HYBRID_CRUNCEP.xts))  

################# Temperature bias correction:
# expand reference period across the entire period for which data exists
t_ref_obs       <- met_loc_obs$tair
t_ref_CRUNCEP   <- met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2006-01-01",]$tair

# divide into two periods, above and = or below 0. obtain sd for both these periods
# use only the above 0 for the bias correction:
sigma_t_obs_a0 <- sd(t_ref_obs[t_ref_obs>273.15])
sigma_t_obs_b0 <- sd(t_ref_obs[t_ref_obs<=273.15])
sigma_t_CRUNCEP_a0  <- sd(t_ref_CRUNCEP[t_ref_obs>273.15])
sigma_t_CRUNCEP_b0  <- sd(t_ref_CRUNCEP[t_ref_obs<=273.15])

#create dataframe for 100+ years that holds both a0 and b0 values:
tair_CRUNCEP_update                               <- rep(0,dim(met_loc_CRUNCEP)[1])
tair_CRUNCEP_update[met_loc_CRUNCEP$tair>273.15]  <- mean(t_ref_obs) +  (sigma_t_obs_a0/sigma_t_CRUNCEP_a0) * (met_loc_CRUNCEP[met_loc_CRUNCEP$tair>273.15,]$tair - mean(t_ref_CRUNCEP))
tair_CRUNCEP_update[met_loc_CRUNCEP$tair<=273.15] <- mean(t_ref_obs) + (sigma_t_obs_b0/sigma_t_CRUNCEP_b0) * (met_loc_CRUNCEP[met_loc_CRUNCEP$tair<=273.15,]$tair - mean(t_ref_CRUNCEP))

  
# plot temp dosnwcaled, CRUNCEP and observed; for comparison 
met_loc_CRUNCEP$tair_update                                   <- tair_CRUNCEP_update  # corrected temperature, for entire century
met_loc_CRUNCEP$tair_obs                                      <- NA
met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2006-01-01",]$tair_obs <- t_ref_obs

plot(met_loc_CRUNCEP$date,met_loc_CRUNCEP$tair-273.15,type="l",ylim=c(-25,30),col = CRUNCEP_col,ylab='Air Temp [°C]',xlab='Year')
# add Bias- corrected dataset
lines(met_loc_CRUNCEP$date,met_loc_CRUNCEP$tair_update-273.15, col= ds_col)
# add observations from 2006-2016
lines(met_loc_CRUNCEP$date,met_loc_CRUNCEP$tair_obs-273.15, col = obs_colour)


################ Precip interpolation

# locally observed temperature and rainfall close to our study site
met_loc_obs         <- read.table('data/MIX_MET_for_HYBRID9_N08_2006_2016.txt',head=TRUE) # rain and tair are from local observations


##########
met_loc_CRUNCEP$month <- lubridate::month(met_loc_CRUNCEP$date)
met_loc_CRUNCEP$year <- lubridate::year(met_loc_CRUNCEP$date)

met_loc_obs$date <-  as.Date(met_loc_obs$day-1, origin = paste(met_loc_obs$year,"-01-01",sep=""))
met_loc_obs$month <- month(met_loc_obs$date)

########### DOWNSCALE  entire century
collect_mon_wet.day_thresh<-rep(NA,12)
 met_loc_CRUNCEP$rain_update <- NA
for( m in 1:12){
 # modprecip     <- met_loc_CRUNCEP[ met_loc_CRUNCEP$month == m,]$rain
  modprecip <- met_loc_CRUNCEP[which(met_loc_CRUNCEP$date < "2006-01-01" & met_loc_CRUNCEP$month == m),]$rain
  obsprecip     <- met_loc_obs[met_loc_obs$month ==m, ]$rain
    fit.mon     <- fitQmap(obsprecip,modprecip,
               method="QUANT",qstep=0.001)
    
      qm1       <- doQmap(modprecip,fit.mon,type="tricub")
      
   met_loc_CRUNCEP[ which(met_loc_CRUNCEP$date < "2006-01-01" & met_loc_CRUNCEP$month == m),]$rain_update <- qm1
}
 
# the last century does not fit very well with the transformation, because of wet bias starting to encroach from west alps ( mostly the  Zurich records) becoming wetter in that decade.. Don't use in quantile mapping as it will distort the past frequency distribution. Sion weather station precipiration distribution has been rather stable in precipitation and we carefully infer that Visp therefore has as well. Must not downscale this bias introduced by later 'zurich-year's in the CRUNCEP data.
# Stations N47.5E7.5
#All CRU TS stations used in construction of the gridcell centred at 46.25N, 7.75E
#pre.1905011321.clean.dtb), 
#Will perform independend transformation on it.
 ### for reference period
# this will be replaced width local observations for the simulations anyways, but will be calculated here for plotting:

 for( m in 1:12){
 # modprecip     <- met_loc_CRUNCEP[ met_loc_CRUNCEP$month == m,]$rain
  modprecip.ref <- met_loc_CRUNCEP[which(met_loc_CRUNCEP$date >="2006-01-01" & met_loc_CRUNCEP$month == m),]$rain
  obsprecip     <- met_loc_obs[met_loc_obs$month ==m, ]$rain
    fit.mon     <- fitQmap(obsprecip,modprecip.ref,
               method="QUANT",qstep=0.001)
    collect_mon_wet.day_thresh[m]<-fit.mon$wet.day # for rainfall cococcurence on days
      qm1       <- doQmap(modprecip.ref,fit.mon,type="tricub")
      
   met_loc_CRUNCEP[which(met_loc_CRUNCEP$date >="2006-01-01" & met_loc_CRUNCEP$month == m),]$rain_update <- qm1
 }
 
####################################

met_loc_CRUNCEP$rain_obs <- NA
met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2006-01-01",]$rain_obs<-met_loc_obs$rain

plot(met_loc_CRUNCEP$date,met_loc_CRUNCEP$rain,type="l",col= CRUNCEP_col,ylab='Precip [mm/day]')
lines( met_loc_CRUNCEP$date,met_loc_CRUNCEP$rain_update,type="h")
lines(met_loc_CRUNCEP[met_loc_CRUNCEP$date >= "2006-01-01",]$date,met_loc_CRUNCEP[met_loc_CRUNCEP$date >= "2006-01-01",]$rain_obs, col = obs_colour)
mtext('Year',side=1,line=2)


legend("bottomright",horiz=TRUE,legend=c("CRUNCEP_d","CRUNCEP","Obs"), fill=c(ds_col,CRUNCEP_col, obs_colour),cex=0.7,inset=c(0.2,-0.25),xpd=NA,bty="n")
dev.off()

``` 

```{r met_data_aggregated, eval= FALSE}
jpeg(filename='figures/APPDX_Downscaling_aggrg_century.jpeg')
# Not used in Paper
par(mfrow=c(2,1),mar=c(1,4,1,0.5),oma=c(4,0,0,0.5))

# plot yearly mean temperatures
met_tmp <- apply.yearly(xts(met_loc_CRUNCEP,met_loc_CRUNCEP$date),FUN = mean)
plot(year(met_tmp$date),met_tmp$tair-273.15, type="l",ylim=c(-1,12),col = CRUNCEP_col, ,xlab="Year", ylab="mean annual temperature (°C)")
lines(year(met_tmp$date),met_tmp$tair_update-273.15,col=ds_col)
lines(year(met_tmp$date),met_tmp$tair_obs-273.15,col=obs_colour)

# used in text in manuscript, section spatial downscaling 
mean(met_tmp[year(met_tmp$date) >= 2006,]$tair_obs)-273.15
mean(met_tmp[year(met_tmp$date) >= 2006,]$tair)-273.15
#####################################################################

tmp.xts <- xts(met_loc_CRUNCEP,met_loc_CRUNCEP$date)
met_tmp$rain_obs     <- tapply(met_loc_CRUNCEP$rain_obs,year(met_loc_CRUNCEP$date), sum)
met_tmp$rain_update  <- tapply(met_loc_CRUNCEP$rain_update,year(met_loc_CRUNCEP$date), sum)
met_tmp$rain         <- tapply(met_loc_CRUNCEP$rain,year(met_loc_CRUNCEP$date), sum)

# plot yearly cumulative rainfall
plot(year(met_tmp$date),met_tmp$rain, type="h",ylim=c(-1,3000),col = CRUNCEP_col,xlab="Year", ylab= "cumulative ann. precipitation (mm/year)")
lines(year(met_tmp$date),met_tmp$rain_update,col=ds_col,type="h")
lines(year(met_tmp$date),met_tmp$rain_obs,col=obs_colour,type="h")

# used in text in manuscript, section spatial downscaling
mean(met_tmp[year(met_tmp$date) >= 2006,]$rain_obs)
mean(met_tmp[year(met_tmp$date) >= 2006,]$rain)
dev.off
```

```{r NORUN_create_downscaled_input_file, eval=FALSE}
# create temperature- biase corrected(updated) HYBRID9 inputfile:

# replace tair and rain with downscaled values for entire century
tmp <- as.data.frame(matrix(NA,nrow=length(met_loc_CRUNCEP$tair),ncol=8))
colnames(tmp) <- c('year','day','tair','lwdown','swdown','qair','press','rain')
tmp$year   <- year(met_loc_CRUNCEP$date)
tmp$day    <- as.numeric(strftime(met_loc_CRUNCEP$date,format="%j"))
tmp$tair   <- met_loc_CRUNCEP$tair_update 
tmp$lwdown <- met_loc_CRUNCEP$lwdown
tmp$swdown <- met_loc_CRUNCEP$swdown
tmp$qair   <- met_loc_CRUNCEP$qair
tmp$press  <- met_loc_CRUNCEP$press
tmp$rain   <- met_loc_CRUNCEP$rain_update 

HYBRID_CRUNCEP <- tmp

# write out new forcing data with downscaled rainfall and air temperature:
HYBRID_CRUNCEP<-HYBRID_CRUNCEP[HYBRID_CRUNCEP$year>=1901,]# previous years are just replicated from that century for HYBRID spin-up
write.table(HYBRID_CRUNCEP,'data/MIX_bias_red_prc_t_100plus_met_for_HYBRID_N08.txt',row.names=FALSE)

####################################################

# use observed rainfall and temeprature where possible
# replace the last 11 years with actually observed rainfall data:
met_loc_CRUNCEP[met_loc_CRUNCEP$date >= "2006-01-01",]$rain_update <-met_loc_obs$rain

# replace the last 11 years with actually observed temperature data:
met_loc_CRUNCEP[met_loc_CRUNCEP$date >= "2006-01-01",]$tair_update <- met_loc_obs$tair

# create temperature- biase corrected(updated) inputfile:
tmp <- as.data.frame(matrix(NA,nrow=length(met_loc_CRUNCEP$tair),ncol=8))
colnames(tmp) <- c('year','day','tair','lwdown','swdown','qair','press','rain')
tmp$year   <- year(met_loc_CRUNCEP$date)
tmp$day    <- as.numeric(strftime(met_loc_CRUNCEP$date,format="%j"))
tmp$tair   <- met_loc_CRUNCEP$tair_update # contains downscaled precipitation up to 2006, then local observations from VISP from 2006 to 2016.
tmp$lwdown <- met_loc_CRUNCEP$lwdown
tmp$swdown <- met_loc_CRUNCEP$swdown
tmp$qair   <- met_loc_CRUNCEP$qair
tmp$press  <- met_loc_CRUNCEP$press
tmp$rain   <- met_loc_CRUNCEP$rain_update # contains downscaled precipitation up to 2006, then local observations from VISP from 2006 to 2016.

rm(HYBRID_CRUNCEP)
HYBRID_CRUNCEP<-tmp

# write out new forcing data with downscaled and where possible observed rainfall and air temperature:
HYBRID_CRUNCEP <- HYBRID_CRUNCEP[HYBRID_CRUNCEP$year>=1901,] # pre 1901 years are just replicated from that century for HYBRID9 spin-up
write.table(HYBRID_CRUNCEP,'data/MIX_bias_red_prc_t_met_for_HYBRID_N08.txt',row.names=FALSE)

rm(HYBRID_CRUNCEP.yd)
```

```{r table_downscaling}
met_tmp <- apply.yearly(xts(met_loc_CRUNCEP,met_loc_CRUNCEP$date),FUN = mean)

table_downscaling<-data.frame(CRUNCEP     = sum(met_loc_CRUNCEP[met_loc_CRUNCEP$date >= "2006-01-01",]$rain)/(11), 
                              CRUNCEP_d   = sum(met_loc_CRUNCEP$rain_update)/(116),
                              Observed    = sum(met_loc_obs$rain)/(11),
                              CRUNCEP_dec = sum(met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2006-01-01",]$rain_update)/11 )

table_downscaling  <- data.frame( Variable   = c("Temperature [°C]","Rainfall [mm/year]"),
                                  CRUNCEP    = c(mean(met_tmp[year(met_tmp$date) >= 2006,]$tair)-273.15       ,mean(met_tmp[year(met_tmp$date) >= 2006,]$rain) ),
                                  CRUNCEP_d  = c(mean(met_tmp[year(met_tmp$date) >= 2006,]$tair_update)-273.15,mean(met_tmp[year(met_tmp$date) >= 2006,]$rain_update)),
                                  Observed   = c(mean(met_tmp[year(met_tmp$date) >= 2006,]$tair_obs)-273.15   ,mean(met_tmp[year(met_tmp$date) >= 2006,]$rain_obs)))

knitr::kable(table_downscaling) #\ref{tbl:Downscaling}
```

```{r plot_Appendix_Monthly_Downscaling_qqplot_ecumprop_QMAP.jpeg}

sqrtquant <- function(x,qstep=0.001){
  qq <- quantile(x,prob=seq(0,1,by=qstep))
}
##########################
### DOWNSCALE last decade
for(i in 1:2){ # split into Jan-June and July -Dec for plotting
  
  if( i==1){
jpeg(filename='data/APPDX_Monthly_Downscaling_qqplot_ecumprop_QMAP_1.jpeg')
par(mfrow=c(6,3),mar=c(1.5,1,1.8,3),oma=c(3,3,0.5,0.5))
start=1
end=6}
  if(i==2){
jpeg(filename='data/APPDX_Monthly_Downscaling_qqplot_ecumprop_QMAP_2.jpeg')
par(mfrow=c(6,3),mar=c(1.5,1,1.8,3),oma=c(3,3,0.5,0.5))
start=7
end=12}
for( m in start:end){

  modprecip     <- met_loc_CRUNCEP[ met_loc_CRUNCEP$month == m,]$rain
  modprecip.ref <- met_loc_CRUNCEP[which(met_loc_CRUNCEP$date >="2006-01-01" & met_loc_CRUNCEP$month == m),]$rain
  obsprecip     <- met_loc_obs[met_loc_obs$month ==m, ]$rain
    fit.mon     <- fitQmap(obsprecip,modprecip.ref,
                           method="QUANT",qstep=0.001)
      qm1       <- doQmap(modprecip.ref,fit.mon,type="tricub")
      
   plot(sqrtquant(modprecip.ref), sqrtquant(obsprecip),type="p", pch=5, col="light grey",ylab='',xlab='')
   lines(sqrtquant(modprecip.ref), sqrtquant(qm1),type="l",col="black")
   mtext(paste(month.abb[m]),side=3,adj=0.05,line=-1.3)
  mtext(expression('P'[m]*'[mm /day]'),side=1, line=2.2,cex=0.9)
  mtext( expression('P'[o]*'[mm /day]'),side=2,line=2,cex=0.9)
   
   
    modprecip.ordered = sort(modprecip.ref)
    plot(modprecip.ordered,type="p" ,pch= 3, col="grey", xaxt='n',ylab='',xlab='')
    obsprecip.ordered = sort(obsprecip)
    lines(obsprecip.ordered,type="p", pch=5, col="light grey")
    qm1.ordered = sort(qm1)
    lines(qm1.ordered,type="l")
    axis(1,labels = c(0,0.5,1) , at=c(0,175,350))
    mtext('Empirical probability',side=1,line=2.2,cex=0.9)
    mtext( 'P [mm /day]',side=2,line=2,cex=0.9)
    
   met_loc_CRUNCEP[which(met_loc_CRUNCEP$date >="2006-01-01" & met_loc_CRUNCEP$month == m),]$rain_update <- qm1
   
  
   plot(cumsum(obsprecip),type="p", pch=5, col="light grey",ylab='',xlab='')
   lines(cumsum(modprecip.ref),type="p" ,pch= 3, col="grey")   
   lines(cumsum(qm1),type="l")
   mtext('days',side=1,line=2.2,cex=0.9)
    mtext( 'P [mm /day]',side=2,line=2,cex=0.9)
    
    if(m==6 |m ==12){
      legend("bottomright", bty="n",lty=c(NA,NA,1), pch=c(5,3,NA), legend=c(expression('P'[o]),expression('P'[m]), expression('P'[o]*'= h(P'[m]*'), dec')),col=c("grey","grey","black"),cex=1.2)
    }
   # modprecip     <- data.frame(rain_update=coredata(qm1),date=date(modprecip))
  # modprecip.ref <- data.frame(rain=coredata(modprecip.ref),date=date(modprecip.ref))
  # obsprecip     <- data.frame(rain_obs=coredata(obsprecip),date=date(obsprecip))
  
}
  dev.off()
} # end split image


```

```{r plot_downscaling_precip_parameterisation}

jpeg(filename='data/APPDX_Downscaling_precip_parameterisation.jpeg')

# omit the last decade, as not used in qqmapping to generate model forcing:
precip.mod           <-  met_loc_CRUNCEP[which(met_loc_CRUNCEP$date < "2006-01-01"),] %>%
                           group_by(month) %>%
                           dplyr::summarize(mod_precip = sum(rain)/105,ds_precip = sum(rain_update)/105)
precip.obs           <- met_loc_obs  %>%
                           group_by(month) %>%
                           dplyr::summarize(obs_precip = sum(rain)/11)
precip <- cbind(precip.mod, precip.obs$obs_precip)


rownames(precip)<-month.abb
precip$month<-NULL
colnames(precip) <-c('default','downscaled','observed')


################## -  PLOTTING - ##################
# precipitation bias correction
par(mfrow=c(2,2),mar=c(3,3,2,2),oma=c(1,1,1,1))

pr_a2011_obs  <- met_loc_obs[met_loc_obs$year>=2011,]$rain
pr_a2011_CRUNCEP <- met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2011-01-01",]$rain
obsprecip <- met_loc_obs$rain
modprecip <- met_loc_CRUNCEP[met_loc_CRUNCEP$date>="2006-01-01",]$rain
modprecip_past<-met_loc_CRUNCEP[met_loc_CRUNCEP$date <"2005-12-31",]$rain


# http://www.imsbio.co.jp/RGM/R_rdfile?f=qmap/man/qmap-package.Rd&d=R_CC
# adjust parameters to get frequency distribution correct:
#par(mfrow=c(2,2),mar=c(0.5,2,0.5,0.5))

# including the reference period increases the yearly average precipiration by 10 mm, but not enough.
# maybe the removal of the maximum value helps adjust the transformation function.
modprecip_past_inclrefP<-met_loc_CRUNCEP$rain

modprecip_past_inclrefP<-met_loc_CRUNCEP$rain
modprecip_past_inclrefP[c(40500)]<-102.97

## call to fitQmapPTF and doQmapPTF
qm1.fit <- fitQmap(obsprecip,modprecip,
              method="PTF",
              transfun="expasympt",
              cost="RSS",wett.day=TRUE)
qm1 <- doQmap(modprecip,qm1.fit)

## call to fitQmapDIST and doQmapDIST
qm2.fit <- fitQmap(sqrt(obsprecip),sqrt(modprecip),
              method="DIST",qstep=0.001,
              transfun="berngamma")
qm2 <- doQmap(sqrt(modprecip),qm2.fit)^2

## call to fitQmapRQUANT and doQmapRQUANT
qm3.fit <- fitQmap(obsprecip,modprecip,
              method="RQUANT",qstep=0.01)
qm3 <- doQmap(modprecip,qm3.fit,type="linear")


## call to fitQmapQUANT and doQmapQUANT
qm4.fit <- fitQmap(obsprecip,modprecip,
              method="QUANT",qstep=0.01)
qm4 <- doQmap(modprecip,qm4.fit,type="tricub")

## call to fitQmapSSPLIN and doQmapSSPLIN
qm5.fit <- fitQmap(obsprecip,modprecip,qstep=0.01,
                   method="SSPLIN")
qm5 <- doQmap(modprecip,qm5.fit)

sqrtquant <- function(x,qstep=0.001){
  qq <- quantile(x,prob=seq(0,1,by=qstep))

}


  plot(sqrtquant(modprecip),
       sqrtquant(obsprecip),pch=19,col="gray",ylab='P_m [mm /day]',xlab= 'P_m [mm /day]')
  lines(sqrtquant(modprecip),
        sqrtquant(qm1),col=1)
  lines(sqrtquant(modprecip),
        sqrtquant(qm2),col=2)
  lines(sqrtquant(modprecip),
        sqrtquant(qm3),col=3)
  lines(sqrtquant(modprecip),
        sqrtquant(qm4),col=4)
  lines(sqrtquant(modprecip),
        sqrtquant(qm5),col=5)
mtext("A)",adj=0.03,line=-2,cex=1.2)
mtext(expression('P'[m]*'[mm /day]'),side=1, line=2.2)
mtext( expression('P'[o]*'[mm /day]'),side=2,line=2)
legend("bottomright",
       legend=c("PTF","DIST","RQUANT","QUANT","SSPLIN"),
       lty=1, col=1:5,bty="n")

# map all successfully adjusted  months

## example month - mapping for within chapter, all other months are in appendix

# quantile-quantile plot 
   m = 5
  modprecip<- met_loc_CRUNCEP[which(met_loc_CRUNCEP$date >="2006-01-01" & met_loc_CRUNCEP$month == m),]$rain
  obsprecip     <- met_loc_obs[met_loc_obs$month ==m, ]$rain
    fit.mon     <- fitQmap(obsprecip,modprecip.ref,
               method="QUANT",qstep=0.001)
      qm1       <- doQmap(modprecip.ref,fit.mon,type="tricub")
    plot(sqrtquant(modprecip), sqrtquant(obsprecip),type="p", pch=5, col="light grey")
    lines(sqrtquant(modprecip), sqrtquant(qm1),type="l",col="green")

   # quatile-quantile plot century - show that the transfpormation function maps 
   # the past century climate on to the relationship between the modelled and observed reference period.
   # plot(sqrtquant(modprecip), sqrtquant(obsprecip),type="p", pch=5, col="light grey")
     #lines(sqrtquant(modprecip_past), sqrtquant(qm2),type="l",col="black")
     legend("bottomright", bty="n", lty=c(NA,1),pch=c(5,NA), legend=c("data", expression('P'[o]*'= h(P'[m]*'), dec')),col=c("grey","green"))
   
   mtext("B)",adj=0.03,line=-2,cex=1.2)
   mtext(expression('P'[m]*'[mm /day]'),side=1, line=2.2)
   mtext( expression('P'[o]*'[mm /day]'),side=2,line=2)


    # empirical cumulative density function for the reference period:
    modprecip.ordered = sort(modprecip)
    plot(modprecip.ordered,type="p" ,pch= 3, col="grey",xaxt= 'n')
    obsprecip.ordered = sort(obsprecip)
    lines(obsprecip.ordered,type="p", pch=5, col="light grey")
    qm1.ordered = sort(qm1)
    lines(qm1.ordered,type="l")
    mtext("C)",adj=0.03,line=-2,cex=1.2)
    legend("topright", bty="n",lty=c(NA,NA,1), pch=c(5,3,NA), legend=c(expression('P'[o]),expression('P'[m]), expression('P'[o]*'= h(P'[m]*')')),col=c("grey","grey","black"))
    axis(1,labels = c(0,0.5,1) , at=c(0,175,350))
    mtext('Empirical probability',side=1, line=2.2)
    mtext( 'P [mm /day]',side=2,line=2)
    

barplot(t(precip), beside=TRUE, cex.axis = 1.2, cex.names = 1.5,
         las=2 ,ylim=c(0,220), col=c(CRUNCEP_col,ds_col,obs_colour))
box(bty="l")

mtext(side=2,"Mean precipitation [mm]", line=3)
mtext("D)",adj=0.03,line=-2,cex=1.2)

dev.off()
```



#HYBRID9 verification


##Soil moisture

```{r import_sw_LT_observations}
file<-'data/OBS_SW_N08.csv'

LT_in<-read.csv(file, header=TRUE)


# create as timestamp
colnames(LT_in)[1]<-"TIMESTAMP"
LT_in$TIMESTAMP<-as.POSIXct(strptime(LT_in$TIMESTAMP, "%d/%m/%Y %H:%M"),tz= "CET")
# filter out timestamps without observations
LT<-LT_in[which(!is.na(LT_in$TIMESTAMP)),]

# observed soil moisture at LT
sw_LT<-LT_in[which(!is.na(LT_in$TIMESTAMP)),]
sw_LT$mean_10cm <- rowMeans(subset(LT, select = c(Depth_10,Depth_10.1,Depth_10.2,Depth_10.3,Depth_10.4)), na.rm = TRUE)/100 # change unit to mm/mm
sw_LT$mean_70cm <- rowMeans(subset(LT, select = c(Depth_70,Depth_70.1,Depth_70.2,Depth_70.3,Depth_70.4)), na.rm = TRUE)/100 # change unit

# turn NANs into NAs in the mean values to avoid issues downstream in the correlation analysis :
sw_LT[is.nan(sw_LT$mean_70cm),]$mean_70cm<-NA
sw_LT[is.nan(sw_LT$mean_10cm),]$mean_10cm<-NA

# Get upper and lower observed soil moisture range for subsequent plotting
sw_LT<-transform(sw_LT, SD70cmmin=apply(subset(sw_LT, select = c(Depth_70,Depth_70.1,Depth_70.2,Depth_70.3,Depth_70.4)),1, min, na.rm = TRUE))
sw_LT<-transform(sw_LT, SD70cmmax=apply(subset(sw_LT, select = c(Depth_70,Depth_70.1,Depth_70.2,Depth_70.3,Depth_70.4)),1, max, na.rm = TRUE))
sw_LT<-transform(sw_LT, SD10cmmin=apply(subset(sw_LT, select = c(Depth_10,Depth_10.1,Depth_10.2,Depth_10.3,Depth_10.4)),1, min, na.rm = TRUE))
sw_LT<-transform(sw_LT, SD10cmmax=apply(subset(sw_LT, select = c(Depth_10,Depth_10.1,Depth_10.2,Depth_10.3,Depth_10.4)),1, max, na.rm = TRUE))

# unit conversion for upper and lower boundary of observations.
sw_LT$SD_70Upper <- sw_LT$SD70cmmax/100
sw_LT$SD_70Lower <- sw_LT$SD70cmmin/100
sw_LT$SD_10Upper <- sw_LT$SD10cmmax/100
sw_LT$SD_10Lower <- sw_LT$SD10cmmin/100
  
```

```{r read_in_soil_vvol_HYBRID9}
# obtain soil mositure for TRW - soil moisture correlations:

# using downscaled CRUNCEP forcing for entire century
sm_ds_cent           <- read.table('data/HYBRID9_SW_DS.dat',header=FALSE)
colnames(sm_ds_cent) <- c("jyear"	, "DOY"	,"iT"	,"NS", 	"b1","b2","b3","b4","b5","b6","b7","b8")
sm_ds_cent           <- create_Timestamp(sm_ds_cent,sub=FALSE)
sm_ds_cent$TIMESTAMP <- as.POSIXct(paste(sm_ds_cent$Date, paste(sm_ds_cent$hour,sm_ds_cent$minute, as.character(00)
                                                                ,sep=":"),sep=" "), format="%Y-%m-%d %H:%M:%S")  # necessary for plotting further down


#using forcing from local met observations, where possible
sm_obs_cent              <- read.table('data/HYBRID9_SW.dat',header = FALSE)
colnames(sm_obs_cent)    <- c("jyear"	, "DOY"	,"iT"	,"NS", 	"b1","b2","b3","b4","b5","b6","b7","b8")
sm_obs_cent              <- create_Timestamp(sm_obs_cent,sub=FALSE)
sm_obs_cent$TIMESTAMP    <- as.POSIXct(paste(sm_obs_cent$Date, paste(sm_obs_cent$hour,sm_obs_cent$minute, as.character(00)  
                                                                     ,sep=":"),sep=" "), format="%Y-%m-%d %H:%M:%S")  # necessary for plotting further down

#subset simulated soil moisture to last decade
soil_moisture_downscaled_sel <- sm_ds_cent[sm_ds_cent$year >=2007,]
soil_moisture_loc_obs_sel    <- sm_obs_cent[sm_obs_cent$year>=2007,]

# observed soil moisture at N08
LT   <- sw_LT

#######simulated soil mositure at N08
soil_moisture_loc_obs    <- soil_moisture_loc_obs_sel[5:18]# decade , 30 minute
soil_moisture_downscaled <- soil_moisture_downscaled_sel[5:18]# decade
sm_obs_cent              <- sm_obs_cent[5:18] # century, daily
sm_ds_cent               <- sm_ds_cent[5:18]  # century

```

```{r plot_sw_with_upper_lower_boundary}
jpeg(filename='figures/Figure_2.jpeg',units="mm", width=80, height=80, res=1000)

par(oma=c(2,2,0.1,0.5),mar=c(0,2,0,0),mfrow=c(2,1),ps = 8, cex = 1)

#make period for plotting flexible
uyear <- 2015
lyear <- 2010
# check how the local and the downscaled sw both make sense for layer 2 and 6
 #sw_LT#<-LT
sw_loc_obs <- soil_moisture_loc_obs # with local forcing

sw_ds <- soil_moisture_downscaled #+ calib # with downscaled grid forcing

# Layer 2:
sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_10Lower[which(is.infinite(sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_10Lower))] <- NA

 sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_10Upper[which(is.infinite(sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_10Upper))] <- NA
 
  sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_70Upper[which(is.infinite(sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_70Upper))] <- NA
 
  sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_70Lower[which(is.infinite(sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_70Lower))] <- NA
  
  
#prepare observations,with upper and lower bound shading:
polys<-polylims(sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_10Upper,sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_10Lower)


plot(sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$mean_10cm,col="grey",type="l",ylim=c(0,0.35),axes=F)

for(i in 1:length(polys)){
	polygon(polys[[i]]$x, polys[[i]]$y, col = rgb(0,0,1,0.5), border = NA)
}


lines(sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$mean_10cm,
      col="purple",
      type="l",ylim=c(0,0.35))
  
lines(sw_loc_obs[year(sw_loc_obs$TIMESTAMP)<=uyear &year(sw_loc_obs$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_loc_obs[year(sw_loc_obs$TIMESTAMP)<=uyear &year(sw_loc_obs$TIMESTAMP)>=lyear ,"b2"],  # matlines
          xaxs="i",
          type = "l", 
          xlab=NA,
          col=obs_colour,
          xaxt="n", 
          ylab = layermiddle[i],
          ylim=c(0,ymax),
           lty=2)

lines(sw_ds[year(sw_ds$TIMESTAMP)<=uyear &year(sw_ds$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_ds[year(sw_ds$TIMESTAMP)<=uyear &year(sw_ds$TIMESTAMP)>=lyear ,"b2"],  # matlines
          xaxs="i",
          type = "l", 
          xlab=NA,
          col=ds_col,
          xaxt="n", 
          ylab = layermiddle[i],
          ylim=c(0,ymax),
           lty=2)

box()
# set ticks
axis.POSIXct(side = 1, tck = -.015,at =c("2010-01-01", "2011-01-01", "2012-01-01", "2013-01-01", "2014-01-01", "2015-01-01", "2016-01-01"), label=NA)
axis(side = 2, tck = -.015, labels = NA)
# set tick labels
axis(side = 2, lwd = 0, line = -0.8, las = 1)
axis.POSIXct(side = 1, at =c("2010-01-01", "2011-01-01", "2012-01-01", "2013-01-01", "2014-01-01", "2015-01-01", "2016-01-01"),lwd=0, line=-1)

# Layer 6:
#prepare observations
  polys<-polylims(sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_70Upper,sw_LT[year(sw_LT$TIMESTAMP)<=uyear&year(sw_LT$TIMESTAMP)>=lyear ,]$SD_70Lower)


plot(sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$mean_70cm,col="grey",type="l",ylim=c(0,0.35),axes=F)

for(i in 1:length(polys)){
	polygon(polys[[i]]$x, polys[[i]]$y, col = rgb(0,0,1,0.5), border = NA)
 }
lines(sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_LT[year(sw_LT$TIMESTAMP)<=uyear &year(sw_LT$TIMESTAMP)>=lyear ,]$mean_70cm,col="purple",type="l",ylim=c(0,0.35))
  
lines(sw_loc_obs[year(sw_loc_obs$TIMESTAMP)<=uyear &year(sw_loc_obs$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_loc_obs[year(sw_loc_obs$TIMESTAMP)<=uyear &year(sw_loc_obs$TIMESTAMP)>=lyear ,"b6"],  # matlines
          xaxs="i",
          type = "l", 
          xlab=NA,
          col=obs_colour,
          xaxt="n", 
          ylab = layermiddle[i],
          ylim=c(0,ymax),
           lty=4)
  
lines(sw_ds[year(sw_ds$TIMESTAMP)<=uyear &year(sw_ds$TIMESTAMP)>=lyear ,]$TIMESTAMP,sw_ds[year(sw_ds$TIMESTAMP)<=uyear &year(sw_ds$TIMESTAMP)>=lyear ,"b6"],  # matlines
          xaxs="i",
          type = "l", 
          xlab=NA,
          col=ds_col,
          xaxt="n", 
          ylab = layermiddle[i],
          ylim=c(0,ymax),
           lty=2)
box()
# set ticks
axis.POSIXct(side = 1, tck = -.015,at =c("2010-01-01", "2011-01-01", "2012-01-01", "2013-01-01", "2014-01-01", "2015-01-01", "2016-01-01"), label=NA)
axis(side = 2, tck = -.015, labels = NA)
# set tick labels
axis(side = 2, lwd = 0, line = -0.8, las = 1)
axis.POSIXct(side = 1, at =c("2010-01-01", "2011-01-01", "2012-01-01", "2013-01-01", "2014-01-01", "2015-01-01", "2016-01-01"),lwd=0, line=-1)

mtext(expression("Volumetric water content"~(mm^3~mm^-3)), side=2,line=2.6,adj=-0.7,cex=1)
mtext(side=2,line=1.5,adj=4.5, "10 cm  below ground")
mtext(side=1,"Year",line=1,cex=1)
mtext(side=2,line=1.5, "70 cm  below ground")
mtext(side=3,line=5.9,adj=0.01, "(a)",font=2)
mtext(side=3,line=-1.7,adj=0.01, "(b)",font=2)
legend("bottomleft",legend=c( expression('obs'[forc]),expression('ds'[forc]),'obs'),fill=c(obs_colour,ds_col,"purple"),bty='n',cex=0.8,horiz = TRUE )


dev.off()
```

```{r soil_moisture_all_layers}
jpeg(filename='figures/APPDX_soil_moisture_all_layers.jpeg')

# Simulated volumetric water content at 8 layers. Lables show the middle depth of the simulated layer.
# Modelled volumetric water content in all soil layers, together with observed soil water content, where possible. The model was either forced with downscaled and locally observed forcing.

##########plot all layers with observations
par(oma=c(3,2,2.5,0.5),mar=c(0,4,0,0),mfrow=c(8,1))

lgth<-dim(soil_moisture_downscaled)[1]

# loop over layers and plot
for (i in 1:7) {
 plot(soil_moisture_downscaled$TIMESTAMP,soil_moisture_downscaled[1:lgth,i],  # matlines
          xaxs="i",
          type = "l", 
          xlab=NA,
          col="black",
          xaxt="n", 
          ylab = layermiddle[i],
          ylim=c(0,ymax))
if(i==3){
    lines(LT$TIMESTAMP,LT$mean_10cm,col=obs_colour,lwd=1.5)

}
  if(i==6){ # layer 6 in HYBRID : midtpoint: 661.0. upper:-493. lower: -829
    lines(LT$TIMESTAMP,LT$mean_70cm,col=obs_colour,lwd=1.5)
   
#plot(LT$TIMESTAMP,LT$mean_70cm/100,col="purple",lwd=1.5)
  }
  #rect(xleft=as.POSIXct("2008-01-01 00:00:00"),xright = as.POSIXct("2012-12-31 00:00:00"),ybottom=0,ytop=0.389, border =NULL, density=20, col="grey", lwd=1.2)
 lines(soil_moisture_loc_obs$TIMESTAMP,soil_moisture_loc_obs[1:lgth,i], col=" dark green")

}
# final plot with time ticks on x-axis
i=8
 plot(soil_moisture_downscaled$TIMESTAMP,soil_moisture_downscaled[1:lgth,i],  # matlines
          xaxs="i",
          type = "l", 
          col="black",
          ylab = layermiddle[i],
          ylim=c(0,ymax))
  
 lines(soil_moisture_loc_obs$TIMESTAMP,soil_moisture_loc_obs[1:lgth,i], col=" dark green")
 
legend("topright",legend=c("forcing_obs","Observed","forcing_ds"), fill=c("dark green",obs_colour,"black"))
dev.off()

###############
## Century:
## used to highlight the coincidence between low values of modelled soil moisture content and the lowest 5 TRWi values.
## plot all layers for th century
jpeg(filename='figures/APPDX_soil_moisture_TScent_lowTRWs.jpeg',units="cm", width=25, height=19, res=300)

par(oma=c(3.2,2,2.5,0.5),mar=c(0,4,0,0),mfrow=c(8,1))

# subset for predawn values each day over century
subs <- subset(sm_obs_cent,format(sm_obs_cent$TIMESTAMP, "%H:%M") == "06:00")
lgth<-dim(subs)[1]

# loop over layers and plot
for (i in 1:7) {
 
 plot(subs$TIMESTAMP,subs[1:lgth,i],  # matlines
          xaxs="i",
          type = "l", 
          xlab=NA,
          col="black",
          xaxt="n", 
          ylab = paste(layermiddle[i]/10,"cm",sep=" "), # convert to cm
          ylim=c(0,ymax))
  
  # the lowest values were identified using: comp_cent_Larch[order(comp_cent_Larch$TRWi)[1:5],]
    rect(xleft=as.POSIXct("1921-01-01 00:00:00"),xright = as.POSIXct("1921-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
      rect(xleft=as.POSIXct("1974-01-01 00:00:00"),xright = as.POSIXct("1974-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
       rect(xleft=as.POSIXct("1947-01-01 00:00:00"),xright = as.POSIXct("1948-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
       rect(xleft=as.POSIXct("1964-01-01 00:00:00"),xright = as.POSIXct("1964-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
              rect(xleft=as.POSIXct("1949-01-01 00:00:00"),xright = as.POSIXct("1949-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
}
              # final plot with time ticks on x-axis
i=8
 plot(subs$TIMESTAMP,subs[1:lgth,i],  # matlines
          xaxs="i",
          type = "l", 
          col="black",
          ylab = paste(layermiddle[i]/10,"cm",sep=" "), # convert to cm
          ylim=c(0,ymax))
 
 # the lowest values were identified using: comp_cent_Larch[order(comp_cent_Larch$TRWi),]
    rect(xleft=as.POSIXct("1921-01-01 00:00:00"),xright = as.POSIXct("1921-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
      rect(xleft=as.POSIXct("1974-01-01 00:00:00"),xright = as.POSIXct("1974-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
       rect(xleft=as.POSIXct("1947-01-01 00:00:00"),xright = as.POSIXct("1948-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
       rect(xleft=as.POSIXct("1964-01-01 00:00:00"),xright = as.POSIXct("1964-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)
              rect(xleft=as.POSIXct("1949-01-01 00:00:00"),xright = as.POSIXct("1949-12-31 00:00:00"),ybottom=0,ytop=0.389, border ="grey", density=0, col="dark grey", lwd=1.2)

  
mtext("Volumetric water content [mm3/mm3] at ", side=2, outer=TRUE,line=0.5)
mtext("Years ", side=1, outer=TRUE,line=2)

legend("topright",legend=c("HYBRID9 soil moisture","Years with lowest observed TRWi"), fill=c("black","dark grey"))
dev.off()

rm(subs)

```

```{r tables_sw_model_comp}

# Comparison of observed and modelled soil moisture_table1: Correlation coefficients (spearman) and root mean square deviations between modelled (L3 and L6) and mean observed (10cm and 70cm) volumetric soil moisture content over different time periods.


# data from 2008
#merge
subs     <- na.omit(merge(soil_moisture_loc_obs[,c("b3","b6","TIMESTAMP")],sw_LT[,c(1,24,25)]))
subs.xts <- xts(subs,subs$TIMESTAMP)

record_metrics<-data.frame(rho=rep(NA,2),RMSE_daily=rep(NA,2),RMSE_monthly=rep(NA,2),RMSE_yearly=rep(NA,2))
rownames(record_metrics)<-c("Layer 3","Layer 6")
#Layer3
record_metrics[1,1]<- as.numeric(cor.test(subs$b3,subs$mean_10cm,method=c("spearman"))$estimate)
record_metrics[1,2]<-rmse(apply.daily(subs.xts$b3,mean),apply.daily(subs.xts$mean_10cm,mean))
record_metrics[1,3]<-rmse(apply.monthly(subs.xts$b3,mean),apply.monthly(subs.xts$mean_10cm,mean))
record_metrics[1,4]<-rmse(apply.yearly(subs.xts$b3,mean),apply.yearly(subs.xts$mean_10cm,mean))

#Layer6 
record_metrics[2,1]<- as.numeric(cor.test(subs$b6,subs$mean_70,method=c("spearman"))$estimate)
record_metrics[2,2]<-rmse(apply.daily(subs.xts$b6,mean),apply.daily(subs.xts$mean_10cm,mean))
record_metrics[2,3]<-rmse(apply.monthly(subs.xts$b6,mean),apply.monthly(subs.xts$mean_10cm,mean))
record_metrics[2,4]<-rmse(apply.yearly(subs.xts$b6,mean),apply.yearly(subs.xts$mean_10cm,mean))

knitr::kable(record_metrics)


```

##Evapotranspiration

```{r read_MODIS_ET}

#var point = ee.Geometry.Point(7.75,46.29);
#// this is not the point from publications but close to where Patrick set the
#// pin down in the google maps image.

ET_Modis_8d <- read.csv('data/MODIS_006_MOD16A2_ET.csv')
# data description: https://lpdaac.usgs.gov/products/mod16a2v006/

# ET: kg/m²/8day
#https://lpdaac.usgs.gov/data/data-citation-and-policies/
colnames(ET_Modis_8d)[2] <- c("ET")
ET_Modis_8d$Date         <- as.Date(gsub("\\_", "-", ET_Modis_8d$imageID))
ET_Modis_8d$ET           <- ET_Modis_8d$ET*0.1 # conversion factor to get to correct unit 


```

```{r read_in_HYBRID_ET}

et_subdaily_in <- read.table(file='data/HYBRID_HYDRO_FLX.dat',header = FALSE)
colnames(et_subdaily_in)<- c("year","DOY", "qflx_surf", "qflx_evap_grnd", "qflx_tran_veg_col", "rsub_top", "qflx_infl", "qflx_rsub_sat","w1","w0","w_soi_liq") # w1= all outgoing water, including total soil water,  w2=all incoming water including total soil water

# leapyears are not considered by HYBRID, so create Date timestamp by hand:
out <- create_Timestamp(et_subdaily_in,sub = FALSE,end = 2016,start=2000)

et_daily <- out[c('Date','qflx_tran_veg_col','qflx_evap_grnd')] #mm/s
#because: in HYBRID9 48 surface timesteps per day, 48 == 1800 s/dt
dt = 1800
et_daily$ET <- (et_daily$qflx_tran_veg_col+et_daily$qflx_evap_grnd)*dt #(mm/day)  

##########
# Obtain 8-day aggregate values comparable to MODIS ET

for (y in 2003:2016){
# use 8 days from previous year to derive first January value:
prev_year_date <- head(et_daily[year(et_daily$Date)==y ,])$Date[1]-7
subs<- et_daily[year(et_daily$Date) == y | et_daily$Date >= prev_year_date ,]

#ep = end point for period which function( mean(x)) should be applied. 
ep <- seq( 8*48, dim(subs)[1], by=8*48 ) # 8*48, to account for the 30-minute time slices in each day
subs <- period.apply( xts(subs$ET,subs$Date), INDEX=ep, FUN=function(x) sum(x) )  
subs <- data.frame(ET_HYBRID = coredata(subs), Date = as.POSIXct(subs))
subs$Date <- as.Date(subs$Date)

if(y==2003){
 out <- subs
}else{
 out  <- rbind(out,subs) # mm/8day
}

}

# final unit conversion
# mm/8 day == kg/m2/8day

## Bring model and data together ##
# Now that HYBRID9 GPP is aggregated for the same dates as we have values from the MODIS product,
# add the MODIS product GPP and the agregated HYBRID9 'product' into the same dataframe:
ET <- merge(out,ET_Modis_8d)
 
# prepare final ET dataframe:
ET <- ET[c('Date','ET_HYBRID','ET')]

```

```{r plot_MODIS_ET_comparison }
#Combined surface evaporation and evapotranspiration of HYBRID9 (8 day aggregate) and the MODIS Evapotranspiration product. MODIS product was taken from (7.77,46.27)"
jpeg(filename='figures/Figure_3.jpeg',units="mm", width=80, height=80, res=1000) # Fig2_MODIS_ET_HYBRID

par(mfrow=c(1,1))
par(oma=c(1,1,0.1,0.1),mar=c(1.1,1.1,0,0),mfrow=c(1,1),ps = 8, cex = 1)

uyear=2016
lyear=2008
plot(ET[year(ET$Date)>=lyear & year(ET$Date)<=uyear,]$Date,ET[year(ET$Date)>=lyear & year(ET$Date)<=uyear,]$ET_HYBRID,type="l",col="green",ylim=c(0,35), ylab="",xlab="",axes=F)
lines(ET[year(ET$Date)>=lyear & year(ET$Date)<=uyear,]$Date,ET[year(ET$Date)>=lyear & year(ET$Date)<=uyear,]$ET,type="l",col="black")
legend('topleft',legend=c("MODIS","HYBRID"),fill=c("black","green"),border = NA,bty="n",ncol=2)
mtext(side=2,cex=1.2,line=1.1,expression("Total evaporation "*'(kg m-2 8days-1)'))
mtext(side=1,cex=1.2,line=1,"Year")

box()
# set ticks
axis.Date(c(2008, 2010, 2012, 2014, 2016) ,side = 1, tck = -.015,at = as.Date(c("2008-01-01", "2010-01-01", "2012-01-01", "2014-01-01","2016-01-01")), label=NA)
axis(side = 2, tck = -.015, labels = NA)
# set tick labels
axis(side = 2, lwd = 0, line = -0.5, las = 1)
axis.Date(c(2008, 2010, 2012, 2014, 2016),side = 1,at = as.Date(c("2008-01-01", "2010-01-01", "2012-01-01", "2014-01-01","2016-01-01")),lwd=0, line=-0.8)



dev.off()

```

```{r, eval=FALSE}
### not included in paper, but maybe the reviewers will ask for it
# correlation analysis of MODIS and HYBRID ET

subs.xts <- xts(ET,ET$Date)
subs.xts <- na.omit(subs.xts)

record_metrics<-data.frame(rho=rep(NA,1),RMSE_daily=rep(NA,1),RMSE_monthly=rep(NA,1),RMSE_yearly=rep(NA,1))


record_metrics[1,1] <- as.numeric(cor.test(ET$ET_HYBRID,ET$ET,method=c("spearman"))$estimate)
record_metrics[1,2] <- rmse(apply.daily(subs.xts$ET_HYBRID,mean),apply.daily(subs.xts$ET,mean))
record_metrics[1,3] <- rmse(apply.monthly(subs.xts$ET_HYBRID,mean),apply.monthly(subs.xts$ET,mean))
record_metrics[1,4] <- rmse(apply.yearly(subs.xts$ET_HYBRID,mean),apply.yearly(subs.xts$ET,mean))

knitr::kable(record_metrics)

rm(subs.xts)


```




#Daily-resolution soil water-TRW correlations

```{r plot_TRW_dynamics_at_NO8}
#this chunk loads in TRW data, builds chronology, plots
#Ring width series (transparent lines) with series median (thick line) (note variable n) and ring width residual chronology for Larch (green) and Spruce (black) around site N08, near Lötschental, Switzerland.


#jpeg(filename='figures/APPDX_TRW_larch_pine.jpeg',units="in", width=6, height=6, res=300)
jpeg(filename='figures/APPDX_TRW_larch_pine.jpeg',units="mm", width=80, height=80, res=300)

#########LARCH
TRW_larch   <- read_excel('data/OBS_TRW_N08.xls',sheet=4,col_types="numeric")
TRW_larch   <- as.data.frame(TRW_larch)

TRW_larch_sub<-TRW_larch[,1:17]
TRW_larch_plot<-TRW_larch_sub[TRW_larch_sub$Year>1796,] # for plotting

TRW_larch_sub_norm<-TRW_larch[,19:35] # normalised values. not currently used
TRW_larch_sub<-TRW_larch_sub[which(TRW_larch_sub$Year>=1901 &TRW_larch_sub$Year<=2015),]

TRW_larch  <-TRW_larch_sub
#spag.plot(TRW_larch_sub_norm[2:17],plot.type="spag")



#########SPRUCE
TRW_spruce  <- read_excel('data/OBS_TRW_N08.xls',sheet=5,col_types="numeric")# read.xls outdated, use read_excel and library(readxl)
TRW_spruce   <- as.data.frame(TRW_spruce)

#TRW_spruce
TRW_spruce_sub<-TRW_spruce[,1:19]
TRW_spruce_plot<-TRW_spruce_sub[TRW_spruce_sub$Year>1853,]


TRW_spruce_sub_norm<-TRW_spruce[,21:39]
#ns<-c("N8DS1"	"N8ADS2",	"N8DS2"	,N8DS3	N8DS4	N8OPS1	N8OPS2	N8OPS3	N8OPS4	N8SHW1	NO8_S1_08	NO8_S1_14	NO8_S2_08	NO8_S2_14	NO8_S3_08	NO8_S3_14	NO8_S4_08	NO8_S4_14)
TRW_spruce_sub <- TRW_spruce_sub[which(TRW_spruce_sub$Year>= 1901 & TRW_spruce_sub$Year<=2015),]

TRW_spruce<-TRW_spruce_sub # normalised values. not currently used
TRW_spruce_sub_norm<-TRW_spruce_sub_norm[which(TRW_spruce_sub_norm$ring.Year>= 1901 & TRW_spruce_sub_norm$ring.Year<=2015),]


# detrend and create chronology, Larch 
TRW_larch_sub_det<-detrend(rwl = TRW_larch_sub, method = "Spline", nyrs=32)#,make.plot = TRUE)#, verbose = TRUE)
TRW_larch_sub_det$Year <- seq(1901,2015,1)
TRW_larch_sub_chr<-chron(TRW_larch_sub_det, prefix = "LAR")

# detrend and create chronology, Spruce
TRW_spruce_sub_det<-detrend(rwl = as.data.frame(TRW_spruce_sub),method = "Spline", nyrs=32)
TRW_spruce_sub_chr<-chron(TRW_spruce_sub_det, prefix = "SPR")

# collect chronology details
TRW_spruce_sub$Chronology <- TRW_spruce_sub_chr$SPRstd
TRW_larch_sub$Chronology  <- TRW_larch_sub_chr$LARstd


for(t in 1:2){
  if(t==1){TRW<-TRW_larch}
  if(t==2){TRW<-TRW_spruce}

ntrees<-dim(TRW)[2]
for (i in 1:(dim(TRW)[1])){
  TRW$max[i]<-max(TRW[i,2:ntrees],na.rm=TRUE)
  TRW$min[i]<-min(TRW[i,2:ntrees],na.rm=TRUE)
  TRW$median[i]<-median(as.numeric(TRW[i,(2:dim(TRW)[2])]),na.rm=TRUE)
}
if(t==1){TRW_larch_sub<-cbind(TRW_larch_sub,TRW$max,TRW$min,TRW$median)}
if(t==2){TRW_spruce_sub<-cbind(TRW_spruce_sub,TRW$max,TRW$min,TRW$median)}
}


########################
# Plotting

# get median values for plotting:
for(t in 1:2){
  if(t==1){TRW<-TRW_larch_plot
  ntrees<-dim(TRW_larch_plot)[2] }
  if(t==2){TRW<-TRW_spruce_plot
  ntrees<-dim(TRW_spruce_plot)[2]}

for (i in 1:(dim(TRW)[1])){
  TRW$max[i]<-max(TRW[i,2:ntrees],na.rm=TRUE)
  TRW$min[i]<-min(TRW[i,2:ntrees],na.rm=TRUE)
  TRW$median[i]<-median(as.numeric(TRW[i,(2:dim(TRW)[2])]),na.rm=TRUE)
}
if(t==1){TRW_larch_plot<-cbind(TRW_larch_plot,TRW$max,TRW$min,TRW$median)}
if(t==2){TRW_spruce_plot<-cbind(TRW_spruce_plot,TRW$max,TRW$min,TRW$median)}
}

par(mfrow=c(2,2),mar=c(1,1.5,0.2,0.5),oma=c(2,1,0,0),cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1,ps=8)
#cex. tipps from http://rfunction.com/archives/2154

#### Plot Larch raw data and median TRW
plot(TRW_larch_plot$Year,TRW_larch_plot$NO8_L1_08...2/100,type="l",ylim=c(0,3.5),xlab='Year',ylab="",col= makeTransparent("dark green",60), axes=F)
for(t in 3:dim(TRW_larch_plot)[2]){
lines(TRW_larch_plot$Year,t(TRW_larch_plot[t]/100),col=makeTransparent("dark green",60))
}
lines(TRW_larch_plot$Year,TRW_larch_plot$`TRW$median`/100,col="dark green")
abline(v=1901,lty=3,col="black",lwd=1)
mtext('Tree ring width [mm]',side=2, line =1.3)


box()
# set ticks
axis(side = 1, tck = -.015, labels = NA)
axis(side = 2, tck = -.015, labels = NA)
# set tick labels
axis(side = 2, lwd = 0, line = -0.8, las = 1)
axis(side = 1, lwd=0,line=-1,las=1)

#### Plot Spruce raw data and median TRW
plot(TRW_spruce_plot$Year,TRW_spruce_plot$N8DS1...2/100,type="l",ylim=c(0,3.5),ylab="",xlab='Year',col=makeTransparent("black",60), axes=F)
for(t in 3:dim(TRW_spruce_plot)[2]){
lines(TRW_spruce_plot$Year,t(TRW_spruce_plot[t])/100,col=makeTransparent("black",60))
}
lines(TRW_spruce_plot$Year,TRW_spruce_plot$`TRW$median`/100,col="black")
abline(v=1901,lty=3,col="black",lwd=1)
mtext('Tree ring width [mm]',side=2, line =1)

box()
# set ticks
axis(side = 1, tck = -.015, labels = NA)
axis(side = 2, tck = -.015, labels = NA)
# set tick labels
axis(side = 2, lwd = 0, line = -0.8, las = 1)
axis(side = 1, lwd=0,line=-1,las=1)

#### Plot resulting Chronology of Larch and Spruce
# to keep in same plot
par(new=TRUE,mfrow=c(2,1),mar=c(1,1,0.2,0.5),oma=c(1,1,0,0),cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1,ps=8)

plot(TRW_spruce_sub$Year,TRW_spruce_sub_chr$SPRstd,type="l",ylim = c(0.4,2), ylab='', axes=F)
lines(TRW_spruce_sub$Year,TRW_larch_sub_chr$LARstd, col="dark green")

mtext('Year', side=1,line=1,cex=1)
legend("topleft",legend=c("Larch","Spruce"),fill=c("dark green","black"),bty='n',cex=0.8)
mtext('TRWi',side=2,line=1)

box()
# set ticks
axis(side = 1, tck = -.015, labels = NA)
axis(side = 2, tck = -.015, labels = NA)
# set tick labels
axis(side = 2, lwd = 0, line = -0.8, las = 1)
axis(side = 1, lwd=0,line=-1,las=1)


dev.off()

######Plot requested by reviewer:

jpeg(filename='figures/APPDX_TRW_larch_aginst_pine.jpeg',units="mm", width=80, height=80, res=300)

par(mfrow=c(1,1),mar=c(1,1,0.2,0.5),oma=c(1,1,0,0),cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1,ps=8)
plot(TRW_larch_sub_chr$LARstd,TRW_spruce_sub_chr$SPRstd,type="p",pch=3,ylim = c(0.4,2),xlim= c(0.4,2), ylab='',axes=F)
mtext(side = 1,"Larch TRWi" , line= 1)
mtext(side = 2, "Spruce TRWi", line= 1.2)
axis(side = 1, tck = -.015, labels = NA)
axis(side = 2, tck = -.015, labels = NA)
# set tick labels
axis(side = 2, lwd = 0, line = -0.5, las = 1)
axis(side = 1, lwd=0,line=-0.9,las=1)
box()

dev.off()

# get TRW for only the time during which soil moisture observations are available
comp_cent_Larch  <- data.frame(TRWi=TRW_larch_sub[TRW_larch_sub$Year >=1901,]$Chronology,
                                 TRW_median=TRW_larch_sub[TRW_larch_sub$Year >=1901,]$`TRW$median`,
                                 Year = TRW_larch_sub[TRW_larch_sub$Year >=1901,]$Year)

comp_cent_Spruce <- data.frame(TRWi=TRW_spruce_sub[TRW_spruce_sub$Year >=1901,]$Chronology,
                                 TRW_median=TRW_spruce_sub[TRW_spruce_sub$Year >=1901,]$`TRW$median`,
                                 Year=TRW_spruce_sub[TRW_spruce_sub$Year >=1901,]$Year)

```

```{r process_and_save_swp, eval=FALSE}
# read in HYBRID9 output and subset dataset for predawn and midday values.

swp_subdaily_in <- read.table('data/HYBRID9_SWP.dat',header=FALSE)

# input file truncated ( HYRBID always spun up from 1805), therefore add header manually
colnames(swp_subdaily_in) <- c("jyear"	, "DOY"	,"iT"	,"NS","swp1",	"swp2",	"swp3",	"swp4",	"swp5",	"swp6",	"swp7",	"swp8","beta")
swp_subdaily     <- create_Timestamp(swp_subdaily_in,1,start=1901,end=2016,sub=FALSE)

#remove unneccessary columns and rename
colnames(swp_subdaily)[1] <- "year"
#swp_subdaily$DOY   <- NULL
swp_subdaily$iT    <- NULL
swp_subdaily$NS    <- NULL

swp.predawn <- swp_subdaily[swp_subdaily$hour == 6 & swp_subdaily$minute==0,]
swp.midday  <- swp_subdaily[swp_subdaily$hour == 12 & swp_subdaily$minute == 0,]

save(swp.midday,file='data/HYBRID9_SWP_mdd.Rdata')
save(swp.predawn,file='data/HYBRID9_SWP_pred.Rdata')
#save(swp_subdaily, file='../../2_data/2_interim/swp_subdaily.RData')

```

```{r read_in_swp}

#load(file='data/HYBRID9_SWP_mdd.Rdata')  # variable name: swp.midday
load(file='data/HYBRID9_SWP_pred.Rdata') # variable name: swp.predawn

```

```{r prep_diag_daily_for_intra_month_corr}

# collect variables used for correlation analysis: temp, precip, vvol, swp

# useful for linear correlation analysis:
# http://www.sthda.com/english/wiki/correlation-test-between-two-variables-in-r#pearson-correlation-formula


##################
# obtain temp and precip
diag_daily          <- read.table('data/MIX_bias_red_prc_t_met_for_HYBRID_N08.txt', head=TRUE)
diag_daily          <- diag_daily[diag_daily$year < 2016,]
diag_daily$tair     <- diag_daily$tair -273.15 # convert to °C , doesn't matter really..


# use all predawn swp values from layer 4: 
diag_daily$swp4  <- swp.predawn[swp.predawn$year<2016 & swp.predawn$year>=1901,][,"swp4"]  #predawn
#diag_daily$swp4 <- swp.midday[swp.midday$year<2016 & swp.midday$year>=1901,][,"swp4"]     #midday
  
# remove stark outliers of nonesnsically low soil water potentials, squeues the correlation analysis
a                <- quantile(diag_daily$swp4,c(0.15))[[1]]
diag_daily[diag_daily$swp4 < a,]$swp4 <- a

## vvol water values 
sm_obs_cent              <- read.table('data/HYBRID9_SW.dat',header = FALSE)
colnames(sm_obs_cent)    <- c("jyear"	, "DOY"	,"iT"	,"NS", 	"b1","b2","b3","b4","b5","b6","b7","b8")
# midday, layer 4
#diag_daily$vvol4     <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],12,sub=TRUE,start=1901,end=2015)[,"b4"] 

#predawn, layer 4
diag_daily$vvol4     <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],t=6,sub=TRUE,start=1901,end=2015)[,"b4"] 

```

```{r NORUN_daily_cor_alg, eval= FALSE}

# files generated in this chunk are already in folder /data/<clim_var>_<sp>_predawn_flex_cor.RData
# analysis can be repeated here, or extended with other clim_vars at daily resolution

######### Perform flexible correlation analysis
years           <- c(c(1901:1950),c(1960:2015))# remove decade where there was tree cutting
climate         <- diag_daily
sp              <- "Spruce_predawn"
#sp              <- "Spruce"   
chrono          <- comp_cent_Spruce[comp_cent_Spruce$Year >= 1960 |comp_cent_Spruce$Year <= 1950,]
clim_vars       <- c("tair" , "rain" , "swp4", 'vvol4' )

max_vals_Spruce <- flex_cor(years = years,climate = climate,chrono = chrono,clim_vars = clim_vars,sp = sp)

sp             <- "Larch_predawn"
#sp              <- "Larch"  
chrono          <- comp_cent_Larch[comp_cent_Larch$Year >= 1960 |comp_cent_Larch$Year <= 1950,]
max_vals_Larch  <- flex_cor(years = years,climate = climate,chrono = chrono,clim_vars = clim_vars,sp = sp)

```

```{r plot_flexible_and_monthly_corr_LarchSpruce}

## daily-resolution correlation results are plotted for each variable.
# in a loop-fashion, the file for each variable an species is read in, 
# the matrix transposed to obtain durations,
# the significant levels are filtered
# the image plotted

# preparations for plotting
names_fig <- c('precipitation', 'soil vvol', 'soil water potential', 'temperature')
days_in_month <- c(31,28,31,30,31,30,31,31,30,31,30,31) # leapyears are not considered
ds_in_month_cumsum <- cumsum(days_in_month)
names_get<-c("rain"  ,   "vvol4" ,   "swp4"  ,   "tair") # called swp4, but is using swp3 for spruce

################################################################
#PLOT correlation matrices from flex_cor output (<clim_var>_<sp>_predawn_flex_cor.RData)

# modify heatmap colours
# solution from mgogol on http://seqanswers.com/forums/showthread.php?t=12022
hmcols<-rev(colorRampPalette(c("red","white","blue"))(20))


jpeg(filename='figures/Figure_4.jpeg',units="mm", width=169, height=80, res=1000)

sel <- c(17,19,18,16) # correct order for plotting
par(mfrow = c(2,4),oma=c(2,5,2,2),mar=c(0,0,0,0),cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1,ps=8)
for (species in 1:2){
  
  if(species==1){
    sp <- "Larch"
    #spec <- "Larch"
    spec <- "Larch_predawn"
  }else{
    sp <- "Spruce"
    #spec <- "Spruce"
    spec <- "Spruce_predawn"
  }

for(v in 1:length(sel)){ # plot all selected variables

  load( file=paste(getwd(),"/data/",names_get[v],"_",spec,"_flex_cor.RData",sep=""))
  cr  <- fcorr$coef
  sig <- fcorr$sig
  
 # rearrange matrix for duration plotting:
  n <- ncol(cr)
for(i in 2:ncol(cr)){
  cr[i:n,i] <- cr[1:(n-i+1),i]
  cr[1:(i-1),i] <- NaN
  cr[1:n,i] <- cr[n:1,i]
}

for(i in 2:ncol(sig)){
  sig[i:n,i] <- sig[1:(n-i+1),i]
  sig[1:(i-1),i] <- NaN
  sig[1:n,i] <- sig[n:1,i]
}
  m                         <- cr
  white                     <- m
  white[which((!is.na(m)))] <- 0
  white[which((is.na(m)))]  <- 1
  white[white==0]           <- NA
  
# image 1
  image(t(apply(m, 2, rev)),breaks=c(seq(-0.5,0.5,0.05)),col=hmcols, axes=FALSE) #xaxt="n", yaxt = "n"
  par(new=TRUE)
  
# prep image 2, filter out the non-significant values:
    if(names_get[v]=="temp"){# relax significance level
    m[sig>=0.1]         <- NA
  }else{
    m[sig>=0.05]       <- NA}
  
  # used in text within the manuscript results section, find the periods and correlation coefficient boundaries
  # set sp=Larch_predawn and v=3 , then print 
  # m[,31] # or 
  # max(m[1:100,31])
  # min(m[1:100,31])
  # max(m,na.rm=TRUE)# max correlation obtained at any period 0.300
  # m[1:2,] # stop after DOY 189
  # m[1:2,] # resues after DOY 231
  
  
  # set sp=Spruce_predawn and v=3 , then print 
  # m[,62] # or 
  # max(m[1:100,62],na.rm=TRUE)  # 78
  # min(m[1:100,31])
  # max(m,na.rm=TRUE)# max correlation obtained at any period 0.350, 
  ## which(m == 0.35234)
  # m[1:2,] # stop after DOY 190

  
   
  m_negative        <- m
  m_negative[which((is.na(m)))]  <- 1
  m_negative[which((!is.na(m)))] <- NA
 
  # image 2 # non-significant areas
  image(t(apply(m_negative,2,rev)), xaxt="n", yaxt = "n",col=makeTransparent("grey",alpha=200))
  par(new=TRUE)

  # image 3 # no values
  image(t(apply(white,2,rev)), xaxt="n", yaxt = "n",col="white")
  par(new=TRUE)
  
  plotlines<-c(0.002739726, 0.002739726 * ds_in_month_cumsum) # plot monthly grid
  
  for( i in 1: length(plotlines)){
    abline(h=plotlines[i],col="white")
    abline(v=plotlines[i],col="white")
  }
  
  if(sp=="Spruce"){# plot DOYs only on bottom of graphs
  # set ticks, avoid overlap in tick labels 
  if(v==1){axis(side = 1, tck = -.015, labels = NA,at= plotlines)
  axis(side = 1, lwd=0,line=-0.8,las=1, at= plotlines[-13], c(1,ds_in_month_cumsum)[-13])
  }
  if(v==3){
   axis(side = 1, tck = -.015, labels = NA,at= plotlines)
  axis(side = 1, lwd=0,line=-0.8,las=1, at= plotlines[-1], c(ds_in_month_cumsum))  
  }
  if(v==4 |v==2){
  axis(side = 1, tck = -.015, labels = NA,at= plotlines)
  axis(side = 1, lwd=0,line=-0.8,las=1, at= plotlines, c(1,ds_in_month_cumsum))
  }
  }
  if(v==1){
    if(sp=="Larch"){
    axis(side = 2, tck = -.015, labels = NA, at= plotlines)
    # set tick labels
    axis(side = 2, lwd = 0, line = -0.6, las = 1,at= plotlines[-1], rev(c(1,ds_in_month_cumsum[-12])))}else{
    axis(side = 2, tck = -.015, labels = NA, at= plotlines)
    # set tick labels
    axis(side = 2, lwd = 0, line = -0.6, las = 1,at= plotlines, rev(c(1,ds_in_month_cumsum)))}
  }
  if(v==4){
     if(sp=="Larch"){
     axis(side = 4, tck = -.015, labels = NA,at= plotlines)
     # set tick labels
     axis(side = 4, lwd = 0, line = -0.6, las = 1,at= plotlines[-1], rev(c(1,ds_in_month_cumsum[-12])))}else
     axis(side = 4, tck = -.015, labels = NA,at= plotlines)
     # set tick labels
     axis(side = 4, lwd = 0, line = -0.6, las = 1,at= plotlines, rev(c(1,ds_in_month_cumsum)))
  }
  ##########?????? 
  ### side=4 still prints overlapping tick label  for DOY 1 and DOY 365
  ##########??????
  
  if(sp=="Larch"){
 
  mtext(paste(names_fig[v]),line=0.2,side=3,cex=0.7,adj=0) # plot variable names on top of graphs
  }
 
 # point_x <- seq(0.001,1,0.002739726)[col_max]
 # point_y <- rev(seq(0.001,1,0.002739726))[row_max]
  point_x=307
  point_y=1
 points(point_x,point_y,pch=8,col="black")
  #mtext()
  
}
  
}  

mtext(side=1,"Start DOY",line = 1,outer=TRUE)
mtext(side=2,"Duration (days)",line = 2.2,adj=0.5,outer=TRUE)
mtext("Larch" ,side=3, line=-2.2, adj=-0.10,outer=TRUE)
mtext("Spruce" ,side=3, line=-18.1,adj=-0.11,outer=TRUE)

mtext("(a)" ,side=3, line=-2.25,adj=0.228,font=2,cex=1,outer=TRUE)
mtext("(b)" ,side=3, line=-2.25,adj=0.48,font=2,outer=TRUE)
mtext("(c)" ,side=3, line=-2.25,adj=0.735,font=2,outer=TRUE)
mtext("(d)" ,side=3, line=-2.25,adj=0.989,font=2,outer=TRUE)

mtext("(e)" ,side=3, line=-12.1,adj=0.228,font=2,outer=TRUE)
mtext("(f)" ,side=3, line=-12.1,adj=0.48,font=2,outer=TRUE)
mtext("(g)" ,side=3, line=-12.1,adj=0.735,font=2,outer=TRUE)
mtext("(h)" ,side=3, line=-12.1,adj=0.989,font=2,outer=TRUE)


#legend
 #par(fig = c(0.85,1, 0.75, 1),new = T)
# rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "white")
 par(fig = c(0.60,0.98, 0.01, 0.40),new = T)
  #par(fig = c(0.5,0.75, 0.75, 1),new = T) 
  image.plot(t(apply(m, 2, rev)), xaxt="n", yaxt = "n",breaks=c(seq(-0.5,0.5,0.05)),col=hmcols,legend.only = TRUE)
################################################################
################################################################
################################################################
dev.off()

```



#Growth period selection

```{r xylogenesis_LT}

# Daily resolution soil-water TRW correlations (above) give an end date for growth (i.e. after which TRW and soil mositure no longer correlate) 
# GDD - modelling for each year of the century provides a plaubsible start date for growth.
# Xylogenesis from 3 years of observations provides shot-term information on recent start and end dates.

#####################
# Xylogenesis:
# find the start and end DOY for growth based on wood growth monitoring data.



collect_CDs <- list()
# this format cannot be easily processed with the CAVIAR package
# so reformat after reading in
##########2008
xylo_2008      <- read.csv('data/OBS_Cell_counts_2008.csv',sep = ";")
colnames(xylo_2008)[7:11] <- c("CZ","EZ","WZ","MZ","PR")

xylo_2008$Date <- as.Date(xylo_2008$Date, tryFormats=c("%d.%m.%Y"))
xylo_2008$Year <- lubridate::year(xylo_2008$Date)
xylo_2008$DY   <- as.integer(strftime(xylo_2008$Date, format = "%j"))

L_xylo_2008   <- xylo_2008[xylo_2008$Species=="L",]
S_xylo_2008   <- xylo_2008[xylo_2008$Species=="S",]

##########2009
xylo_2009<-read.csv('data/OBS_Cell_counts_2009.csv',sep = ";")[1:12]
colnames(xylo_2009)[7:11] <- c("CZ","EZ","WZ","MZ","PR")

xylo_2009$Date <- as.Date(xylo_2009$Date, tryFormats=c("%d.%m.%Y"))
xylo_2009$Year <- lubridate::year(xylo_2009$Date)
xylo_2009$DY   <- as.integer(strftime(xylo_2009$Date, format = "%j"))

end_Larch_2009   = min(xylo_2009[which(xylo_2009$Species=="L" & xylo_2009$EZ == 0 & xylo_2009$DY >= 150),]$DY) # first day no record of EZ

L_xylo_2009   <- xylo_2009[xylo_2009$Species=="L",]
S_xylo_2009   <- xylo_2009[xylo_2009$Species=="S",]


##########2010
xylo_2010<-read.csv('data/OBS_Cell_counts_2010.csv',sep = ";")
colnames(xylo_2010)[7:11] <- c("CZ","EZ","WZ","MZ","PR")

xylo_2010$Date <- as.Date(xylo_2010$Date, tryFormats=c("%d.%m.%Y"))
xylo_2010$Year <- lubridate::year(xylo_2010$Date)
xylo_2010$DY   <- as.integer(strftime(xylo_2010$Date, format = "%j"))

end_Larch_2010   = min(xylo_2010[which(xylo_2010$Species=="L" & xylo_2010$EZ == 0 & xylo_2010$DY >= 150),]$DY) # first day no record of EZ

L_xylo_2010   <- xylo_2010[xylo_2010$Species=="L",]

S_xylo_2010   <- xylo_2010[xylo_2010$Species=="S",]

# combine
Larch_Xylo <- rbind(L_xylo_2010,L_xylo_2008,L_xylo_2009)
Spruce_Xylo <-rbind(S_xylo_2010,S_xylo_2008,S_xylo_2009)
Both_Xylo <- rbind(Larch_Xylo,Spruce_Xylo)

# Prepare df for critical date function:
Both_Xylo$RF <- NA
Both_Xylo    <- Both_Xylo[rep(seq_len(nrow(Both_Xylo)), each = 3), ]
Both_Xylo$RF <- rep(c(1,2,3),times=416)
Both_Xylo    <- Both_Xylo[order(Both_Xylo$DY),]

# To plot values of cell dynamics of all three years 2008-10 in one graph:
i=1
Both_Xylo$Sample <-NA
for(d in unique(Both_Xylo$DY)){ #seem to be different amount of samples, so use loop
Both_Xylo[Both_Xylo$DY==d,]$Sample <-i
 i=i+1
}
Both_Xylo$Year   <- '2008 to 10'
Both_Xylo$Date   <- NULL
Both_Xylo$Slope  <- NULL
LX               <- data.frame(Site=Both_Xylo$Site, Year= Both_Xylo$Year, Species=Both_Xylo$Species,Tree = Both_Xylo$Tree, Sample=Both_Xylo$Sample,DY=Both_Xylo$DY,RF=Both_Xylo$RF,CZ=Both_Xylo$CZ,EZ=Both_Xylo$EZ,WZ=Both_Xylo$WZ,MZ=Both_Xylo$MZ,PR=Both_Xylo$PR)

# to identify plausible growth period:
#CD_L <- computeCriticalDates(na.omit(LX[LX$Species=="L",]))
#save(CD_L,file = 'data/OBS_Critical_Dates_Larch.RData')
load(file='data/OBS_Critical_Dates_Larch.RData')
start_Larch <- round(mean(CD_L$bE))  - 7 # -7 days of cambial activity
end_Larch   = min(LX[which(LX$Species=="L" & LX$EZ == 0 & LX$DY >= 150),]$DY) # first day no record of EZ

#CD_S <- computeCriticalDates(na.omit(LX[LX$Species=="S",]),plot=FALSE)
#save(CD_L,file = 'data/OBS_Critical_Dates_Spruce.RData')
load(file='data/OBS_Critical_Dates_Spruce.RData')
start_Spruce <- round(mean(CD_S$bE)) - 7 # -7 days of cambial activity
end_Spruce   = min(LX[which(LX$Species=="S" & LX$EZ == 0 & LX$DY >= 150),]$DY) # first day no record of EZ


################### Use Growing degree days function from CLM to see when growth might have started
#https://www.rdocumentation.org/packages/chillR/versions/0.70.21/topics/GDD

library(lattice)
daily_met <- read.table('data/MIX_bias_red_prc_t_met_for_HYBRID_N08.txt', header=TRUE)
#colnames(subdaily_met) <- c('year','day','nth','tair','lwdown','swdown','qair','press','rain')
daily_met$tair <- daily_met$tair-273.15 # to °C
Date <-  as.Date(daily_met$day-1, origin = paste(daily_met$year,"-01-01",sep=""))
daily_met.xts <- (xts(daily_met,Date))
yearly <- apply.yearly(daily_met.xts$rain,sum)

# equivalent to MAXVAL in FORTRAN
check_max<-function(z){
  if(z <= 0){
    o=0 }else{
      o=z}
  return(o)
} 
############################
jpeg(filename = 'figures/APPDX_CRITICALDATES_growth_start_density.jpeg',units="in", width=5, height=5, res=300)
# growth start:
# Growing degree days:
 # based on CLM
# https://escomp.github.io/ctsm-docs/doc/build/html/tech_note/DGVM/CLM50_Tech_Note_DGVM.html
#2.29.12.1.1. Cold Leaf-out timing
# first day above 5 °C
# sum of mean daily temperatures (T_{day} ^{o}C) above a given threshold T_{g} (0 ^{o}C).
Tbase=5
i=1
DOY_growth_start<-rep(NA,length(1901:2016))
for (y in 1901:2016){
   DOY_growth_start[i] <-  which(cumsum(data.frame(t(sapply(daily_met[daily_met$year == y,]$tair-Tbase,check_max))) )>=5)[1]
   i=i+1
}
d <- density(DOY_growth_start)

maxd<-modes(d)

plot(d$x,d$y, type="l",xlab='',ylab='')
# used dwws answer on https://stackoverflow.com/questions/36948624/shade-area-under-a-curve
polygon(c(d$x, 0), c(d$y, 0), col=makeTransparent("light blue",alpha=70))
abline(v=maxd$x[3],lwd=1.4, col="dark grey")

abline(v=start_Larch,col=Larch_col,lwd=1.4)
abline(v=start_Spruce,col=Spruce_col,lwd=1.4)
mtext("DOY",side=1,line=2)
mtext("probability density",side=2,line=2.2)

dev.off()

```

```{r plot_Woodformation_Dynamics, eval=FALSE}
#To reproduce this plot, the CaviaR package must be installed. This must be manually installed, as it is no longer available on CRAN

jpeg('figures/APPDX_Xylogenesis_Spruce_08-10_N08.jpeg')
plotWoodFormationDynamics(LX[LX$Species == "S",],level="MvM")
dev.off()

jpeg('figures/APPDX_Xylogenesis_Larch_08-10_N08.jpeg')
plotWoodFormationDynamics(LX[LX$Species == "L",],level="MvM")
dev.off()

```


```{r plot_growth_period}

# figure summarising the evidence used for growing degree day definition

jpeg(filename='figures/Fig1.jpeg',units="mm", width=80, height=80, res=1000)

par(mfrow=c(1,1),mar=c(2,2.5,1,0.2),cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1,ps=8)

# GDD model:
plot(d$x[which(d$x>=1)],d$y[which(d$x>=1)], type="l",xlab='',ylab='',xlim=c(0,365),axes=F)
polygon(c(d$x, 0), c(d$y, 0), col=makeTransparent("light blue",alpha=65))


# Xylogenesis: 
lines(start_Larch:end_Larch, rep(0.007,66),col=Larch_col,lwd=2)
lines(start_Spruce:end_Spruce, rep(0.009,81),col=Spruce_col,lwd=2)

lines(rep(start_Larch,3),seq((0.007-0.001),(0.007+0.001),0.001),col="red",lwd=2)
lines(rep(end_Larch,3),seq((0.007-0.001),(0.007+0.001),0.001),col="red",lwd=2)

lines(rep(start_Spruce,3),seq((0.009-0.001),(0.009+0.001),0.001),col="red",lwd=2)
lines(rep(end_Spruce,3),seq((0.009-0.001),(0.009+0.001),0.001),col="red",lwd=2)

abline(v=maxd$x[3],lwd=1.4, col="dark grey")

#########################
# Daily corr analysis:

#rect(xleft=29,xright = 190,ybottom=0.01810658-(0.01810658/5*2),ytop=0.01810658-0.01810658/5-0.0005, border =makeTransparent(Larch_col,alpha=170), density=100, col=makeTransparent(Larch_col,alpha=20), lwd=1.2)

lines(rep(190,4),seq((0.01810658-(0.01810658/5*2)),(0.01810658-0.01810658/5-0.0005),0.001),col=Larch_col,lwd=2)

lines(rep(198,4),seq((0.01810658-0.01810658/5),(0.01810658),0.001),col=Spruce_col,lwd=2)

#rect(xleft=62,xright = 198,ybottom=0.01810658-0.01810658/5,ytop=0.01810658, border =makeTransparent(Spruce_col,alpha=170), density=100, col=makeTransparent(Spruce_col,alpha=20), lwd=1.2)




rect(xleft=100,xright = 125,ybottom=0.0001,ytop=0.0185, border =makeTransparent("red",alpha=170), density=100, col=makeTransparent("red",alpha=10), lwd=1.2,lty=3)


rect(xleft=175,xright = 205,ybottom=0.0001,ytop=0.0185, border =makeTransparent("red",alpha=170), density=100, col=makeTransparent("red",alpha=10), lwd=1.2,lty=3)

mtext("Daily-resolution \ncorrelation analysis",side=1, line=-10,adj=0.9)
mtext("GDD-model", side=1, line=-2.2,adj=0.05)
mtext("Xylogenesis \nObservations",side=1, line =-5.7, adj=0.9)

mtext("Growth start",side=3,adj=0.22,col="red")
mtext("Growth end",side=3,adj=0.57,col="red")
mtext("Probability density",side=2,line=1.8)
mtext("Day of the year",side=1,line=1)
mtext(expression('t'[1]),side=1, line=-6.3,adj= 0.25,cex=1.3) 
#mtext(side=2, expression('beta'[w[crit]] *' [-]'), line=1)
mtext(expression('t'[2]),side=1, line=-6.3,adj= 0.585,cex=1.3)


axis(side = 1, tck = -.015, labels = NA)
axis(side = 1, lwd=0,line=-0.8,las=1)
axis(side = 2, tck = -.015, labels = NA)
axis(side = 2, lwd=0,line=-0.7,las=1)
box()
dev.off()

```



# Tree growth soil moisture response.

```{r binarrangement_comparison_plot}

jpeg("figures/APPDX_Quantile_arrangement.jpeg")
  par(mfrow=c(1,2),oma=c(1,1,1,1))
    quantiles <- get_bin_bounds(comp_cent_Larch$TRWi,10)
    plot(sort(comp_cent_Larch$TRWi),ylab="",xlab="")
    mtext("A)",side=3, line=1, adj=-0.1)
    for(q in quantiles){
    lines(c(0,120),c(q,q))
    }
    quantiles <- quantile(comp_cent_Larch$TRWi,seq(0,1,0.1))
    plot(sort(comp_cent_Larch$TRWi), xlab="",ylab="")
    mtext("B)",side=3, line=1, adj=-0.1)
    for(q in quantiles){
      lines(c(0,120),c(q,q))
    }
    mtext("ordered TRWi samples",side=1, line =-3,outer=TRUE)
    mtext("TRWi",side=2, line =-0.9,outer=TRUE)
    
dev.off()
        
```

##Optimising for w_crit

```{r optimise_TRW_growth_response}


#######################################################
n_bins <- 11  
Larch_col <- "dark green"
Spruce_col <- "black"
swp.sel <- swp.predawn # swp.sel
#####################################################

######################################################
# decrease w_crit sequentially, then search for 'incr' that maximises R2
# use incr to simulate MEF. (the MEF only helps to visualise a good fit)
# plot MEF, R2 and MEF


for(sp in 1:2){
  if(sp==1){sp<-"Larch"}
  if(sp==2){sp<-"Spruce"}


if(sp == "Larch"){
  #df <- df_Larch
  TRWi_sel  <- data.frame(TRWi=comp_cent_Larch$TRWi, Year=comp_cent_Larch$Year)
  # remove decade where there was tree cutting
  years    <- c(c(1901:1950),c(1960:2015))
  TRWi_sel <- TRWi_sel[TRWi_sel$Year >= 1960 | TRWi_sel$Year <= 1950,]
  start = start_Larch
  end   = end_Larch
}
  
if(sp =="Spruce"){
  #df        <- df_Spruce
  TRWi_sel  <-  data.frame(TRWi=comp_cent_Spruce$TRWi,Year=comp_cent_Spruce$Year)
  # remove decade where there was tree cutting
  years    <- c(c(1901:1950),c(1960:2015))
  TRWi_sel <- TRWi_sel[TRWi_sel$Year >= 1960 | TRWi_sel$Year <= 1950,]
  start    = start_Spruce
  end      = end_Spruce
}

w_crit = c(seq(0,-1.4,-0.01),seq(-1.4,-4.6,-0.2))
  
obj_fun <- rep(NaN,length(w_crit ))
slopes <- rep(NaN,length(w_crit ))


for(i in 1:length(w_crit)){
  
  # for a given w_crit, simulate tree growth for each year in each bin
  output <- simulate_growth_in_bins(start,end,w_crit[i],n_bins=n_bins,TRWi_sel=TRWi_sel,swp_sel=swp.sel,species=sp)
  
  # create mean simulated growth per bin
  x_y   <- simulate_mean_bin_growth(output,sp=sp)
  
  # compare with TRWi derived from observations. We want to maximise the R2
  obj   <- compare_growth(x_y)
  
  obj_fun[i] <- obj[[1]] 
  slopes[i]  <- obj[[2]]
}

cutoff_level <- max(obj_fun)
w_crit_optim <- w_crit[min(which(obj_fun >= cutoff_level))]
slope_optim  <- slopes[min(which(obj_fun >= cutoff_level))]

##########

# now that incr is found, set as fixed, check again for good fit using MEF:
df    <- data.frame(parameter = w_crit, MEF=rep(NA,length(w_crit)))


if(sp=="Larch"){
  # save obj_fun for later plotting in main publication graph
  obj_fun_Larch <-obj_fun
  MEF_Larch <- df
  w_crit_optim_L <-w_crit_optim
  slope_optim_L   <- slope_optim
  for(i in 1:length(w_crit)){
   MEF_Larch[i,2]   <-  c(MEF_optim(start,end,w_crit[i],n_bins=n_bins,TRWi_sel,swp_sel=swp.sel,incr=slope_optim_L,sp="Larch")) # equivalent to finding R2, but now with fixed incr
  }
  
}

if(sp=="Spruce"){
  # save obj_fun for later plotting in main graph
  obj_fun_Spruce <-obj_fun
  MEF_Spruce <- df
  w_crit_optim_S <- w_crit_optim
  slope_optim_S   <- slope_optim
  for(i in 1:length(w_crit)){
    MEF_Spruce[i,2]   <-  c(MEF_optim(start,end,w_crit[i],n_bins=n_bins,TRWi_sel,swp_sel=swp.sel,incr=slope_optim_S,sp="Spruce")) # equivalent to finding R2, but now with fixed incr
  }
 
}
}
# now both incr (=a in manuscript) and w_crit are optimised.

############
```


```{r plot_optimised_TRW_growth_response}

###########
#Figure 5:
jpeg(filename='figures/Figure_5.jpeg',units="mm",width=180, height=80, res=700)

par(mfrow=c(2,4),mar=c(0.5,1.5,0,1),oma=c(2,2.8,1,0),cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1,ps=8)

for(species in 1:2){
  
  if(species==1){
    sp<-"Larch"
    df <- MEF_Larch
    TRWi_sel  <- data.frame(TRWi=comp_cent_Larch$TRWi, Year=comp_cent_Larch$Year)
    # remove decade where there was tree cutting
    years      <- c(c(1901:1950),c(1960:2015))
    TRWi_sel   <- TRWi_sel[TRWi_sel$Year >= 1960 | TRWi_sel$Year <= 1950,]
    w_crit_opt <- w_crit_optim_L
    slope_optim = slope_optim_L # equivalent to 'increment' ( a in manuscript)
    # select for critical period
    start = start_Larch
    end   = end_Larch
    obj_fun <-obj_fun_Larch
    
    ########################    
    quantiles <- get_bin_bounds(TRWi_sel$TRWi,n_bins=n_bins)
    plot(sort(TRWi_sel$TRWi),ylab="",xlab="",xaxt ="n",yaxt="n")
     axis(2, mgp=c(3, .6, 0), tck=-.015, labels=NA)
     axis(side = 2, lwd = 0, line = -0.8, las = 0.1)

    for(q in quantiles){
      lines(c(0,120),c(q,q),col="grey",lty=2)
    }
    mtext(expression('TRW'[i]),side=2, line =0.9)
    
    
    swp <- w_crit
    # To create a continuous colour palette, I used answer from joran on
    # https://stackoverflow.com/questions/9946630/colour-points-in-a-plot-differently-depending-on-a-vector-of-values
    
    #Create a function to generate a continuous color palette
    rbPal <- colorRampPalette(c('orange','red'))
    
    #subset for relevant soil water potential range, around which MEF performs well:
    # split sample for better colouring, to make the useful MEF values more visible
    subs1  <-   df[which(df$MEF<= 0.5 & df$MEF > 0),]
    subs2  <-   df[which(df$MEF>= 0.5),]
    
    subs2$Col <-NA
    # This adds a column of color values
    subs1$Col <- "black" #"#FF0000" # red, for not significant ?!
    subs2[order(subs2$MEF),]$Col  <- rbPal(length(subs2$MEF))[as.numeric(cut(sort(subs2$MEF),breaks = length(subs2$MEF)))]
    subs<-rbind(subs1,subs2)
    
    
    
    ########################    
    plot(w_crit,obj_fun,ylim=c(0.0,max(obj_fun)*1.1),
         xlim=c(0,-3),
         xlab = "",ylab = "",col="light grey",type="p",xaxt ="n",yaxt="n")
    axis(2, mgp=c(3, .6, 0), tck=-.015, labels=NA)
    axis(side = 2, lwd = 0, line = -0.8, las = 0.5)
    # we set a visual cutoff at 0.5 above which  MEF is relevant 
    
    lines(subs$parameter,subs$MEF,ylim=c(0,0.8),xlim=c(0,-3),col=subs$Col, ylab="n", xlab="",type="p")
    points(subs[subs$MEF==max(subs$MEF),]$parameter,max(subs$MEF),col="yellow")
    abline(h=0.5,lty=4)
  
    mtext("Variance explained", side =2 ,line=1)
    numb <- round(max(obj_fun),3)
    numa <- round(slope_optim,3)
    mtext(bquote(R^2~' = '~.(numb)),side=3,line=-1.2, adj=0.9)#line=-2,adj=9.5)
    mtext(bquote(a^'*'~ ' = '~.(numa)),side=3,line=-2.2, adj=0.9)

    ########################
    plot(swp,btc(swp,subs$parameter[1] ),ylim=c(0,1),type="l", xlim=c(0,-3),col=subs$Col[1],xaxt ="n",yaxt="n")
    mtext(side=2, expression('beta'[w[crit]] *' [-]'), line=1)
    axis(2, mgp=c(3, .6, 0), tck=-.015, labels=NA)
    axis(side = 2, lwd = 0, line = -0.8, las = 0.5)
    for (l in 2:length(subs$MEF)){
      lines(swp,btc(swp,subs$parameter[l] ),ylim=c(0,1),type="l", xlim=c(-3,0),col=subs$Col[l])
    }
    lines(swp,btc(swp,w_crit_opt),ylim=c(0,1),type="l", xlim=c(0,-3),col="yellow")
    mtext(side=3,bquote( Psi[crit(growth)] *'='~.(w_crit_opt)),line=-1.9,adj=0.9)
    
    ########################
    if(sp == "Larch"){
      #  collect_samples_per_bin_L <- plot_all_update(start=start,end=end,w_crit=w_crit_optim,n_bins=n_bins,TRWi_sel=TRWi_sel,swp_sel=swp.sel,type='corr',species=sp,incr=slope_optim) #type 
      collect_samples_per_bin_L <- plot_all_update(start=start,end=end,w_crit=w_crit_opt,n_bins=n_bins,TRWi_sel=TRWi_sel,swp_sel=swp.sel,type='growth_model',species=sp,incr=slope_optim) 
      #  collect_samples_per_bin_L <- plot_all_update(start=start,end=end,w_crit=w_crit_optim,n_bins=n_bins,TRWi_sel=TRWi_sel,swp_sel=swp.sel,type='exp_var',species=sp,incr=slope_optim) 

      mtext(expression('Mean binned TRW'[i]),side=2, line =1)
       }
    }
  if(species==2){
    sp<-"Spruce"
   
    #df        <- df_Spruce_predawn
    df        <- MEF_Spruce
    TRWi_sel  <-  data.frame(TRWi=comp_cent_Spruce$TRWi,Year=comp_cent_Spruce$Year)
    # remove decade where there was tree cutting
    years    <- c(c(1901:1950),c(1960:2015))
    TRWi_sel <- TRWi_sel[TRWi_sel$Year >= 1960 | TRWi_sel$Year <= 1950,]
    w_crit_opt <- w_crit_optim_S  
    slope_optim       = slope_optim_S  # equivalent to 'increment'
    # select for critical period
    start    = start_Spruce
    end      = end_Spruce
    obj_fun <- obj_fun_Spruce
  

  quantiles <- get_bin_bounds(TRWi_sel$TRWi,n_bins=n_bins)
  plot(sort(TRWi_sel$TRWi),ylab="",xlab="",xaxt ="n",yaxt="n")
  axis(2, mgp=c(3, .6, 0), tck =-.015, labels = NA)
  axis(side = 2, lwd = 0, line = -0.8, las = 0.1)
  axis(1, mgp=c(3, .6, 0), tck=-0.015, labels = NA)
  axis(side = 1, lwd=0,line=-1,las=1)
  #mtext("A)",side=3, line=-0.7 , adj=-0.22)
  
  for(q in quantiles){
    lines(c(0,120),c(q,q),col="grey",lty=2)
  }
  mtext("Ordered TRWi samples",side=1, line = 1.2)
  mtext(expression('TRW'[i]),  side=2, line = 0.9)
  
  
  swp <- w_crit
  
  # To create a continuous colour palette, I used answer from joran on
  # https://stackoverflow.com/questions/9946630/colour-points-in-a-plot-differently-depending-on-a-vector-of-values
  
  #Create a function to generate a continuous color palette
  rbPal <- colorRampPalette(c('orange','red'))
  
  #subset for relevant soil water potential range, around which MEF performs well:
  # split sample for better colouring, to make the useful MEF values more visible
  subs1  <-   df[which(df$MEF<= 0.5 & df$MEF > 0),]
  subs2  <-   df[which(df$MEF>= 0.5),]
  
  subs2$Col <-NA
  # This adds a column of color values
  subs1$Col <- "black"
  subs2[order(subs2$MEF),]$Col  <- rbPal(length(subs2$MEF))[as.numeric(cut(sort(subs2$MEF),breaks = length(subs2$MEF)))]
  subs<-rbind(subs1,subs2)
  
  plot(w_crit,obj_fun,ylim=c(0.0,max(obj_fun)*1.1),
       xlim=c(0,-3),
       xlab = "",ylab = "",col="light grey",type="p",xaxt ="n",yaxt="n")
  axis(2, mgp=c(3, .6, 0), tck=-.015, labels=NA)
  axis(side = 2, lwd = 0, line = -0.8, las = 0.5)
  axis(1, mgp=c(3, .6, 0), tck=-.015, labels = NA)
  axis(side = 1, lwd=0,line=-1,las=1)
  # we set a visual cutoff at 0.5 above which  MEF is relevant 

  lines(subs$parameter,subs$MEF,ylim=c(0,0.8),xlim=c(0,-3),col=subs$Col, ylab="n", xlab="",type="p")
  points(subs[subs$MEF==max(subs$MEF),]$parameter,max(subs$MEF),col="yellow")
  abline(h=0.5,lty=4)
  mtext(side=1,expression( Psi[crit(growth)] *' [MPa]'),line=1.5)

  mtext("Variance explained", side =2 ,line=1)
  numb <- round(max(obj_fun),3)
  numa <- round(slope_optim,3)
  mtext(bquote(R^2~' = '~.(numb)),side=3,line=-1.2, adj=0.9)#line=-2,adj=9.5)
  mtext(bquote(a^'*'~ ' = '~.(numa)),side=3,line=-2.2, adj=0.9)

  plot(swp,btc(swp,subs$parameter[1] ),ylim=c(0,1),type="l", xlim=c(0,-3),col=subs$Col[1],xaxt ="n",yaxt="n")
  axis(2, mgp=c(3, .6, 0), tck=-.015, labels=NA)
  axis(side = 2, lwd = 0, line = -0.8, las = 0.5)
  axis(1, mgp=c(3, .6, 0), tck=-.015, labels = NA)
  axis(side = 1, lwd=0,line=-1,las=1)
  mtext(side=2, expression('beta'[w[crit]] *' [-]'), line=1)
  for (l in 2:length(subs$MEF)){
    lines(swp,btc(swp,subs$parameter[l] ),ylim=c(0,1),type="l", xlim=c(-3,0),col=subs$Col[l])
  }
  lines(swp,btc(swp,w_crit_opt),ylim=c(0,1),type="l", xlim=c(0,-3),col="yellow")

  mtext( bquote(Psi[s]~'[MPa]'), side= 1,line= 1.5)
  
  mtext(side=3,bquote( Psi[crit(growth)] *'='~.(w_crit_opt)),line=-1.9,adj=0.9)
  
  
  ###########
 
  if(sp == "Spruce"){
    collect_samples_per_bin_S <- plot_all_update(start=start,end=end,w_crit=w_crit_opt,n_bins=n_bins,TRWi_sel=TRWi_sel,swp_sel=swp.sel,type='growth_model',species=sp,incr=slope_optim)
  }
  mtext(expression('Mean binned TRW'[i]),side=2, line =1)
  mtext(expression('TRW'[s]), side =1 ,line=1.5)
  }
}
#####Larch
par(fig = c(0.031,0.18, 0.77, 1),new = T) 
rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "white")
par(fig = c(0.031,0.18, 0.77, 1),new = T) 

plot(collect_samples_per_bin_L,type="h",lwd = 3,ylim=c(0,(max(collect_samples_per_bin_L)+1)),cex.lab=0.96,cex.axis=0.87,ylab= "# samples in bin" ,mgp=c(0.65,0.21,0), tck=-0.05, col="dark grey",xaxt ="n")
axis(1, mgp=c(3, .6, 0), tck=-.035, labels = NA)
axis(side = 1, lwd=0,line=-1.08,las=1,cex.axis=0.9)
mtext(side=1,'bin',line = 0.32, cex= 0.67)

#####Spruce
par(fig = c(0.03,0.18, 0.27, 0.5),new = T) 
rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "white")
par(fig = c(0.03,0.18, 0.27, 0.5),new = T) 

plot(collect_samples_per_bin_S,type="h",lwd = 3,ylim=c(0,(max(collect_samples_per_bin_S)+1)),cex.lab=0.96,cex.axis=0.89,ylab= "# samples in bin" ,mgp=c(0.65,0.21,0), tck=-0.05, col="dark grey",xaxt ="n")
axis(1, mgp=c(3, .6, 0), tck=-.035, labels = NA)
axis(side = 1, lwd=0,line=-1.08,las=1,cex.axis=0.89)
mtext(side=1,'bin',line = 0.32, cex= 0.67)

par(fig = c(0,1,0,1),new = T) 
mtext("(a)",side=3, line=-0.75, adj=-0.32,   font= 2)
mtext("(b)",side=3, line=-0.75, adj= 0.065,  font= 2)
mtext("(c)",side=3, line=-0.75, adj= 0.45,  font= 2)
mtext("(d)",side=3, line=-0.75, adj= 0.832,  font= 2)

mtext("(e)",side=3, line=-11.72,  adj=-0.3,   font= 2)
mtext("(f)",side=3, line=-11.72,  adj= 0.065,  font= 2)
mtext("(g)",side=3, line=-11.72,  adj= 0.45,  font= 2)
mtext("(h)",side=3, line=-11.72,  adj= 0.832,  font= 2)


mtext("Larch"  ,side=3, line=-2,    adj= -0.380, cex= 1)
mtext("Spruce" ,side=3, line=-18.6, adj= -0.389, cex= 1)

dev.off()

#############################################


```



#Comparison with other source and sink functions and observations

```{r plot_other_response_functions}
jpeg(filename='figures/Figure_6.jpeg',units="mm",width=80, height=80, res=700)

insertRow <- function(existingDF, newrow, r) {
  existingDF[seq(r+1,nrow(existingDF)+1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}
#https://stackoverflow.com/questions/11561856/add-new-row-to-dataframe-at-specific-row-index-not-appended
# from Ari B. Friedman
Env2=insertRow(Env2,c(0.002,20,0.005,0.0000000,21000,40,101325),1)
################

# Modelled sources and sinks 
# all models are driven by Env2 environment dataframe from supplementary materials in Eller et al 2020

source_col = "dark green"
Growth_col = "brown"

CLM45_pch = 6
SOX_pch   = 5
CLM5_pch  = 3

Type1_lty = 1
Type2_lty = 2
Type3_lty = 3
Type4_lty = 4

col_SOX     = source_col
lty_SOX     = SOX_pch
col_CLM45   = source_col
lty_CLM45   = CLM45_pch
col_HYBRID9 = source_col
lty_HYBRID9 = 1
HYBRID_pch  = 1
col_CLM5    = source_col
lty_CLM5    = CLM5_pch 

cex_set    = 0.6 # size of organism specific pch-es
lwd_set    = 1
lwd_SA     = 0.9 # line width for (grey) graphs which differ from default parameterisation
par(mar=c(0,0,0,5),mfrow=c(2,1),oma=c(2,2,0.5,0),cex.axis=1, cex.lab=1, cex.main=1, cex.sub=1,ps=8)#,mar = par()$mar + c(0,0,0,7))

##################           
# SOX MODEL 
krc_max  =  0.01 #1/1000
He       = 25 # tree height

# Rosner et al 2018 Table 1 - Hydraulic traits 
# Species                     P50
# Picea abies trunk wood     -2.2
# Larix decidua sapling      -4.4


outPA   <- SOX(P50 = -2.2, 
               krc_max = krc_max, 
               He = He, 
               Photo = Photo, 
               Env = Env2)
outLD   <- SOX(P50 = -4.4, 
               krc_max = krc_max, 
               He = He, 
               Photo = Photo, 
               Env = Env2)

# normalise output
PA_norm <- max(outPA$Aopt)
LD_norm <- max(outLD$Aopt)

# Ps = soil water potential
plot(outPA$Aopt/PA_norm ~ Env2$Ps, 
     type = "l",xlim=c(0,-4),ylim=c(0,1.2),col=col_SOX, xlab="",ylab="",xaxt="n",yaxt="n", yaxs="i", xaxs="i", lty=lty_SOX,lwd=lwd_set)

lines(outLD$Aopt/LD_norm ~ Env2$Ps,
      type = "l",col=col_SOX, xlab="",ylab="",xaxt="n",yaxt="n", yaxs="i", xaxs="i", lty=2,lwd=lwd_set)


##################           
# CLM4.5 

# Parameters for relevant plant functional types from CLM technical note 4.5 Oleson et al 2014
Psi_o    = c(66000,   66000, 1019.72, 101972) * 9.81 / -1.0E6 # convert to MPa
Psi_crit = c(255000, 152958,  203944, 305916) * 9.81 / -1.0E6 # convert to MPa
Sensitivity   <- data.frame(PFT=c("NET_Temperate","Change1","Change2","Change3"), Psi_o = Psi_o,Psi_crit=Psi_crit)

out_CLM   <- list()
# fill beta_soil with different response curves, to simulate stomata and photosynthesis response curve
for(b in 1:4){
beta_soil = btcCLM45(Env2$Ps,
                     Sensitivity$Psi_crit[b],#[pft],
                     Sensitivity$Psi_o[b])#[pft])  # CLM4.5 soil moisture response curve

# run CLM45 Photosynthesis and Stomatal Conductance function
out_CLM[[b]] <- CLM45(forc_pbot = Env2$Pa,
                 Ta = Env2$Ta,
                 cair = Env2$Ca,
                 oair = Env2$Oa,
                 btran = beta_soil)
}      

names(out_CLM) <- Sensitivity$PFT

CLM_normp      <- max(out_CLM$NET_Temperate$Anet)# default NET_ temperate parameterisation
            
# Photosynthesis - soil water response
lines(out_CLM$NET_Temperate$Anet/CLM_normp ~ Env2$Ps, 
      type = "l", xlim=c(0,-4), ylim=c(0,1.2), col=col_CLM45, xlab="",ylab="",xaxt="n",yaxt="n", yaxs="i",
      xaxs="i", lty=lty_CLM45,lwd=lwd_set)

##################           
# HYBRID9 

LAI_ind    = 2 # Needle leaf
crown_area = 
  mean( c(18.85065    ,   18.79010      , 12.03321   ,    19.00207      , 18.63200,
          18.71060 ,      18.60301  ,     18.96607   ,    18.97120)) # mean crown area from plot at Swiss alps

beta_forcingstom = btc_plot(Env2$Ps,-1) # soil moisture effect on stomata
beta_forcingphot = beta_forcingstom     # soil moisture effect on photosynthesis (Vcmax and Jmax)

out_HYBRID9 <- HYBRID9_photosynthesis( Wfac_phot = beta_forcingphot,  # stomatal and non-stomatal soil moisture effects currently have the same parameterisation.
                                       beta = beta_forcingstom,
                                       Nfol = 0.012,
                                       sla = 3.9,
                                       Vcmax_max = 94.0,
                                       Jmax_max = 143.0,
                                       LAI_ind = LAI_ind,
                                       CO2 = Env2$Ca,
                                       tair = Env2$Ta,
                                       Iapar = Env2$IPAR,
                                       VDD = Env2$D,
                                       crown_area = crown_area,
                                       Pa = Env2$Pa,
                                       Oi = Env2$Oa)


lines(Env2$Ps,out_HYBRID9$Anet/max(out_HYBRID9$Anet),
      ylim=c(0,1.2),xlim=c(0,-4),col=col_HYBRID9,lty=lty_HYBRID9,lwd=lwd_set,yaxs="i", xaxs="i")


################
# empirically derived soil moisture effect on ring width increment - this study
#Larch
swp=seq(0.00,-2.04,-0.02)
beta_swp_growth <- w_crit_optim_L # from this study's optimisation algorithm
lines(swp,btc_plot(swp,beta_swp_growth),
      ylim=c(0,1.2),xlim=c(0,-4),
      col=Growth_col,type="l",lwd=lwd_set,yaxs="i", xaxs="i")



#Spruce
beta_swp_growth <-w_crit_optim_S # from this study's optimisation algorithm
lines(swp,btc_plot(swp,beta_swp_growth),
      ylim=c(0,1.2),xlim=c(0,-4),
      col=Growth_col,type="l",lwd=lwd_set,yaxs="i", xaxs="i")


mtext(side=3, line =-0.82, adj= 0.015,"(a)",font=2)
mtext(side=3,line =-0.82,adj= 0.3, "Model response curves")

mtext("1-relative effect", side=2,outer=TRUE,line=1) #"normalised Source/Sink response"

axis(side=1, tck=.022, label=NA)
axis(side=1, tck=-.022, label=NA)
axis(side=2, mgp=c(3, .6, 0), tck=-.022, label=NA)
axis(side=2,labels=c(0.2,0.4,0.6,0.8,1),at=c(0.2,0.4,0.6,0.8,1),lwd=0,line=-0.8)
 

#####################################################################################
#####################################################################################
#####################################################################################
#####################################################################################

# observed sources and sinks


sigmoid <- function(params, x){params[1] / (1 + exp(-params[2] * (x - params[3])))}

##### Boyer et al 1970
#leaf water potential (bars)

x <- sunflower_enl1 <- c(-1.9,-2,-2.5,-3.1,-4,-6)*0.1 +0.60 # to MPa
y <- sunflower_enl2 <- c(32,31,20,3,0.0002,0.0001)  /32 # normalise

#fitmodel
fitmodel <- nls(y~a/(1 + exp(-b * (x-c))), start=list(a=1,b=10,c=0.33))
params   =  coef(fitmodel)
y2       <- sigmoid(params,seq(0,0.4,0.01))    
plot(seq(0,0.4,0.01)-0.60,y2,
     type="l",xlim=c(0,-4),ylim=c(0,1.2), xlab="",ylab="",xaxt="n",yaxt="n", yaxs="i", xaxs="i",lwd=lwd_set,col=Growth_col) # plot rescaled 
points(x-0.60,y,col=makeTransparent(Growth_col,alpha=90),cex=cex_set)

x <- sunflower_phot1<-c(-4,-8,-12,-16)*0.1 + 1.6 # to MPa
y <- sunflower_phot2<-c(32,32,22,10) /32

# fitmodel
fitmodel <- nls(y~(a/(1 + exp(-b * (x-c)))), start=list(a=1,b=2,c=0.2) )
params   =  coef(fitmodel)
y2       <- sigmoid(params,seq(0,1.6,0.01)  )  
lines((seq(0,1.6,0.01)-1.6),y2,
      type="l",xlim=c(0,-4),ylim=c(0,1),col=source_col,lwd=lwd_set) # plot recaled into -
points((x-1.60),y,col=makeTransparent(source_col,alpha=90),cex=cex_set)


##############Corn
#
x <- Corn_elong1 <- c(-2,-2,-2,-4,-4.4,-10.9,-12) * 0.1 + 1.2
y <- Corn_elong2 <- c(40,30,21,5,4,0,0)/40

#fitmodel
#fitmodel <- nls(y~(x^a),start=c(a=8.6))
#params   =  coef(fitmodel)
lines(seq(0,1,0.01)-1.2,seq(0,1,0.01)^8.6,
      xlim=c(0,-4),ylim=c(0,1),lwd=lwd_set,col=Growth_col) # plot recaled
points((x-1.2),y,pch=4,col= makeTransparent(Growth_col,alpha=90), cex=cex_set)

x <- corn_phot1<- c(-4,-5.6,-8.5,-11,-13,-16)* 0.1 +1.6
y <- corn_phot2<- c(55,50,45,40,27,10)/55

fitmodel <- nls(y ~(x^a), start=c(a=0.5))
params   =  coef(fitmodel)
lines(seq(0,1.2,0.01)-1.2,seq(0,1.2,0.01)^params,xlim=c(0,-4),ylim=c(0,1),col=source_col,lwd=lwd_set) # plot rescaled into -
points(x-1.2,y,pch=4,col= makeTransparent(source_col,alpha=90),cex=cex_set)


#Lempereuer et al 2016
#fig 5 (b)
GPP     <- c(6.1, 5, 4.5, 4.3, 4, 3.9, 3.2, 3, 2.2, 2.1, 2, 0.8, 0.8, 0.2, 0.1) / 6.1
swpLemp <- c(-0.5, -0.7, -1, -1.1, -1.4, -1.5, -1.6, -1.7, -2.4, -2.5, -2.6, -4, -4.2, -4.4, -4.7)
lines(swpLemp,GPP,
      type="l",col=source_col,lwd=lwd_set,xlim=c(0,-4)) # not fitted, just line plotted
lines(swpLemp,GPP,
      pch=5,type="p",col=makeTransparent(source_col,alpha=90),cex=cex_set,xlim=c(0,-4))#Lemp_col)
mtext(side=3, line =-0.82, adj= 0.015,"(b)",font=2)
mtext(side=3, line =-0.82, adj= 0.43, "Observed response curves")

swpLemp <-c(-0.5,-0.7,-1,-1.1)
deltaBA <- c(1.5,1,0.2,0) /1.5
lines(swpLemp,deltaBA,col=Growth_col)
lines(swpLemp,deltaBA,pch=5,type="p",col=makeTransparent(Growth_col,alpha=90),cex=cex_set)


##Lebourgeois et al 1998
#Pooled pine species see figure 5 bottom graphs.
Predlwp_pine_calp <- -c(0.4,0.45,0.49,0.95,1.1,1.35,1.5,2.15,2.5,2.65)
A_pine_calp <- c(2.7,3.4,3.2,2.9,1.3,0.8,0.5,0.5,0.3,0.5)

Predlwd_pine_corp <- -c(0.3,0.3,0.35,0.9,1.12,1.35,1.6,2.3,2.5)
A_pine_corp <- c(2.9,3.05,4,4,1.2,1,0.8,0.3,0.35)

x <- c(Predlwp_pine_calp,Predlwd_pine_corp) +2.6
y <- c(A_pine_calp,A_pine_corp) / max( c(A_pine_calp,A_pine_corp))
fitmodel <- nls(y~a/(1 + exp(-b * (x-c))), start=list(a=1,b=10,c=0.33) )
params   =  coef(fitmodel)
y2       <- sigmoid(params,seq(0,2.6,0.01))  

 
lines((seq(0,2.6,0.01)-2.6),y2,
      type="l",xlim=c(0,-4),ylim=c(0,1),col=source_col,lwd=lwd_set) # plot recaled into -
points((x-2.60),y,col=makeTransparent(source_col,alpha=90),cex=cex_set)

axis(side=1, tck=.022, label=NA)
axis(side=1,labels=c(0,-1,-2,-3,-4),at=c(0,-1,-2,-3,-4),lwd=0,line=-1)
axis(side=2, mgp=c(3, .6, 0), tck=-.022, label=NA)
axis(side=2,labels=c(0.2,0.4,0.6,0.8,1),at=c(0.2,0.4,0.6,0.8,1),lwd=0,line=-0.8)
 


################################################################
################################################################
################################################################
# axis title x-axis

mtext(expression(Psi[s]), side = 1, line = 0.8,adj = 0.3, outer = TRUE,cex=0.9)
mtext(paste('[MPa]'), side = 1, adj = 0.38, line = 0.8,outer =TRUE,cex=0.9)

# Legend placements
 par(fig = c(0.2,1, 0.001, 1),new = T,xpd=TRUE)
 par(mar=c(0,0,0,0),xpd=TRUE)
 

legend(-3.68,1.245,xpd=TRUE, legend= c('This study','HYBRID9','CLM4.5','CLM5',expression('SOX'[S]),expression('SOX'[L])),
       col=c(source_col,source_col,source_col,source_col,source_col,source_col),
       lty = c(Type1_lty,HYBRID_pch,CLM45_pch,CLM5_pch,SOX_pch,2),
       lwd=lwd_set,
       cex=0.8, bty="n", title="Origin",ncol=1)

legend(-3.79,0.65,xpd=TRUE,legend=c('Sources','Sinks'),fill=c(source_col,"brown"),col=c(source_col,"brown"),
       border= c(source_col,"brown"),cex=0.9)

#legend(-1,0.6,legend=c('Sources','Sinks'),legend=c('Sources','Sinks'))


legend(-3.89,0.385,xpd=TRUE,legend = c("sunflower","maize","oak",'pine'),
       pch  = c(1,4,5),
       col  = makeTransparent(c(1,1,1)),
       bty  = "n",cex=0.8, ncol=1,
      title= "Organisms")

legend(-1.1,0.59, legend = "Photosynthesis",fill=c(source_col),bty="n", cex=0.8,border=NA, x.intersp = 0.3)
legend(-1.1,1.187, legend = "Photosynthesis",fill=c(source_col),bty="n", cex=0.8,border=NA, x.intersp = 0.3)
legend(-1.9,1.03, legend = "Stomata",fill=c(source_col),bty="n", cex=0.8,border=NA, x.intersp = 0.25)


################################################################
################################################################
################################################################
#STOMATAL CONDUCTANCE RESPONSE TO SWP - Figure box C)
 par(fig = c(0.39,0.624, 0.551, 0.82),new = T,xpd=TRUE)
 par(mar=c(0,0,0,0))
#################
# SOX MODEL

# normalise output
PA_norm <- max(outPA$Gc)
LD_norm <- max(outLD$Gc)

# Ps = soil water potential
plot(outPA$Gc/PA_norm ~ Env2$Ps, 
     type = "l",xlim=c(0,-4),ylim=c(0,1.2),col=col_SOX, xlab="",ylab="",xaxt="n",yaxt="n", yaxs="i", xaxs="i",
     lty=lty_SOX,lwd=lwd_set)

lines(outLD$Gc/LD_norm ~ Env2$Ps,
      type = "l",col=col_SOX, xlab="",ylab="",xaxt="n",yaxt="n", yaxs="i", xaxs="i", lty=2,lwd=lwd_set)


#################
# CLM4.5 soil moisture effect on stomatal conductance

# Ps = soil water potential
CLM_norm <- max(out_CLM$NET_Temperate$Gc)
#stomatal conductance- soil water response
lines(out_CLM$NET_Temperate$Gc/CLM_norm ~ Env2$Ps, 
     type = "l",xlim = c(0,-4), ylim = c(0,1.2), col = col_CLM45, xlab="swp", xaxt="n", yaxt="n", yaxs="i",
     xaxs = "i",lty = lty_CLM45, lwd = lwd_set)

##################
# HYBRID9
# this is the critical threshold in HYBRID9 for stomatal openness
H9_beta_stoma=-1e+05 * 9.81 / 1.0E6 # convert to MPa for plotting
lines(btc_plot(Env2$Ps,H9_beta_stoma)/max(btc_plot(Env2$Ps,H9_beta_stoma))~Env2$Ps,
      xlim=c(0,-5),ylim=c(0,1.2),type="l",lty=lty_HYBRID9,lwd=lwd_set,col=col_HYBRID9)


################
# CLM5, Kennedy et al 2019
# in Kennedy et al 2019: "Fisher et al., 2006, report average 𝜓le of −1.71 MPa during the wet season and −2.47 MPa during the dry season"
# parameters from Table 1 Kennedy et al 2019
ck = 2.95   # shape fitting parameter
p50 = -1.75 #(MPa) # potential at 50% loss of stomatal conductance

psileaf <- seq(0,-5,-0.3) # note: unnaturally low for swp leaf, in the lower values
# eq 9, sigmoidal function: 
lines(psileaf,2^(-(psileaf/p50)^ck),col=col_CLM5,lty=lty_CLM5,lwd=lwd_set)

mtext(side=3, line =-0.7, adj= 0.08005,"(c)",font=2)
#mtext(side=3, line =-0.85, adj= 0.8505,"Stomata")


axis(side=1, tck=.022, label=NA)
axis(side=1, tck=-.022, label=NA)
axis(side=2, mgp=c(3, .6, 0), tck=-.022, label=NA)
axis(side=1,labels=c(-1,-2,-3,-4),at=c(-1,-2,-3,-4),lwd=0,line=-1.32)
axis(side=2,labels=c(0.2,0.6,1),at=c(0.2,0.6,1),lwd=0,line=-0.94)

          
dev.off()

```


#Supplementaries only

```{r plot_map}
# plot study site location within Europe

jpeg('figures/APPDX_LocationN08.jpeg',units = "cm",width = 18, height=12,res=300)
# site coordinates
Site<-data.frame(SITE=c("N08"),lat =c(46.3026),lon=c(7.7411)) 

# prepare map dimensions for plotting 
lat.low=44
lat.high=49

lon.low=5
lon.high=18

# plot map
newmap<-getMap(resolution ="low")
maps::map(newmap,ylim=c(lat.low,lat.high),xlim=c(lon.low,lon.high),fill = TRUE,col="grey")
map.axes() 

abline(h=seq(lat.low,lat.high,0.5),col="light grey")
abline(v=seq(lon.low,lon.high,0.5),col="light grey")
lines(Site$lon,Site$lat, type="p",col="red",cex=0.8)

mtext(side=2,"Latitude", line=3)
mtext(side=1,"Longitude",line=3)
mtext(side=1,"Italy",    line = -5,cex=0.8)
mtext(side=1,"France",   line = -10,adj=0.02,cex=0.8)
mtext(side=3,"Germany",  line=-2,cex=0.8)
mtext(side=3,"Austria",  line=-4,cex=0.8,adj=0.7)
mtext(side=3,"Slovenia", line=-8.8,cex=0.8,adj=0.75)
mtext(side=3,"Switzerland", line=-6.2,cex=0.8,adj=0.2)

dev.off()
```

```{r NORUNdata_soil_water_retenetion_curves, eval=FALSE}
# chung rearrangse HYBRID9 hydroogical output for storage and quick retrieval for soil water retentioncurve plotting
    

#Soil–Water Retention Curves Derived as a Function of Soil Dry Density(Chen 2018)
  
 vvol <- c(0.44,0.43,0.2,0.12,0.1,0.09,0.06)
 swc<-  c(0.1,1,5,10,50,100,200)
 
 vvol2 <- c(0.41,0.41,0.28,0.18,0.12,0.11,0.08)
 plot(-(swc/1000),vvol)#  ,xlim=c(0,-1),ylim=c(0,0.4))
 
# comparison with modelled soil-water retention curves
collect_layers <- list() 
mdd.swp <-  data.frame(swp=rep(NA,length(diag_daily$vvol4)),vvol=rep(NA,length(diag_daily$vvol4)))

swp_subdaily_in <- swp_subdaily
#1
mdd.swp$swp  <- swp.midday[swp.midday$year< 2016,"swp1"]
mdd.swp$vvol <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],12,sub=TRUE,start=1901,end=2015)[,"b1"]
collect_layers[[1]] <- mdd.swp

#2
mdd.swp$swp <-  swp.midday[swp.midday$year< 2016,"swp2"]
mdd.swp$vvol <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],12,sub=TRUE,start=1901,end=2015)[,"b2"]
collect_layers[[2]] <-mdd.swp

#3
mdd.swp$swp <-  swp.midday[swp.midday$year< 2016,"swp3"]
mdd.swp$vvol <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],12,sub=TRUE,start=1901,end=2015)[,"b3"]
collect_layers[[3]] <-mdd.swp
#4
mdd.swp$swp <- diag_daily$swp4 
mdd.swp$vvol <- diag_daily$vvol4
collect_layers[[4]] <-mdd.swp
#5
mdd.swp$swp <- swp.midday[swp.midday$year< 2016,"swp5"]
mdd.swp$vvol <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],12,sub=TRUE,start=1901,end=2015)[,"b5"]
collect_layers[[5]] <-mdd.swp
#6
mdd.swp$swp <- swp.midday[swp.midday$year< 2016,"swp6"]
mdd.swp$vvol <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],12,sub=TRUE,start=1901,end=2015)[,"b6"]
collect_layers[[6]] <-mdd.swp
#7
mdd.swp$swp <- swp.midday[swp.midday$year< 2016,"swp7"]
mdd.swp$vvol <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],12,sub=TRUE,start=1901,end=2015)[,"b7"]
collect_layers[[7]] <-mdd.swp
#8
mdd.swp$swp <- swp.midday[swp.midday$year< 2016,"swp8"]
mdd.swp$vvol <- create_Timestamp(sm_obs_cent[sm_obs_cent$jyear<2016,],12,sub=TRUE,start=1901,end=2015)[,"b8"]
collect_layers[[8]] <-mdd.swp
 swp_vvol_midday_all_layers <- collect_layers
 save(swp_vvol_midday_all_layers,file='data/HYBRID9_SWP_VVOL_all_layers_APPDX.RData')
 
```

```{r plot_root_dist}
jpeg(filename='figures/APPDX_rooting_distribution.jpeg')

# turn parametersiation into a loop:
for (i in 1:3){
if(i==1){ #tropical parameters
a_root  = 0.15 / 0.91 #! Ratio of z50 to z95 (m/m)
b_root  = 0.91 / 35.0 #! Ratio of z95 to height (m/m)
}

if(i==2){ # temperate parameters}
a_root  = 0.23 / 1.21 #! Ratio of z50 to z95 (m/m)
b_root  = 1.21 / 35.0 #! Ratio of z95 to height (m/m)
}
  
  if(i==3){ # temperate parameters, smaller trees}, used in this publication
a_root  = 0.23#0.14/ 0.76 #0.23 / 1.21 #! Ratio of z50 to z95 (m/m)
b_root  = 0.05496#0.76 /22 #1.21 / 22.0 #! Ratio of z95 to height (m/m)
  }
  #if the same ratio is supposed to be maintained, but the height changing, the above parameters must be recalculated


if(i==1){
  z95 = 1000 * b_root* 34.49 # *1000 to change to mm, to correspond to soildepths.
z50 = a_root * z95
c_root =- 1.27875 / (log10(z50) - log10(z95))  
}
if(i==2 |i==3){
  if(i==3 ){
     z95 = 1000 * b_root * 22#34.49 # *1000 to change to mm, to correspond to soildepths.
  }else{
      z95 = 1000 * b_root * 34.49
  }
z50 = a_root * z95
c_root = -1.88 / (log10(z50) - log10(z95))  # -1.88
print(paste(i,c_root,sep=" "))
}
#preparing for column-plot, to plot water content in middle of each soil layer
layer<-c(-45,-91,-166,-289,-493,-829,-1383,-2296)
dz<-vector()# soil layer thickness
zc<- vector() # depth of that soil layer's middle
layer_all<-c(0,layer) # soil layer boundaries, include upper boundary
for (KL in 2: length(layer_all)){
dz[KL-1]<-layer_all[KL]-layer_all[KL-1]
}
for (KL in 1: length(layer_all)){
zc[KL]<- -(layer_all[KL]-(-dz[KL])/2)
}

root_col<-vector()
for (KL in 1: length(layer_all)){
root_col[KL] = -dz[KL] * (c_root / z50) * 
                            (zc[KL] / z50) ^ (c_root - 1.0) * 
                            (1 + (zc[KL] / z50)^(c_root)) ^
                              (-2.0)
}
if(i==1){
root_col_tropics <- root_col
}
if(i==2){
root_col_temperate <- root_col
}
if(i==3){
root_col_temperate_22m <- root_col
}

}



#plot(root_distr_l,layer,type="l")
#  plot(root_col_temperate_22m[1:8],layer,col="red",type="l",ylab="soil depth [mm]",xlim=c(0,0.4),ylim=c(-3000,0), xlab= "Relative amount of roots [%]")#,
#       main="Rooting distribution Comparison Schenk et Jackson (2002)" )

#lines(root_col_temperate[1:8],layer,col="black")
#lines(root_col_temperate_22m[1:8],layer, col="orange")
#plot(cumsum(root_distr_l),layer,type="l",xaxt="n")
#axis(1, at=c(0.05,0.1,0.15,0.2), labels=c(0.05,0.1,0.15,0.4)*100)

    
#fit<-lm(layer ~ (poly (root_distr,2,raw=TRUE)) )
#lines(fit,layer)

# temperate parameters, smaller trees- used to basically recover the temperature forest rooting distribution at a tree height of 35 m.
a_root  = 0.23#0.14/ 0.76 #0.23 / 1.21 #! Ratio of z50 to z95 (m/m)
b_root  = 0.05496

z95 = 1.0E3 * b_root * 22.65347 # scaled by height 
z50 = a_root * z95

c_root = 1.88  / (log10(z50) - log10(z95))

rootr_col<-vector()
#rootr_col[1] = 1/ (1 + (-dz[1] / z50) ** c_root)
#root_col[1] = 1/ (1 + (-dz[1] / z50) ** c_root)
          for( KL in 1: length(layer)){
   rootr_col[KL] =  1 / (1 + (-dz[KL]   / z50) ** c_root) 
                            - 1 / (1 + (-dz[KL-1] / z50) ** c_root)
   
   root_col[KL] = dz[KL] * (c_root / z50) * 
                            (zc[KL] / z50) ^ (c_root - 1.0) * 
                            (1 + (zc[KL] / z50)^(c_root)) ^
                              (-2.0)
          }


root_col<- c(    4.2859907E-03,  2.8834227E-02,  0.1343957,      0.3399178,      0.3249828,
  0.1258425 ,     3.2185972E-02,  9.5550418E-03)# HYBRID output, with parameters:  aroot= 0.23 broot=, 0.05496,  and 1.88 in c_root equation
par(mfrow=c(1,2))
# layer/10 <- from mm to cm
 plot(root_col[1:8]*100,layer/10,col="dark orange",type="l",ylab="Soil depth [cm]",xlim=c(0,40),ylim=c(-300,0), xlab= "Relative amount of roots [%]")
 abline(h=-10, col="dark grey", lty=3)
 abline(h=-70, col="dark grey", lty=3)
 mtext("A)",side=3,line=-1,adj= -0.4)
 lines(root_col_temperate[1:8]*100,layer/10,col= "black")
 legend("bottomright",legend=c("temperate forest", "this site"),fill=c("black","dark orange"),bty="n")
  
 barplot(rev(root_col[1:8]*100),xlim=c(0,40),horiz=TRUE,names.arg = rev(c("L1","L2","L3","L4","L5","L6","L7","L8")),cex.names = 0.7,las=1,ylab="Soil layer", xlab= "Relative amount of roots [%]")
 mtext("B)",side=3,line=-1,adj= -0.4)

dev.off()
```

```{r plot_soil_water_retention_curves}

load(file='data/APPDX_HYBRID9_SWP_VVOL_all_layers.RData')



########
collect_layers <- swp_vvol_midday_all_layers

jpeg(filename='figures/APPENDIX_soil_water_retention_curves_all_with_data.jpeg',units="in", width=10, height=10, res=300)
l=1
  mdd.swp<-collect_layers[[l]]
  plot(mdd.swp$swp,mdd.swp$vvol,xlim=c(0,-1),ylim=c(0,0.4),col=makeTransparent(I(brewer.pal(8, name = 'Dark2'))[l],alpha = 30),xlab="swp[MPa]",ylab="vvol [mm3/mm3]",lwd=0.5)

for(l in 2:8){
  mdd.swp<-collect_layers[[l]]
  lines(mdd.swp$swp,mdd.swp$vvol,xlim=c(0,-1),ylim=c(0,0.4),col=makeTransparent(I(brewer.pal(8, name = 'Dark2'))[l],alpha = 30),type="p",lwd=0.5)
}
 lines(-(swc/1000),vvol, col="magenta",lwd=1.5,xlim=c(0,-1),ylim=c(0,0.4))
 lines(-(swc/1000),vvol2,xlim=c(0,-1),lwd=1.5,ylim=c(0,0.4))
 

MyLab <- c("layer 1","layer 2","layer 3", "layer 4","layer 5","layer 6","layer 7","layer 8","Chen, 2018 1.3 g/cm3", "Chen, 2018 1.4 g/cm3")
MyCol <- c(I(brewer.pal(8, name = 'Dark2'))[1:8],"magenta","black")

MyOrder = matrix(1:10,nrow = 5, ncol = 2)
legend("topright",MyLab[MyOrder], fill = MyCol[MyOrder],ncol=2)
 
dev.off()
############

```



# References

To be added
