diff --git a/DESCRIPTION b/DESCRIPTION
index 69e42e13db4587f8cec92cdae2e1db58f22b86f0..b8f7a63bef9c97263fa79de68d103a896b6acd54 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.0.5.15
+Version: 1.0.5.16
Date: 2017-01-23
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl")),
diff --git a/R/DataAltiExtrapolation_Valery.R b/R/DataAltiExtrapolation_Valery.R
index fce85cb8128820e9bb903f6fab695728a26ff5f5..397b2de5f4c06d27a18bf6c9799dedd9db1ab54e 100644
--- a/R/DataAltiExtrapolation_Valery.R
+++ b/R/DataAltiExtrapolation_Valery.R
@@ -1,14 +1,17 @@
-DataAltiExtrapolation_Valery <- function(DatesR, Precip, PrecipScale = TRUE, TempMean, TempMin = NULL, TempMax = NULL, ZInputs, HypsoData, NLayers, verbose = TRUE){
+DataAltiExtrapolation_Valery <- function(DatesR,
+ Precip, PrecipScale = TRUE,
+ TempMean, TempMin = NULL, TempMax = NULL,
+ ZInputs, HypsoData, NLayers,
+ verbose = TRUE) {
-
- ##Altitudinal_gradient_functions_______________________________________________________________
- ##unique_gradient_for_precipitation
- GradP_Valery2010 <- function(){
- return(0.00041); ### value from Val? PhD thesis page 126
- }
- ##daily_gradients_for_mean_min_and_max_air_temperature
- GradT_Valery2010 <- function(){
- RESULT <- matrix(c(
+ ##Altitudinal_gradient_functions_______________________________________________________________
+ ##unique_gradient_for_precipitation
+ GradP_Valery2010 <- function() {
+ return(0.00041) ### value from Valerie PhD thesis page 126
+ }
+ ##daily_gradients_for_mean_min_and_max_air_temperature
+ GradT_Valery2010 <- function() {
+ RESULT <- matrix(c(
01, 01, 0.434, 0.366, 0.498,
02, 01, 0.434, 0.366, 0.500,
03, 01, 0.435, 0.367, 0.501,
@@ -374,130 +377,170 @@ DataAltiExtrapolation_Valery <- function(DatesR, Precip, PrecipScale = TRUE, Tem
28, 12, 0.431, 0.366, 0.495,
29, 12, 0.431, 0.366, 0.495,
30, 12, 0.432, 0.366, 0.496,
- 31, 12, 0.433, 0.366, 0.497),ncol=5,byrow=TRUE);
- dimnames(RESULT) <- list(1:366,c("day","month","grad_Tmean","grad_Tmin","grad_Tmax"));
- return(RESULT);
- }
+ 31, 12, 0.433, 0.366, 0.497), ncol = 5, byrow = TRUE)
+
+ dimnames(RESULT) <- list(1:366, c("day", "month", "grad_Tmean", "grad_Tmin", "grad_Tmax"))
+
+ return(RESULT)
+
+ }
##Format_______________________________________________________________________________________
- HypsoData <- as.double(HypsoData);
- ZInputs <- as.double(ZInputs);
+ HypsoData <- as.double(HypsoData)
+ ZInputs <- as.double(ZInputs)
+
##ElevationLayers_Creation_____________________________________________________________________
- ZLayers <- as.double(rep(NA,NLayers));
- if(!identical(HypsoData,as.double(rep(NA,101)))){
- nmoy <- 100 %/% NLayers;
- nreste <- 100 %% NLayers;
- ncont <- 0;
- for(iLayer in 1:NLayers){
- if(nreste > 0){ nn <- nmoy+1; nreste <- nreste-1; } else { nn <- nmoy; }
- if(nn==1){ ZLayers[iLayer] <- HypsoData[ncont+1]; }
- if(nn==2){ ZLayers[iLayer] <- 0.5 * (HypsoData[ncont+1] + HypsoData[ncont+2]); }
- if(nn>2 ){ ZLayers[iLayer] <- HypsoData[ncont+nn/2]; }
- ncont <- ncont+nn;
+ ZLayers <- as.double(rep(NA, NLayers))
+
+ if (!identical(HypsoData, as.double(rep(NA, 101)))) {
+ nmoy <- 100 %/% NLayers
+ nreste <- 100 %% NLayers
+ ncont <- 0
+
+ for (iLayer in 1:NLayers) {
+ if (nreste > 0) {
+ nn <- nmoy + 1
+ nreste <- nreste - 1
+ } else {
+ nn <- nmoy
+ }
+ if (nn == 1) {
+ ZLayers[iLayer] <- HypsoData[ncont + 1]
+ }
+ if (nn == 2) {
+ ZLayers[iLayer] <- 0.5 * (HypsoData[ncont + 1] + HypsoData[ncont + 2])
+ }
+ if (nn > 2) {
+ ZLayers[iLayer] <- HypsoData[ncont + nn / 2]
+ }
+ ncont <- ncont + nn
+ }
}
- }
##Precipitation_extrapolation__________________________________________________________________
##Initialisation
- if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
- LayerPrecip <- list(as.double(Precip))
- } else {
- ##Elevation_gradients_for_daily_mean_precipitation
- GradP <- GradP_Valery2010() ### single value
- TabGradP <- rep(GradP, length(Precip))
- ##Extrapolation
- ##Thresold_of_inputs_median_elevation
- Zthreshold <- 4000
- LayerPrecip_df <- sapply(1:NLayers, function(iLayer) {
- ##If_layer_elevation_smaller_than_Zthreshold
- if (ZLayers[iLayer] <= Zthreshold) {
- prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs)))
- ##If_layer_elevation_greater_than_Zthreshold
- } else {
- ##If_inputs_median_elevation_smaller_than_Zthreshold
- if (ZInputs <= Zthreshold) {
- prcp <- as.double(Precip * exp(TabGradP*(Zthreshold - ZInputs)))
- ##If_inputs_median_elevation_greater_then_Zthreshold
+ if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
+ LayerPrecip <- list(as.double(Precip))
+ } else {
+ ##Elevation_gradients_for_daily_mean_precipitation
+ GradP <- GradP_Valery2010() ### single value
+ TabGradP <- rep(GradP, length(Precip))
+ ##Extrapolation
+ ##Thresold_of_inputs_median_elevation
+ Zthreshold <- 4000
+ LayerPrecip_df <- sapply(1:NLayers, function(iLayer) {
+ ##If_layer_elevation_smaller_than_Zthreshold
+ if (ZLayers[iLayer] <= Zthreshold) {
+ prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs)))
+ ##If_layer_elevation_greater_than_Zthreshold
} else {
- prcp <- as.double(Precip)
+ ##If_inputs_median_elevation_smaller_than_Zthreshold
+ if (ZInputs <= Zthreshold) {
+ prcp <- as.double(Precip * exp(TabGradP * (Zthreshold - ZInputs)))
+ ##If_inputs_median_elevation_greater_then_Zthreshold
+ } else {
+ prcp <- as.double(Precip)
+ }
}
+ return(prcp)
+ })
+ if (PrecipScale) {
+ LayerPrecip_df <- LayerPrecip_df / rowMeans(LayerPrecip_df) * Precip
+ LayerPrecip_df[is.nan(LayerPrecip_df)] <- 0
}
- return(prcp)
- }
- )
- if (PrecipScale) {
- LayerPrecip_df <- LayerPrecip_df / rowMeans(LayerPrecip_df) * Precip
- LayerPrecip_df[is.nan(LayerPrecip_df)] <- 0
+ LayerPrecip <-
+ lapply(seq_len(ncol(LayerPrecip_df)), function(x) LayerPrecip_df[, x])
}
- LayerPrecip <- lapply(seq_len(ncol(LayerPrecip_df)), function(x) LayerPrecip_df[, x])
- }
##Temperature_extrapolation____________________________________________________________________
##Initialisation
- LayerTempMean <- list(); LayerTempMin <- list(); LayerTempMax <- list();
- if(identical(ZInputs,HypsoData[51]) & NLayers==1) {
- LayerTempMean[[1]] <- as.double(TempMean);
- if(!is.null(TempMin) & !is.null(TempMax)){ LayerTempMin[[1]] <- as.double(TempMin); LayerTempMax[[1]] <- as.double(TempMax); }
- } else {
- ##Elevation_gradients_for_daily_mean_min_and_max_temperature
- GradT <- as.data.frame(GradT_Valery2010())
- iday <- match(format(DatesR,format="%d%m"), sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
- TabGradT <- GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
- ##Extrapolation
- ##On_each_elevation_layer...
- for(iLayer in 1:NLayers){
- LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs-ZLayers[iLayer])*abs(TabGradT[, "grad_Tmean"])/100);
- if(!is.null(TempMin) & !is.null(TempMax)){
- LayerTempMin[[iLayer]] <- as.double(TempMin + (ZInputs-ZLayers[iLayer])*abs(TabGradT[, "grad_Tmin"])/100);
- LayerTempMax[[iLayer]] <- as.double(TempMax + (ZInputs-ZLayers[iLayer])*abs(TabGradT[, "grad_Tmax"])/100);
+ LayerTempMean <-
+ list()
+ LayerTempMin <- list()
+ LayerTempMax <- list()
+
+ if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
+ LayerTempMean[[1]] <- as.double(TempMean)
+
+ if (!is.null(TempMin) & !is.null(TempMax)) {
+ LayerTempMin[[1]] <- as.double(TempMin)
+ LayerTempMax[[1]] <- as.double(TempMax)
+ }
+ } else {
+ ##Elevation_gradients_for_daily_mean_min_and_max_temperature
+ GradT <- as.data.frame(GradT_Valery2010())
+ iday <- match(format(DatesR, format = "%d%m"),
+ sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
+ TabGradT <-
+ GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
+ ##Extrapolation
+ ##On_each_elevation_layer...
+ for (iLayer in 1:NLayers) {
+ LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmean"]) / 100)
+ if (!is.null(TempMin) & !is.null(TempMax)) {
+ LayerTempMin[[iLayer]] <- as.double(TempMin + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmin"]) / 100)
+ LayerTempMax[[iLayer]] <- as.double(TempMax + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmax"]) / 100)
+
+ }
}
}
- }
-
+
##Solid_Fraction_for_each_elevation_layer______________________________________________________
- LayerFracSolidPrecip <- list();
- ##Thresold_of_inputs_median_elevation
- Zthreshold <- 1500;
- ##On_each_elevation_layer...
- for(iLayer in 1:NLayers){
- Option <- "USACE";
- if(!is.na(ZInputs)){ if(ZInputs < Zthreshold & !is.null(TempMin) & !is.null(TempMax)){ Option <- "Hydrotel"; } }
- ##Turcotte_formula_from_Hydrotel
- if(Option=="Hydrotel"){
- TempMin <- LayerTempMin[[iLayer]];
- TempMax <- LayerTempMax[[iLayer]];
- SolidFraction <- 1 - TempMax/(TempMax - TempMin);
- SolidFraction[TempMin >= 0] <- 0;
- SolidFraction[TempMax <= 0] <- 1;
- }
- ##USACE_formula
- if(Option=="USACE"){
- USACE_Tmin <- -1.0;
- USACE_Tmax <- 3.0;
- TempMean <- LayerTempMean[[iLayer]];
- SolidFraction <- 1- (TempMean - USACE_Tmin)/(USACE_Tmax - USACE_Tmin);
- SolidFraction[TempMean > USACE_Tmax] <- 0;
- SolidFraction[TempMean < USACE_Tmin] <- 1;
+ LayerFracSolidPrecip <- list()
+
+ ##Thresold_of_inputs_median_elevation
+ Zthreshold <- 1500
+
+ ##On_each_elevation_layer...
+ for (iLayer in 1:NLayers) {
+ Option <- "USACE"
+
+ if (!is.na(ZInputs)) {
+ if (ZInputs < Zthreshold & !is.null(TempMin) & !is.null(TempMax)) {
+ Option <- "Hydrotel"
+ }
+ }
+ ##Turcotte_formula_from_Hydrotel
+ if (Option == "Hydrotel") {
+ TempMin <- LayerTempMin[[iLayer]]
+ TempMax <- LayerTempMax[[iLayer]]
+ SolidFraction <- 1 - TempMax / (TempMax - TempMin)
+ SolidFraction[TempMin >= 0] <- 0
+ SolidFraction[TempMax <= 0] <- 1
+ }
+ ##USACE_formula
+ if (Option == "USACE") {
+ USACE_Tmin <- -1.0
+ USACE_Tmax <- 3.0
+ TempMean <- LayerTempMean[[iLayer]]
+ SolidFraction <- 1 - (TempMean - USACE_Tmin) / (USACE_Tmax - USACE_Tmin)
+ SolidFraction[TempMean > USACE_Tmax] <- 0
+ SolidFraction[TempMean < USACE_Tmin] <- 1
+ }
+ LayerFracSolidPrecip[[iLayer]] <- as.double(SolidFraction)
}
- LayerFracSolidPrecip[[iLayer]] <- as.double(SolidFraction);
- }
##END__________________________________________________________________________________________
- return(list(LayerPrecip=LayerPrecip,LayerTempMean=LayerTempMean,LayerTempMin=LayerTempMin,LayerTempMax=LayerTempMax,
- LayerFracSolidPrecip=LayerFracSolidPrecip,ZLayers=ZLayers));
-
+ return(list(LayerPrecip = LayerPrecip,
+ LayerTempMean = LayerTempMean,
+ LayerTempMin = LayerTempMin,
+ LayerTempMax = LayerTempMax,
+ LayerFracSolidPrecip = LayerFracSolidPrecip,
+ ZLayers = ZLayers))
+
+
}