Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 1576 additions and 1123 deletions
+1576 -1123
R/ErrorCrit.R
View file @ 4bccfee2
ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbose = TRUE) {
ErrorCrit <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## ---------- Arguments check
if (!inherits(InputsCrit, "InputsCrit")) {
stop("InputsCrit must be of class 'InputsCrit'")
}
}
if (!inherits(OutputsModel, "OutputsModel")) {
stop("OutputsModel must be of class 'OutputsModel'")
}
if (!missing(FUN_CRIT)) {
warning("deprecated 'FUN_CRIT' argument. The error criterion function is now automatically get from the 'InputsCrit' object", call. = FALSE)
}
## ---------- Criterion computation
## ----- Single criterion
if (inherits(InputsCrit, "Single")) {
FUN_CRIT <- match.fun(InputsCrit$FUN_CRIT)
......@@ -25,10 +21,10 @@ ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbo
warnings = warnings,
verbose = verbose)
}
## ----- Multiple criteria or Composite criterion
if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
listOutputsCrit <- lapply(InputsCrit, FUN = function(iInputsCrit) {
FUN_CRIT <- match.fun(iInputsCrit$FUN_CRIT)
......@@ -37,12 +33,12 @@ ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbo
warnings = warnings,
verbose = verbose)
})
listValCrit <- sapply(listOutputsCrit, function(x) x[["CritValue"]])
listNameCrit <- sapply(listOutputsCrit, function(x) x[["CritName"]])
listweights <- unlist(lapply(InputsCrit, function(x) x[["Weights"]]))
listweights <- listweights / sum(listweights)
listweights <- listweights / sum(listweights)
if (inherits(InputsCrit, "Compo")) {
CritValue <- sum(listValCrit * listweights)
OutputsCritCompo <- list(MultiCritValues = listValCrit,
......@@ -61,7 +57,7 @@ ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbo
message("Crit. Composite = ", sprintf("%.4f", CritValue))
msgForm <- paste(sprintf("%.2f", listweights), listNameCrit, sep = " * ", collapse = ", ")
msgForm <- unlist(strsplit(msgForm, split = ","))
msgFormSep <- rep(c(",", ",", ",\n\t\t "), times = ceiling(length(msgForm)/3))[1: length(msgForm)]
msgFormSep <- rep(c(",", ",", ",\n\t\t "), times = ceiling(length(msgForm)/3))[1:length(msgForm)]
msgForm <- paste(msgForm, msgFormSep, sep = "", collapse = "")
msgForm <- gsub("\\,\\\n\\\t\\\t $|\\,$", "", msgForm)
message("\tFormula: sum(", msgForm, ")\n")
......@@ -70,10 +66,10 @@ ErrorCrit <- function(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbo
OutputsCrit <- listOutputsCrit
class(OutputsCrit) <- c("Multi", "ErrorCrit")
}
}
return(OutputsCrit)
}
R/ErrorCrit_KGE.R
View file @ 4bccfee2
ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "KGE", OutputsModel = OutputsModel, warnings = warnings)
CritValue <- NA
SubCritValues <- rep(NA, 3)
SubCritNames <- c("r", "alpha", "beta")
SubCritPrint <- rep(NA, 3)
if (EC$CritCompute) {
## Other variables preparation
meanVarObs <- mean(EC$VarObs[!EC$TS_ignore])
meanVarSim <- mean(EC$VarSim[!EC$TS_ignore])
## SubErrorCrit KGE rPearson
SubCritPrint[1L] <- paste0(EC$CritName, " cor(sim, obs, \"pearson\") =")
Numer <- sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) * (EC$VarSim[!EC$TS_ignore] - meanVarSim))
Deno1 <- sqrt(sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) ^ 2))
Deno2 <- sqrt(sum((EC$VarSim[!EC$TS_ignore] - meanVarSim) ^ 2))
if (Numer == 0) {
if (Deno1 == 0 & Deno2 == 0) {
Crit <- 1
......@@ -31,18 +31,18 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
Crit <- 0
}
} else {
Crit <- Numer / (Deno1 * Deno2)
Crit <- Numer / (Deno1 * Deno2)
}
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[1L] <- Crit
}
## SubErrorCrit KGE alpha
SubCritPrint[2L] <- paste0(EC$CritName, " sd(sim)/sd(obs) =")
Numer <- sd(EC$VarSim[!EC$TS_ignore])
Denom <- sd(EC$VarObs[!EC$TS_ignore])
if (Numer == 0 & Denom == 0) {
Crit <- 1
} else {
......@@ -51,10 +51,10 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[2L] <- Crit
}
## SubErrorCrit KGE beta
SubCritPrint[3L] <- paste0(EC$CritName, " mean(sim)/mean(obs) =")
if (meanVarSim == 0 & meanVarObs == 0) {
Crit <- 1
} else {
......@@ -63,20 +63,20 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[3L] <- Crit
}
## ErrorCrit
if (sum(is.na(SubCritValues)) == 0) {
CritValue <- (1 - sqrt((SubCritValues[1L] - 1)^2 + (SubCritValues[2L] - 1)^2 + (SubCritValues[3L] - 1)^2))
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " "))
}
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
......@@ -85,8 +85,11 @@ ErrorCrit_KGE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("KGE", "ErrorCrit")
return(OutputsCrit)
}
class(ErrorCrit_KGE) <- c("FUN_CRIT", class(ErrorCrit_KGE))
R/ErrorCrit_KGE2.R
View file @ 4bccfee2
ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "KGE2", OutputsModel = OutputsModel, warnings = warnings)
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "KGE2", OutputsModel = OutputsModel, warnings = warnings)
CritValue <- NA
SubCritValues <- rep(NA, 3)
SubCritNames <- c("r", "gamma", "beta")
SubCritPrint <- rep(NA, 3)
if (EC$CritCompute) {
## Other variables preparation
meanVarObs <- mean(EC$VarObs[!EC$TS_ignore])
meanVarSim <- mean(EC$VarSim[!EC$TS_ignore])
## SubErrorCrit KGE rPearson
SubCritPrint[1L] <- paste0(EC$CritName, " cor(sim, obs, \"pearson\") =")
Numer <- sum((EC$VarObs[!EC$TS_ignore] - meanVarObs) * (EC$VarSim[!EC$TS_ignore] - meanVarSim))
Deno1 <- sqrt(sum((EC$VarObs[!EC$TS_ignore] - meanVarObs)^2))
Deno2 <- sqrt(sum((EC$VarSim[!EC$TS_ignore] - meanVarSim)^2))
if (Numer == 0) {
if (Deno1 == 0 & Deno2 == 0) {
Crit <- 1
......@@ -31,15 +31,15 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
Crit <- 0
}
} else {
Crit <- Numer / (Deno1 * Deno2)
Crit <- Numer / (Deno1 * Deno2)
}
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[1L] <- Crit
}
## SubErrorCrit KGE gamma
SubCritPrint[2L] <- paste0(EC$CritName, " cv(sim)/cv(obs) =")
if (meanVarSim == 0) {
if (sd(EC$VarSim[!EC$TS_ignore]) == 0) {
CVsim <- 1
......@@ -48,7 +48,7 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
}
} else {
CVsim <- sd(EC$VarSim[!EC$TS_ignore]) / meanVarSim
}
if (meanVarObs == 0) {
if (sd(EC$VarObs[!EC$TS_ignore]) == 0) {
......@@ -68,10 +68,10 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[2L] <- Crit
}
## SubErrorCrit KGE beta
SubCritPrint[3L] <- paste0(EC$CritName, " mean(sim)/mean(obs) =")
if (meanVarSim == 0 & meanVarObs == 0) {
Crit <- 1
} else {
......@@ -80,20 +80,20 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
if (is.numeric(Crit) & is.finite(Crit)) {
SubCritValues[3L] <- Crit
}
## ErrorCrit
if (sum(is.na(SubCritValues)) == 0) {
CritValue <- (1 - sqrt((SubCritValues[1L] - 1)^2 + (SubCritValues[2L] - 1)^2 + (SubCritValues[3L] - 1)^2))
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " "))
}
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
......@@ -102,8 +102,10 @@ ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("KGE2", "ErrorCrit")
return(OutputsCrit)
}
class(ErrorCrit_KGE2) <- c("FUN_CRIT", class(ErrorCrit_KGE2))
R/ErrorCrit_NSE.R
View file @ 4bccfee2
ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "NSE", OutputsModel = OutputsModel, warnings = warnings)
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "NSE", OutputsModel = OutputsModel, warnings = warnings)
CritValue <- NA
if (EC$CritCompute) {
## Other variables preparation
meanVarObs <- mean(EC$VarObs[!EC$TS_ignore])
meanVarSim <- mean(EC$VarSim[!EC$TS_ignore])
## ErrorCrit
Emod <- sum((EC$VarSim[!EC$TS_ignore] - EC$VarObs[!EC$TS_ignore])^2)
Eref <- sum((EC$VarObs[!EC$TS_ignore] - mean(EC$VarObs[!EC$TS_ignore]))^2)
if (Emod == 0 & Eref == 0) {
Crit <- 0
} else {
......@@ -26,22 +22,24 @@ ErrorCrit_NSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = T
if (is.numeric(Crit) & is.finite(Crit)) {
CritValue <- Crit
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
}
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("NSE", "ErrorCrit")
return(OutputsCrit)
}
class(ErrorCrit_NSE) <- c("FUN_CRIT", class(ErrorCrit_NSE))
R/ErrorCrit_RMSE.R
View file @ 4bccfee2
ErrorCrit_RMSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) {
## Arguments check
if (!inherits(OutputsModel, "OutputsModel")) {
stop("'OutputsModel' must be of class 'OutputsModel'")
}
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "RMSE", OutputsModel = OutputsModel, warnings = warnings)
EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "RMSE", OutputsModel = OutputsModel, warnings = warnings)
CritValue <- NA
if (EC$CritCompute) {
## ErrorCrit
Numer <- sum((EC$VarSim - EC$VarObs)^2, na.rm = TRUE)
Denom <- sum(!is.na(EC$VarObs))
if (Numer == 0) {
Crit <- 0
} else {
......@@ -22,22 +22,24 @@ ErrorCrit_RMSE <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose =
if (is.numeric(Crit) & is.finite(Crit)) {
CritValue <- Crit
}
## Verbose
if (verbose) {
message(sprintf("Crit. %s = %.4f", EC$CritName, CritValue))
}
}
## Output
OutputsCrit <- list(CritValue = CritValue,
CritName = EC$CritName,
CritBestValue = EC$CritBestValue,
Multiplier = EC$Multiplier,
Ind_notcomputed = EC$Ind_TS_ignore)
class(OutputsCrit) <- c("RMSE", "ErrorCrit")
return(OutputsCrit)
}
class(ErrorCrit_RMSE) <- c("FUN_CRIT", class(ErrorCrit_RMSE))
R/Imax.R 0 → 100644
View file @ 4bccfee2
Imax <- function(InputsModel,
IndPeriod_Run,
TestedValues = seq(from = 0.1, to = 3, by = 0.1)) {
## ---------- check arguments
## InputsModel
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
if (!inherits(InputsModel, "hourly")) {
stop("'InputsModel' must be of class 'hourly'")
}
## IndPeriod_Run
if (!is.vector(IndPeriod_Run)) {
stop("'IndPeriod_Run' must be a vector of numeric values")
}
if (!inherits(IndPeriod_Run, "integer")) {
stop("'IndPeriod_Run' must be of type integer")
}
if (!identical(as.integer(IndPeriod_Run), IndPeriod_Run[1]:IndPeriod_Run[length(IndPeriod_Run)])) {
stop("'IndPeriod_Run' must be a continuous sequence of integers")
}
## TestedValues
if (!(is.numeric(TestedValues))) {
stop("'TestedValues' must be 'numeric'")
}
## ---------- hourly inputs aggregation
## aggregate data at the daily time step
daily_data <- SeriesAggreg(InputsModel[IndPeriod_Run], Format = "%Y%m%d")
## ---------- calculate interception
## calculate total interception of daily GR models on the period
cum_daily <- sum(pmin(daily_data$Precip, daily_data$PotEvap))
if (anyNA(cum_daily)) {
stop("'IndPeriod_Run' must be set to select 24 hours by day")
}
## calculate total interception of the GR5H interception store on the period
## and compute difference with daily values
differences <- array(NA, c(length(TestedValues)))
for (Imax in TestedValues) {
C0 <- 0
cum_hourly <- 0
for (i in IndPeriod_Run) {
Ec <- min(InputsModel$PotEvap[i], InputsModel$Precip[i] + C0)
Pth <- max(0, InputsModel$Precip[i] - (Imax-C0) - Ec)
C0 <- C0 + InputsModel$Precip[i] - Ec - Pth
cum_hourly <- cum_hourly + Ec
}
differences[which(Imax == TestedValues)] <- abs(cum_hourly - cum_daily)
}
## return the Imax value that minimises the difference
return(TestedValues[which.min(differences)])
}
R/PE_Oudin.R 0 → 100644
View file @ 4bccfee2
PE_Oudin <- function(JD, Temp,
Lat, LatUnit = c("rad", "deg"),
TimeStepIn = "daily", TimeStepOut = "daily",
RunFortran = FALSE) {
## ---------- check arguments
if (!(inherits(JD, "numeric") | inherits(JD, "integer"))) {
stop("'JD' must be of class 'numeric'")
}
if (!(inherits(Temp, "numeric") | inherits(Temp, "integer"))) {
stop("'Temp' must be of class 'numeric'")
}
if (length(JD) != length(Temp)) {
stop("'JD' and 'Temp' must have the same length")
}
if (!RunFortran & (!inherits(Lat, "numeric") | length(Lat) != 1)) {
stop("'Lat' must be a 'numeric' of length one")
}
if (RunFortran & (!inherits(Lat, "numeric") | (!length(Lat) %in% c(1, length(Temp))))) {
stop("'Lat' must be a 'numeric' of length one or of the same length as 'Temp'")
}
LatUnit <- match.arg(LatUnit, choices = c("rad", "deg"))
if (LatUnit[1L] == "rad" & (all(Lat >= pi/2) | all(Lat <= -pi/2))) {
stop("'Lat' must be comprised between -pi/2 and +pi/2 degrees")
}
if (LatUnit[1L] == "deg" & (all(Lat >= 90) | all(Lat <= -90))) {
stop("'Lat' must be comprised between -90 and +90 degrees")
}
if (LatUnit[1L] == "rad") {
FI <- Lat
}
if (LatUnit[1L] == "deg") {
FI <- Lat / (180 / pi)
}
if (any(JD < 0) | any(JD > 366)) {
stop("'JD' must only contain integers from 1 to 366")
}
TimeStep <- c("daily", "hourly")
TimeStepIn <- match.arg(TimeStepIn , choices = TimeStep)
TimeStepOut <- match.arg(TimeStepOut, choices = TimeStep)
rleJD <- rle(JD)
msgDaliy <- "each day should have only one identical value of Julian days. The time series is not sorted, or contains duplicate or missing dates"
msgHourly <- "each day must have 24 identical values of Julian days (one for each hour). The time series is not sorted, or contains duplicate or missing dates"
if (TimeStepIn == "daily" & any(rleJD$lengths != 1)) {
warning(msgDaliy)
}
## ---------- hourly inputs aggregation
if (TimeStepIn == "hourly") {
JD <- rleJD$values
idJD <- rep(seq_along(JD), each = rleJD$lengths[1L])
if (length(Temp) != length(idJD)) {
stop(msgHourly)
} else {
Temp <- as.vector(tapply(X = Temp, INDEX = idJD, FUN = mean))
}
}
if (TimeStepIn == "hourly" & any(rleJD$lengths != 24)) {
warning(msgHourly)
}
## ---------- Oudin's formula
if (RunFortran) {
LInputs = as.integer(length(Temp))
if (length(FI) == 1) {
FI <- rep(FI, LInputs)
}
RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR",
##inputs
LInputs = LInputs,
InputsLAT = as.double(FI),
InputsTT = as.double(Temp),
InputsJJ = as.double(JD),
##outputs
PE_Oudin_D = rep(as.double(-99e9), LInputs)
)
PE_Oudin_D = RESULTS$PE_Oudin_D
} else {
PE_Oudin_D <- rep(NA, length(Temp))
COSFI <- cos(FI)
for (k in seq_along(Temp)) {
TETA <- 0.4093 * sin(JD[k] / 58.1 - 1.405)
COSTETA <- cos(TETA)
COSGZ <- max(0.001, cos(FI - TETA))
GZ <- acos(COSGZ)
COSOM <- 1 - COSGZ / COSFI / COSTETA
if (COSOM < -1) {
COSOM <- -1
}
if (COSOM > 1) {
COSOM <- 1
}
COSOM2 <- COSOM * COSOM
if (COSOM2 >= 1) {
SINOM <- 0
} else {
SINOM <- sqrt(1 - COSOM2)
}
OM <- acos(COSOM)
COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1)
if (COSPZ < 0.001) {
COSPZ <- 0.001
}
ETA <- 1 + cos(JD[k] / 58.1) / 30
GE <- 446 * OM * COSPZ * ETA
if (is.na(Temp[k])) {
PE_Oudin_D[k] <- NA
} else {
if (Temp[k] >= -5.0) {
PE_Oudin_D[k] <- GE * (Temp[k] + 5) / 100 / 28.5
} else {
PE_Oudin_D[k] <- 0
}
}
}
}
## ---------- disaggregate PE from daily to hourly
if (TimeStepOut == "hourly") {
parab_D2H <- c(0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
0.035, 0.062, 0.079, 0.097, 0.110, 0.117,
0.117, 0.110, 0.097, 0.079, 0.062, 0.035,
0.000, 0.000, 0.000, 0.000, 0.000, 0.000)
PE_Oudin_H <- rep(PE_Oudin_D, each = 24) * rep(parab_D2H, times = length(PE_Oudin_D))
}
## ---------- outputs warnings
if (any(is.na(Temp))) {
if (any(is.na(PE_Oudin_D))) {
warning("'Temp' time series, and therefore the returned 'PE' time series, contain missing value(s)")
} else {
warning("'Temp' time series contains missing value(s)")
}
}
if (!any(is.na(Temp)) & any(is.na(PE_Oudin_D))) {
warning("returned 'PE' time series contains missing value(s)")
}
if (TimeStepOut == "daily") {
PE_Oudin <- PE_Oudin_D
} else {
PE_Oudin <- PE_Oudin_H
}
return(PE_Oudin)
}
R/PEdaily_Oudin.R deleted 100644 → 0
View file @ 0967defb
PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg")) {
## ---------- check arguments
if (!missing(LatRad)) {
warning("Deprecated \"LatRad\" argument. Please, use \"Lat\" instead.")
if (missing(Lat)) {
Lat <- LatRad
}
}
if (!(inherits(JD, "numeric") | inherits(JD, "integer"))) {
stop("'JD' must be of class 'numeric'")
}
if (!(inherits(Temp, "numeric") | inherits(Temp, "integer"))) {
stop("'Temp' must be of class 'numeric'")
}
if (length(JD) != length(Temp)) {
stop("'Temp' and 'LatUnit' must have the same length")
}
if (!any(LatUnit %in% c("rad", "deg"))) {
stop("'LatUnit' must be one of \"rad\" or \"deg\"")
}
if (!inherits(Lat, "numeric") | length(Lat) != 1) {
stop("'Lat' must be a 'numeric' of length one")
}
if (LatUnit[1L] == "rad" & ((Lat >= pi/2) | (Lat <= -pi/2))) {
stop("'Lat' must be comprised between -pi/2 and +pi/2 degrees")
}
if (LatUnit[1L] == "deg" & ((Lat >= 90) | (Lat <= -90))) {
stop("'Lat' must be comprised between -90 and +90 degrees")
}
if (LatUnit[1L] == "rad") {
FI <- Lat
}
if (LatUnit[1L] == "deg") {
FI <- Lat / (180 / pi)
}
if (any(JD < 0) | any(JD > 366)) {
stop("'JD' must only contain integers from 1 to 366")
}
## ---------- Oudin's formula
PE_Oudin_D <- rep(NA, length(Temp))
COSFI <- cos(FI)
AFI <- abs(FI / 42)
for (k in seq_along(Temp)) {
TETA <- 0.4093 * sin(JD[k] / 58.1 - 1.405)
COSTETA <- cos(TETA)
COSGZ <- max(0.001, cos(FI - TETA))
GZ <- acos(COSGZ)
COSOM <- 1 - COSGZ / COSFI / COSTETA
if (COSOM < -1) {
COSOM <- -1
}
if (COSOM > 1) {
COSOM <- 1
}
COSOM2 <- COSOM * COSOM
if (COSOM2 >= 1) {
SINOM <- 0
} else {
SINOM <- sqrt(1 - COSOM2)
}
OM <- acos(COSOM)
COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1)
if (COSPZ < 0.001) {
COSPZ <- 0.001
}
ETA <- 1 + cos(JD[k] / 58.1) / 30
GE <- 446 * OM * COSPZ * ETA
if (is.na(Temp[k])) {
PE_Oudin_D[k] <- NA
} else {
if (Temp[k] >= -5.0) {
PE_Oudin_D[k] <- GE * (Temp[k] + 5) / 100 / 28.5
} else {
PE_Oudin_D[k] <- 0
}
}
}
if (any(is.na(Temp))) {
if (any(is.na(PE_Oudin_D))) {
warning("'Temp' time series, and therefore the returned 'PE' time series, contain missing value(s)")
} else {
warning("'Temp' time series contains missing value(s)")
}
}
if (!any(is.na(Temp)) & any(is.na(PE_Oudin_D))) {
warning("returned 'PE' time series contains missing value(s)")
}
return(PE_Oudin_D)
}
R/RunModel.R
View file @ 4bccfee2
RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD) {
RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, ...) {
FUN_MOD <- match.fun(FUN_MOD)
return(FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param))
if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
# Lag model take one parameter at the beginning of the vector
iFirstParamRunOffModel <- 2
RunOptions$FeatFUN_MOD$NbParam <- RunOptions$FeatFUN_MOD$NbParam - 1
} else {
# All parameters
iFirstParamRunOffModel <- 1
}
OutputsModel <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param[iFirstParamRunOffModel:length(Param)], ...)
if (inherits(InputsModel, "SD") && !identical(FUN_MOD, RunModel_Lag)) {
OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel)
}
return(OutputsModel)
}
R/RunModel_CemaNeige.R
View file @ 4bccfee2
......@@ -3,17 +3,17 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
NParam <- ifelse(test = IsHyst, yes = 4L, no = 2L)
NParamCN <- ifelse(test = IsHyst, yes = 4L, no = 2L)
NStates <- 4L
FortranOutputsCemaNeige <- .FortranOutputs(GR = NULL, isCN = TRUE)$CN
## Arguments_check
## Arguments check
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
if (!inherits(InputsModel, "daily")) {
stop("'InputsModel' must be of class 'daily'")
if (!inherits(InputsModel, "daily") & !inherits(InputsModel, "hourly")) {
stop("'InputsModel' must be of class 'daily' or 'hourly'")
}
if (!inherits(InputsModel, "CemaNeige")) {
stop("'InputsModel' must be of class 'CemaNeige'")
......@@ -27,24 +27,30 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
if (!is.vector(Param) | !is.numeric(Param)) {
stop("'Param' must be a numeric vector")
}
if (sum(!is.na(Param)) != NParam) {
stop(sprintf("'Param' must be a vector of length %i and contain no NA", NParam))
if (sum(!is.na(Param)) != NParamCN) {
stop(sprintf("'Param' must be a vector of length %i and contain no NA", NParamCN))
}
## Input_data_preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0)) {
if (inherits(InputsModel, "daily")) {
time_step <- "daily"
time_mult <- 1
}
if (inherits(InputsModel, "hourly")) {
time_step <- "hourly"
time_mult <- 24
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp) + 1):length(IndPeriod1)
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):length(IndPeriod1)
## Output data preparation
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## SNOW_MODULE________________________________________________________________________________
## CemaNeige________________________________________________________________________________
ParamCemaNeige <- Param
NLayers <- length(InputsModel$LayerPrecip)
......@@ -57,55 +63,48 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
} else {
IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call_DLL_CemaNeige_________________________
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*time_mult + 40*time_mult) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*time_mult + 40*time_mult) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_CemaNeige", PACKAGE = "airGR",
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = as.integer(length(IndPeriod1)), ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/time step]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParam), ### number of model parameter
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(-999.999, ### output series [mm]
nrow = length(IndPeriod1),
ncol = length(IndOutputsCemaNeige)),
StateEnd = rep(-999.999, NStates) ### state variables at the end of the model run (reservoir levels [mm] and HU)
Outputs = matrix(as.double(-99e9), nrow = length(IndPeriod1), ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/time step or degC]
StateEnd = rep(as.double(-99e9), NStates) ### state variables at the end of the model run
)
RESULTS$Outputs[round(RESULTS$Outputs , 3) == -999.999] <- NA
RESULTS$StateEnd[round(RESULTS$StateEnd, 3) == -999.999] <- NA
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data_storage
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige]
if (ExportStateEnd) {
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
rm(RESULTS)
} ### ENDFOR_iLayer
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
......@@ -114,28 +113,32 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
CemaNeigeStateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeige, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = NULL, RoutStore = NULL, ExpStore = NULL,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
## Output_data_preparation
## Output data preparation
## OutputsModel only
if (!ExportDatesR & !ExportStateEnd) {
OutputsModel <- list(CemaNeigeLayers)
names(OutputsModel) <- NameCemaNeigeLayers
}
## DatesR and OutputsModel only
if (ExportDatesR & !ExportStateEnd) {
OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
list(CemaNeigeLayers))
names(OutputsModel) <- c("DatesR", NameCemaNeigeLayers)
}
## OutputsModel and StateEnd only
if (!ExportDatesR & ExportStateEnd) {
OutputsModel <- c(list(CemaNeigeLayers),
list(CemaNeigeStateEnd))
names(OutputsModel) <- c(NameCemaNeigeLayers, "StateEnd")
}
## DatesR and OutputsModel and StateEnd
if (ExportDatesR & ExportStateEnd) {
OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
list(CemaNeigeLayers),
......@@ -143,11 +146,11 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
names(OutputsModel) <- c("DatesR", NameCemaNeigeLayers, "StateEnd")
}
## End
class(OutputsModel) <- c("OutputsModel", "daily", "CemaNeige")
if(IsHyst) {
class(OutputsModel) <- c("OutputsModel", time_step, "CemaNeige")
if (IsHyst) {
class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
}
return(OutputsModel)
}
R/RunModel_CemaNeigeGR4H.R 0 → 100644
View file @ 4bccfee2
RunModel_CemaNeigeGR4H <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 4L
NStates <- 4L
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
## Output data preparation
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*24 + 40*24) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*24 + 40*24) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/h]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/h or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (iLayer > 1) {
CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (ExportStateEnd) {
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
rm(RESULTS)
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ### ENDIF RunSnowModule
if (!inherits(RunOptions, "CemaNeige")) {
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- NULL
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
}
## GR model
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr4h", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/h]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam = NParamMod, ### number of model parameter
Param = ParamMod, ### parameter set
NStates = NStatesMod, ### number of state variables used for model initialising
StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs = IndOutputsMod, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/h]
StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4H, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:(20*24)) + 7],
UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
}
R/RunModel_CemaNeigeGR4J.R
View file @ 4bccfee2
RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
RunModel_CemaNeigeGR4J <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
NParam <- ifelse(test = IsHyst, yes = 8L, no = 6L)
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 4L
NStates <- 4L
FortranOutputs <- .FortranOutputs(GR = "GR4J", isCN = TRUE)
##Arguments_check
if(!inherits(InputsModel,"InputsModel")){ stop("'InputsModel' must be of class 'InputsModel'") }
if(!inherits(InputsModel,"daily" )){ stop("'InputsModel' must be of class 'daily' ") }
if(!inherits(InputsModel,"GR" )){ stop("'InputsModel' must be of class 'GR' ") }
if(!inherits(InputsModel,"CemaNeige" )){ stop("'InputsModel' must be of class 'CemaNeige' ") }
if(!inherits(RunOptions,"RunOptions" )){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(!inherits(RunOptions,"GR" )){ stop("'RunOptions' must be of class 'GR' ") }
if(!inherits(RunOptions,"CemaNeige" )){ stop("'RunOptions' must be of class 'CemaNeige' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
## Output data preparation
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/d or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ParamCemaNeige <- Param[(length(Param)-1-2*as.integer(IsHyst)):length(Param)];
NParamMod <- as.integer(length(Param)-(2+2*as.integer(IsHyst)));
ParamMod <- Param[1:NParamMod];
NLayers <- length(InputsModel$LayerPrecip);
NStatesMod <- as.integer(length(RunOptions$IniStates)-NStates*NLayers);
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##SNOW_MODULE________________________________________________________________________________##
if(inherits(RunOptions,"CemaNeige")){
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN));
} else { IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam=as.integer(NParam), ### number of model parameter = 2
Param=as.double(ParamCemaNeige), ### parameter set
NStates=as.integer(NStates), ### number of state variables used for model initialising = 2
StateStart=StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs=as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs=IndOutputsCemaNeige, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
StateEnd=rep(as.double(-999.999),as.integer(NStates)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Data_storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige];
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ###ENDIF_RunSnowModule
if(!inherits(RunOptions,"CemaNeige")){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]; }
##MODEL______________________________________________________________________________________##
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputs$GR));
} else { IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*ParamMod[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*ParamMod[3]; ### routing store level (mm)
if (iLayer > 1) {
CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
##Call_fortan
RESULTS <- .Fortran("frun_GR4J",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=CatchMeltAndPliq, ### input series of total precipitation [mm/d]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=NParamMod, ### number of model parameter
Param=ParamMod, ### parameter set
NStates=NStatesMod, ### number of state variables used for model initialising
StateStart=RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs=IndOutputsMod, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsMod)), ### output series [mm]
StateEnd=rep(as.double(-999.999),NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
if(inherits(RunOptions,"CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputs$GR[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
##Output_data_preparation
##OutputsModel_only
if(!ExportDatesR & !ExportStateEnd){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers) );
names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
##DatesR_and_OutputsModel_only
if( ExportDatesR & !ExportStateEnd){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers) );
names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
##OutputsModel_and_SateEnd_only
if(!ExportDatesR & ExportStateEnd){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if( ExportDatesR & ExportStateEnd){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); }
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige");
if(IsHyst) {
class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
}
return(OutputsModel);
rm(RESULTS)
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ### ENDIF RunSnowModule
if (!inherits(RunOptions, "CemaNeige")) {
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- NULL
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
}
## GR model______________________________________________________________________________________
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr4j", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/d]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam = NParamMod, ### number of model parameter
Param = ParamMod, ### parameter set
NStates = NStatesMod, ### number of state variables used for model initialising
StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs = IndOutputsMod, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/d]
StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:20) + 7],
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
}
R/RunModel_CemaNeigeGR5H.R 0 → 100644
View file @ 4bccfee2
RunModel_CemaNeigeGR5H <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 5L
NStates <- 4L
IsIntStore <- inherits(RunOptions, "interception")
if (IsIntStore) {
Imax <- RunOptions$Imax
} else {
Imax <- -99
}
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
## Output data preparation
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*24 + 40*24) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20*24 + 40*24) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/h]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/h or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (iLayer > 1) {
CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if (ExportStateEnd) {
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
rm(RESULTS)
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ### ENDIF RunSnowModule
if (!inherits(RunOptions, "CemaNeige")) {
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- NULL
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
}
## GR model
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
if (IsIntStore) {
RunOptions$IniStates[4] <- RunOptions$IniResLevels[4] * Imax ### interception store level (mm)
}
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr5h", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/h]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam = NParamMod, ### number of model parameter
Param = ParamMod, ### parameter set
NStates = NStatesMod, ### number of state variables used for model initialising
StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
Imax = Imax, ### maximal capacity of interception store
NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs = IndOutputsMod, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/h]
StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR5H, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
IntStore = RESULTS$StateEnd[4L],
UH1 = NULL,
UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
}
R/RunModel_CemaNeigeGR5J.R
View file @ 4bccfee2
RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
RunModel_CemaNeigeGR5J <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
NParam <- ifelse(test = IsHyst, yes = 9L, no = 7L)
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 5L
NStates <- 4L
FortranOutputs <- .FortranOutputs(GR = "GR5J", isCN = TRUE)
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }
if(inherits(InputsModel,"daily" )==FALSE){ stop("'InputsModel' must be of class 'daily' ") }
if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") }
if(inherits(InputsModel,"CemaNeige" )==FALSE){ stop("'InputsModel' must be of class 'CemaNeige' ") }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") }
if(inherits(RunOptions,"CemaNeige" )==FALSE){ stop("'RunOptions' must be of class 'CemaNeige' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ParamCemaNeige <- Param[(length(Param)-1-2*as.integer(IsHyst)):length(Param)];
NParamMod <- as.integer(length(Param)-(2+2*as.integer(IsHyst)));
ParamMod <- Param[1:NParamMod];
NLayers <- length(InputsModel$LayerPrecip);
NStatesMod <- as.integer(length(RunOptions$IniStates)-NStates*NLayers);
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##SNOW_MODULE________________________________________________________________________________##
if(inherits(RunOptions,"CemaNeige")==TRUE){
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN));
} else { IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam=as.integer(NParam), ### number of model parameter = 2
Param=as.double(ParamCemaNeige), ### parameter set
NStates=as.integer(NStates), ### number of state variables used for model initialising = 2
StateStart=StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs=as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs=IndOutputsCemaNeige, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
StateEnd=rep(as.double(-999.999),as.integer(NStates)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Data_storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige];
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ###ENDIF_RunSnowModule
if(inherits(RunOptions,"CemaNeige")==FALSE){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]; }
##MODEL______________________________________________________________________________________##
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputs$GR));
} else { IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*ParamMod[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*ParamMod[3]; ### routing store level (mm)
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/d or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
##Call_fortan
RESULTS <- .Fortran("frun_GR5J",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=CatchMeltAndPliq, ### input series of total precipitation [mm/d]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=NParamMod, ### number of model parameter
Param=ParamMod, ### parameter set
NStates=NStatesMod, ### number of state variables used for model initialising
StateStart=RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs=IndOutputsMod, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsMod)), ### output series [mm]
StateEnd=rep(as.double(-999.999),NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR5J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL, UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputs$GR[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
##Output_data_preparation
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers) );
names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers) );
names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if(ExportDatesR==TRUE & ExportStateEnd==TRUE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); }
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige");
if(IsHyst) {
class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
if (iLayer > 1) {
CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
return(OutputsModel);
if (ExportStateEnd) {
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
rm(RESULTS)
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ### ENDIF RunSnowModule
if (!inherits(RunOptions, "CemaNeige")) {
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- NULL
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
}
## GR model______________________________________________________________________________________
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr5j", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/d]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam = NParamMod, ### number of model parameter
Param = ParamMod, ### parameter set
NStates = NStatesMod, ### number of state variables used for model initialising
StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs = IndOutputsMod, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/d]
StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR5J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL,
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
}
R/RunModel_CemaNeigeGR6J.R
View file @ 4bccfee2
RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
RunModel_CemaNeigeGR6J <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsHyst <- inherits(RunOptions, "hysteresis")
NParam <- ifelse(test = IsHyst, yes = 10L, no = 8L)
NParamCN <- RunOptions$FeatFUN_MOD$NbParam - 6L
NStates <- 4L
FortranOutputs <- .FortranOutputs(GR = "GR6J", isCN = TRUE)
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }
if(inherits(InputsModel,"daily" )==FALSE){ stop("'InputsModel' must be of class 'daily' ") }
if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") }
if(inherits(InputsModel,"CemaNeige" )==FALSE){ stop("'InputsModel' must be of class 'CemaNeige' ") }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") }
if(inherits(RunOptions,"CemaNeige" )==FALSE){ stop("'RunOptions' must be of class 'CemaNeige' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
Param_X1X3X6_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[1L] <- Param_X1X3X6_threshold
}
if (Param[3L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[3L] <- Param_X1X3X6_threshold
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3X6_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[1L] <- Param_X1X3X6_threshold
}
if (Param[3L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[3L] <- Param_X1X3X6_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
if (Param[6L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[6L] <- Param_X1X3X6_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ParamCemaNeige <- Param[(length(Param) - 1 - 2 * as.integer(IsHyst)):length(Param)]
NParamMod <- as.integer(length(Param) - (2 + 2 * as.integer(IsHyst)))
ParamMod <- Param[1:NParamMod]
NLayers <- length(InputsModel$LayerPrecip)
NStatesMod <- as.integer(length(RunOptions$IniStates) - NStates * NLayers)
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## CemaNeige________________________________________________________________________________
if (inherits(RunOptions, "CemaNeige")) {
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsCemaNeige <- as.integer(1:length(RunOptions$FortranOutputs$CN))
} else {
IndOutputsCemaNeige <- which(RunOptions$FortranOutputs$CN %in% RunOptions$Outputs_Sim)
}
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- "CemaNeigeLayers"
## Call CemaNeige Fortran_________________________
for (iLayer in 1:NLayers) {
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
if (Param[6L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[6L] <- Param_X1X3X6_threshold
}
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ParamCemaNeige <- Param[(length(Param)-1-2*as.integer(IsHyst)):length(Param)];
NParamMod <- as.integer(length(Param)-(2+2*as.integer(IsHyst)));
ParamMod <- Param[1:NParamMod];
NLayers <- length(InputsModel$LayerPrecip);
NStatesMod <- as.integer(length(RunOptions$IniStates)-NStates*NLayers);
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##SNOW_MODULE________________________________________________________________________________##
if(inherits(RunOptions,"CemaNeige")==TRUE){
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN));
} else { IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
if (!IsHyst) {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers)]
} else {
StateStartCemaNeige <- RunOptions$IniStates[(7 + 20 + 40) + c(iLayer, iLayer+NLayers, iLayer+2*NLayers, iLayer+3*NLayers)]
}
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam=as.integer(NParam), ### number of model parameter = 2
Param=as.double(ParamCemaNeige), ### parameter set
NStates=as.integer(NStates), ### number of state variables used for model initialising = 2
StateStart=StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs=as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs=IndOutputsCemaNeige, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
StateEnd=rep(as.double(-999.999),as.integer(NStates)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
##Data_storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige];
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ###ENDIF_RunSnowModule
if(inherits(RunOptions,"CemaNeige")==FALSE){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]; }
##MODEL______________________________________________________________________________________##
if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputs$GR));
} else { IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
RunOptions$IniStates[3] <- RunOptions$IniResLevels[3] ### exponential store level (mm)
RESULTS <- .Fortran("frun_cemaneige", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$LayerPrecip[[iLayer]][IndPeriod1], ### input series of total precipitation [mm/d]
InputsFracSolidPrecip = InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
InputsTemp = InputsModel$LayerTemp[[iLayer]][IndPeriod1], ### input series of air mean temperature [degC]
MeanAnSolidPrecip = RunOptions$MeanAnSolidPrecip[iLayer], ### value of annual mean solid precip [mm/y]
NParam = as.integer(NParamCN), ### number of model parameters = 2 or 4
Param = as.double(ParamCemaNeige), ### parameter set
NStates = as.integer(NStates), ### number of state variables used for model initialising = 4
StateStart = StateStartCemaNeige, ### state variables used when the model run starts
IsHyst = as.integer(IsHyst), ### use of hysteresis
NOutputs = as.integer(length(IndOutputsCemaNeige)), ### number of output series
IndOutputs = IndOutputsCemaNeige, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsCemaNeige)), ### output series [mm, mm/d or degC]
StateEnd = rep(as.double(-99e9), as.integer(NStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## Data storage
CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
names(CemaNeigeLayers[[iLayer]]) <- RunOptions$FortranOutputs$CN[IndOutputsCemaNeige]
IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt")
if (iLayer == 1) {
CatchMeltAndPliq <- RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
##Call_fortan
RESULTS <- .Fortran("frun_GR6J",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=CatchMeltAndPliq, ### input series of total precipitation [mm/d]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=NParamMod, ### number of model parameter
Param=ParamMod, ### parameter set
NStates=NStatesMod, ### number of state variables used for model initialising
StateStart=RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs=IndOutputsMod, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsMod)), ### output series [mm]
StateEnd=rep(as.double(-999.999),NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
if (iLayer > 1) {
CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[, IndPliqAndMelt] / NLayers
}
if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputs$GR[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
##Output_data_preparation
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers) );
names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers) );
names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if(ExportDatesR==TRUE & ExportStateEnd==TRUE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); }
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige");
if(IsHyst) {
class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
if (ExportStateEnd) {
CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd)
}
return(OutputsModel);
rm(RESULTS)
} ### ENDFOR iLayer
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ### ENDIF RunSnowModule
if (!inherits(RunOptions, "CemaNeige")) {
CemaNeigeLayers <- list()
CemaNeigeStateEnd <- NULL
NameCemaNeigeLayers <- NULL
CatchMeltAndPliq <- InputsModel$Precip[IndPeriod1]
}
## GR model______________________________________________________________________________________
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputsMod <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputsMod <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
RunOptions$IniStates[3] <- RunOptions$IniResLevels[3] ### exponential store level (mm)
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr6j", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = CatchMeltAndPliq, ### input series of total precipitation [mm/d]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam = NParamMod, ### number of model parameter
Param = ParamMod, ### parameter set
NStates = NStatesMod, ### number of state variables used for model initialising
StateStart = RunOptions$IniStates[1:NStatesMod], ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputsMod)), ### number of output series
IndOutputs = IndOutputsMod, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputsMod)), ### output series [mm or mm/d]
StateEnd = rep(as.double(-99e9), NStatesMod) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel, IsHyst = IsHyst,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
UH1 = RESULTS$StateEnd[(1:20) + 7],
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 1]],
eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 2]],
GthrCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 3]],
GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
verbose = FALSE)
}
if (inherits(RunOptions, "CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim) {
RESULTS$Outputs[, which(RunOptions$FortranOutputs$GR[IndOutputsMod] == "Precip")] <-
InputsModel$Precip[IndPeriod1]
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers)
}
R/RunModel_GR1A.R
View file @ 4bccfee2
RunModel_GR1A <- function(InputsModel,RunOptions,Param){
NParam <- 1;
FortranOutputs <- .FortranOutputs(GR = "GR1A")$GR
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }
if(inherits(InputsModel,"yearly" )==FALSE){ stop("'InputsModel' must be of class 'yearly' ") }
if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs));
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); }
##Output_data_preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
BOOL_Fortran <- FALSE; if(BOOL_Fortran){
##Call_fortan
RESULTS <- .Fortran("frun_gr1a",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/y]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/y]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
} else {
##R_version
L <- length(IndPeriod1)
P0 <- InputsModel$Precip[ IndPeriod1][1:(L-1)]
P1 <- InputsModel$Precip[ IndPeriod1][2: L ]
E1 <- InputsModel$PotEvap[IndPeriod1][2: L ]
Q1 <- P1*(1.-1./(1.+((0.7*P1+0.3*P0)/Param[1]/E1)^2.0)^0.5)
PEQ <- rbind(c(NA,NA,NA),cbind(P1,E1,Q1))
Outputs <- PEQ[,IndOutputs]
if(is.vector(Outputs)){ Outputs <- cbind(Outputs); }
RESULTS <- list(NOutputs=length(IndOutputs),IndOutputs=IndOutputs,Outputs=Outputs,StatesEnd=NA)
}
##Output_data_preparation
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","yearly","GR");
return(OutputsModel);
RunModel_GR1A <- function(InputsModel, RunOptions, Param) {
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr1a", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/y]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/y]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm/y]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}
R/RunModel_GR2M.R
View file @ 4bccfee2
RunModel_GR2M <- function(InputsModel,RunOptions,Param){
RunModel_GR2M <- function(InputsModel, RunOptions, Param) {
NParam <- 2;
FortranOutputs <- .FortranOutputs(GR = "GR2M")$GR
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }
if(inherits(InputsModel,"monthly" )==FALSE){ stop("'InputsModel' must be of class 'monthly' ") }
if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
Param_X1X2_threshold <- 1e-2
if (Param[1L] < Param_X1X2_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X2_threshold, Param_X1X2_threshold))
Param[1L] <- Param_X1X2_threshold
}
if (Param[2L] < Param_X1X2_threshold) {
warning(sprintf("Param[2] (X2: routing store capacity [mm]) < %.2f\n X2 set to %.2f", Param_X1X2_threshold, Param_X1X2_threshold))
Param[2L] <- Param_X1X2_threshold
}
Param <- as.double(Param)
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs));
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); }
Param_X1X2_threshold <- 1e-2
if (Param[1L] < Param_X1X2_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X2_threshold, Param_X1X2_threshold))
Param[1L] <- Param_X1X2_threshold
}
if (Param[2L] < Param_X1X2_threshold) {
warning(sprintf("Param[2] (X2: routing store capacity [mm]) < %.2f\n X2 set to %.2f", Param_X1X2_threshold, Param_X1X2_threshold))
Param[2L] <- Param_X1X2_threshold
}
##Output_data_preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##Use_of_IniResLevels
if(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[2]; ### routing store level (mm)
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Call_fortan
RESULTS <- .Fortran("frun_GR2M",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/month]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/month]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs [round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
##Output_data_preparation
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[2] ### routing store level (mm)
}
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","monthly","GR");
return(OutputsModel);
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr2m", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/month]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/month]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/month]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL,
UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## Output data preparation
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}
R/RunModel_GR4H.R
View file @ 4bccfee2
RunModel_GR4H <- function(InputsModel,RunOptions,Param){
RunModel_GR4H <- function(InputsModel, RunOptions, Param) {
NParam <- 4;
FortranOutputs <- .FortranOutputs(GR = "GR4H")$GR
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }
if(inherits(InputsModel,"hourly" )==FALSE){ stop("'InputsModel' must be of class 'hourly' ") }
if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
Param <- as.double(Param)
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs));
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); }
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Output_data_preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##Use_of_IniResLevels
if(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm)
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Call_fortan
RESULTS <- .Fortran("frun_GR4H",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4H, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:(20*24))+7], UH2 = RESULTS$StateEnd[(1:(40*24))+(7+20*24)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
##Output_data_preparation
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
}
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","hourly","GR");
return(OutputsModel);
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr4h", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/h]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4H, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:(20*24)) + 7],
UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}
R/RunModel_GR4J.R
View file @ 4bccfee2
RunModel_GR4J <- function(InputsModel,RunOptions,Param){
RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
NParam <- 4;
FortranOutputs <- .FortranOutputs(GR = "GR4J")$GR
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }
if(inherits(InputsModel,"daily" )==FALSE){ stop("'InputsModel' must be of class 'daily' ") }
if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
Param <- as.double(Param)
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs));
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); }
##Input_data_preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Use_of_IniResLevels
if(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm)
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Call_fortan
RESULTS <- .Fortran("frun_GR4J",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
}
##Output_data_preparation
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_StateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_StateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR");
return(OutputsModel);
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr4j", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/d]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = RESULTS$StateEnd[(1:20) + 7],
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}
R/RunModel_GR5H.R 0 → 100644
View file @ 4bccfee2
RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
## Initialization of variables
IsIntStore <- inherits(RunOptions, "interception")
if (IsIntStore) {
Imax <- RunOptions$Imax
} else {
Imax <- -99
}
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
Param <- as.double(Param)
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
if (IsIntStore) {
RunOptions$IniStates[4] <- RunOptions$IniResLevels[4] * Imax ### interception store level (mm)
}
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr5h", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
Imax = Imax, ### maximal capacity of interception store
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/h]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
IntStore = RESULTS$StateEnd[4L],
UH1 = NULL,
UH2 = RESULTS$StateEnd[(1:(40*24)) + (7+20*24)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}