Newer
Older
##################################### FUNCTION TO COMPUTE CLIMATIC DATA FOR FRANCE
## function to compute sum of degree days above 5.56
##' .. Compute teh sum of degree days above 5.56C from monthly data with a spline ..
##'
##' .. A spline is used to smooth the variation of T over each day from mothly data. Require package 'season' ..
##' @title
##' @param temp monthly temperature in C
##' @param threshold.sgdd threshold of temperature to compte sum of degree day default 5.56
##' @return
##' @author Kunstler
fun.sgdd <- function(temp, threshold.sgdd = 5.56) {
require(season)
temp <- unlist(temp)
ndays.month <- flagleap(data.frame(year = 2013, month = 1:12), F)$ndaysmonth
ndays.year <- sum(ndays.month)
x <- c(-ndays.month[12]/2, (ndays.month/2) + c(0, cumsum(ndays.month[1:11])),
ndays.year + ndays.month[1]/2)
## plot(x,c(temp[12],temp,temp[1]),xlim=c(0,365),type='b')
myfit <- loess(y ~ x, data = data.frame(y = c(temp[12], temp, temp[1]), x = x),
span = 0.4, degree = 2)
mypred <- predict(myfit, 1:ndays.year)
## lines( 1:ndays.year,mypred,col='red')
sgdd <- sum(mypred[mypred >= threshold.sgdd] - threshold.sgdd)
return(sgdd)
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
################################################################################################# function to compute max soil water content based on Piedallu 2010 Geoderma
fun.swhc <- function(affroc, cailloux, text2, text1, prof2, prof1, codeprof, codecaillou,
data.texture) {
## browser()
### transform code of prof= soil depth in cm and code of percentage of rock in %
prof2.t <- prof2
prof2.t[is.na(prof2)] <- 0
prof2.cm.temp <- codeprof[prof2 + 1]
prof2.cm.temp[is.na(prof2)] <- 0
prof1.t <- prof1
prof1.t[is.na(prof1)] <- 0
text1[is.na(text1)] <- 0
prof1.cm <- codeprof[prof1.t + 1]
prof1.cm[is.na(prof1)] <- 0
cailloux.t <- cailloux
cailloux.t[is.na(cailloux)] <- 0
cailloux.perc <- codecaillou[cailloux.t + 1]
cailloux.perc[is.na(cailloux)] <- 0
## compute depth of second horizon
prof2.cm <- prof2.cm.temp - prof1.cm
perc.top1 <- apply(cbind(10/prof1.cm, rep(1, length(prof1.cm))), MARGIN = 1,
FUN = min)
perc.top2 <- apply(cbind(apply(cbind(rep(0, length(prof1.cm)), (10 - prof1.cm)/prof2.cm.temp),
MARGIN = 1, FUN = max), rep(1, length(prof1.cm))), MARGIN = 1, FUN = min)
caill <- (1 - codecaillou[affroc + 1]/100) * (1 - (sqrt(cailloux.perc/100))^3)
## print(c(prof2.cm,perc.top2))
top.hor <- (prof1.cm * perc.top1 * data.texture[text1 + 1, "U_almajou_top"] +
prof2.cm * perc.top2 * data.texture[text2 + 1, "U_almajou_top"])
sub.hor <- (prof1.cm * (1 - perc.top1) * data.texture[text1 + 1, "U_almajou_sub"] +
prof2.cm * (1 - perc.top2) * data.texture[text2 + 1, "U_almajou_sub"])
swhc <- caill * (top.hor + sub.hor)
return(swhc)
################################################################# function to compute PET
fun.PET <- function(i, rad, temp) {
## conversion data
require(season)
ndays.month <- flagleap(data.frame(year = 2013, month = 1:12), F)$ndaysmonth
Rs <- unlist(rad[i, ])/100/30 * 1000
Ta <- unlist(temp[i, ])
# compute TURC FORMULA
vec.pet <- 0.4 * ((0.0239001 * Rs) + 50) * (Ta/(Ta + 15))/30 * ndays.month
## with Rs daily global (total) solar radiation (kJ/m2/day) and Ta monthly mean
## temp in C and ndays.month number of day in each month
return(vec.pet)
}
#################################################################################### FUNCTION TO COMPUTE WATER BUDGET based on Bugmann & Cramer 1998
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
fun.WaterBudget <- function(i, prcp.m, PET.m, Ta.m, SWHC.v, n = 2) {
WATER <- rep(NA, 12 * n)
WATER[1] <- unlist(SWHC.v)[i]
SWHC <- unlist(SWHC.v)[i]
AET <- rep(NA, 12 * n)
D <- rep(NA, 12 * n)
prcp <- rep(unlist(prcp.m[i, ]), n)
PET <- rep(unlist(PET.m[i, ]), n)
Ta <- rep(unlist(Ta.m[i, ]), n)
for (i in 2:(12 * n)) {
Pi <- min(0.3 * prcp[i], PET[i])
Ps <- prcp[i] - Pi
S <- 120 * WATER[i - 1]/SWHC
D[i] <- PET[i] - Pi
AET[i] <- max(0, min(D[i], S))
WATER[i] <- max(0, min(WATER[i - 1] + Ps - AET[i], SWHC))
}
## plot(1:(12*n),WATER,type='b')
WATER <- WATER[13:24]
Ta <- Ta[13:24]
AET <- AET[13:24]
D <- D[13:24]
## NEED to compute mean water availability of the growing season
WB.s <- mean(WATER[Ta > 5.56], na.rm = T)
## NEED to compute mean water availability of the year
WB.y <- mean(WATER, na.rm = T)
## NEED to compute water stress index of the growing season
WS.s <- sum(na.omit(AET[Ta > 5.56 & D > 0]))/sum(na.omit(D[Ta > 5.56 & D > 0]))
## NEED to compute water stress index of the year
WS.y <- sum(na.omit(AET[D > 0]))/sum(na.omit(D[D > 0]))
return(c(WB.s = WB.s, WB.y = WB.y, WS.s = WS.s, WS.y = WS.y))
}