##################################### 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)
}



################################################################################################# 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

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))
}