diff --git a/DESCRIPTION b/DESCRIPTION
index f857e81d07104e2536e2349f8e6f2f7a155d9b94..ea9b15ef120b2bd281d9fac9058535ce4b4dcaaa 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.1.3.11
-Date: 2019-02-21
+Version: 1.2.6.0
+Date: 2019-02-25
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
diff --git a/NAMESPACE b/NAMESPACE
index d5fab141d7ddcbae2a3827463337b265b578c7e7..031f782463a365c03c27be6f04c2bdd97d228b39 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -43,6 +43,7 @@ export(RunModel_GR6J)
export(SeriesAggreg)
export(TransfoParam)
export(TransfoParam_CemaNeige)
+export(TransfoParam_CemaNeigeHyst)
export(TransfoParam_GR1A)
export(TransfoParam_GR2M)
export(TransfoParam_GR4H)
diff --git a/NEWS.rmd b/NEWS.rmd
index 3ae207dc73d6fafe38f3f50ebb6340cfed763449..e8c94884a0e04f5d8c6cbff166c7fd556c9dd749 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -13,7 +13,7 @@ output:
-### 1.1.3.11 Release Notes (2019-02-21)
+### 1.2.6.0 Release Notes (2019-02-25)
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index 4b58ed75c0f32f05250e7e2be52816d891be006b..b21367c0538ba8a3c6835c32ab43ecc3d96e0483 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -54,19 +54,25 @@ Calibration_Michel <- function(InputsModel, RunOptions, InputsCrit, CalibOptions
FUN_TRANSFO <- TransfoParam_GR1A
}
if (identical(FUN_MOD, RunModel_CemaNeige )) {
- FUN_TRANSFO <- TransfoParam_CemaNeige
+ if (inherits(FUN_MOD, "hysteresis")) {
+ FUN_TRANSFO <- TransfoParam_CemaNeigeHyst
+ } else {
+ FUN_TRANSFO <- TransfoParam_CemaNeige
+ }
}
if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
if (identical(FUN_MOD, RunModel_CemaNeigeGR4J)) {
FUN1 <- TransfoParam_GR4J
- FUN2 <- TransfoParam_CemaNeige
}
if (identical(FUN_MOD, RunModel_CemaNeigeGR5J)) {
FUN1 <- TransfoParam_GR5J
- FUN2 <- TransfoParam_CemaNeige
}
if (identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
FUN1 <- TransfoParam_GR6J
+ }
+ if (inherits(FUN_MOD, "hysteresis")) {
+ FUN2 <- TransfoParam_CemaNeigeHyst
+ } else {
FUN2 <- TransfoParam_CemaNeige
}
FUN_TRANSFO <- function(ParamIn, Direction) {
diff --git a/R/CreateIniStates.R b/R/CreateIniStates.R
index 384f4ad8f90b37eb118b79f875531a879dbc0882..a5bc0d095ccf847f09fd3c7a3a51ab891bd60176 100644
--- a/R/CreateIniStates.R
+++ b/R/CreateIniStates.R
@@ -2,6 +2,7 @@ CreateIniStates <- function(FUN_MOD, InputsModel,
ProdStore = 350, RoutStore = 90, ExpStore = NULL,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+ GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
verbose = TRUE) {
@@ -181,7 +182,13 @@ CreateIniStates <- function(FUN_MOD, InputsModel,
if (is.null(eTGCemaNeigeLayers)) {
eTGCemaNeigeLayers <- rep(Inf, NLayers)
}
-
+ if (is.null(GthrCemaNeigeLayers)) {
+ GthrCemaNeigeLayers <- rep(Inf, NLayers)
+ }
+ if (is.null(GlocmaxCemaNeigeLayers)) {
+ GlocmaxCemaNeigeLayers <- rep(Inf, NLayers)
+ }
+ cat(GCemaNeigeLayers)
# check negative values
if (any(ProdStore < 0) | any(RoutStore < 0) |
@@ -224,7 +231,8 @@ CreateIniStates <- function(FUN_MOD, InputsModel,
## format output
IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore),
UH = list(UH1 = UH1, UH2 = UH2),
- CemaNeigeLayers = list(G = GCemaNeigeLayers, eTG = eTGCemaNeigeLayers))
+ CemaNeigeLayers = list(G = GCemaNeigeLayers, eTG = eTGCemaNeigeLayers,
+ Gthr = GthrCemaNeigeLayers, Glocmax = GlocmaxCemaNeigeLayers))
IniStatesNA <- unlist(IniStates)
IniStatesNA[is.infinite(IniStatesNA)] <- NA
IniStatesNA <- relist(IniStatesNA, skeleton = IniStates)
diff --git a/R/CreateInputsCrit.R b/R/CreateInputsCrit.R
index 269eb92b49718af6ab45a0d0cb2950e181a6b8e2..15ef5c538ac29a6928c14de3a305c3d5c9df79eb 100644
--- a/R/CreateInputsCrit.R
+++ b/R/CreateInputsCrit.R
@@ -148,12 +148,6 @@ CreateInputsCrit <- function(FUN_CRIT,
}
## check 'obs'
- # lapply(iListArgs2$obs, function(iObs) {
- # if (!is.vector(iObs) | length(iObs) != LLL | !is.numeric(iObs)) {
- # stop(sprintf("'obs' must be a (list of) vector(s) of numeric values of length %i \n", LLL), call. = FALSE)
- # return(NULL)
- # }
- # })
if (!is.vector(iListArgs2$obs) | length(iListArgs2$obs) != LLL | !is.numeric(iListArgs2$obs)) {
stop(sprintf("'obs' must be a (list of) vector(s) of numeric values of length %i \n", LLL), call. = FALSE)
return(NULL)
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index bf3d5ebbb44668423356a39043be45bc92bb7d63..b01041181219206e56feccffb86dd79d8e7a026f 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -208,7 +208,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
NState <- 7 + 3 * 24 * 20
}
if ("daily" %in% ObjectClass) {
- NState <- 7 + 3 * 20 + 2 * NLayers
+ NState <- 7 + 3 * 20 + 4 * NLayers
}
if ("monthly" %in% ObjectClass) {
NState <- 2
@@ -290,7 +290,7 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Qsim")
}
if ("CemaNeige" %in% ObjectClass) {
- Outputs_all <- c(Outputs_all,"Pliq", "Psol", "SnowPack", "ThermalState", "Gratio", "PotMelt", "Melt", "PliqAndMelt", "Temp")
+ Outputs_all <- c(Outputs_all,"Pliq", "Psol", "SnowPack", "ThermalState", "Gratio", "PotMelt", "Melt", "PliqAndMelt", "Temp", "Gthreshold", "Glocalmax")
}
##check_Outputs_Sim
diff --git a/R/RunModel_CemaNeige.R b/R/RunModel_CemaNeige.R
index 4eec9007db5c75e7af63b04c4eebf5e5b7dcba67..13e6321b4c79b3e9dd143ae518a63a5a6e94c5b2 100644
--- a/R/RunModel_CemaNeige.R
+++ b/R/RunModel_CemaNeige.R
@@ -1,10 +1,20 @@
-RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
+RunModel_CemaNeige <- function(InputsModel, RunOptions, Param, IsHyst = FALSE) {
+
+ ## Arguments_check
+ if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
+ stop("'IsHyst' must be a 'logical' of length 1")
+ return(NULL)
+ }
- NParam <- 2L
+
+ ## Initialization of variables
+ NParam <- ifelse(IsHyst, 4L, 2L)
+ NStates <- 4L
FortranOutputsCemaNeige <- c("Pliq", "Psol",
"SnowPack", "ThermalState", "Gratio",
- "PotMelt", "Melt", "PliqAndMelt", "Temp")
+ "PotMelt", "Melt", "PliqAndMelt", "Temp",
+ "Gthreshold", "Glocalmax")
@@ -47,7 +57,7 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-
+
@@ -81,17 +91,18 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
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(2), ### number of model parameter = 2
+ NParam = as.integer(NParam), ### number of model parameter
Param = as.double(ParamCemaNeige), ### parameter set
- NStates = as.integer(2), ### number of state variables used for model initialising = 2
+ NStates = as.integer(NStates), ### number of state variables used for model initialising
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, 2) ### state variables at the end of the model run (reservoir levels [mm] and HU)
+ StateEnd = rep(-999.999, 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
@@ -112,11 +123,14 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
if (ExportStateEnd) {
+ idNStates <- seq_len(NStates*NLayers) %% NStates
CemaNeigeStateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeige, InputsModel = InputsModel,
ProdStore = NULL, RoutStore = NULL, ExpStore = NULL,
UH1 = NULL, UH2 = NULL,
- GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 1]],
- eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 0]],
+ 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)
}
@@ -145,6 +159,9 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) {
## End
class(OutputsModel) <- c("OutputsModel", "daily", "CemaNeige")
+ if(IsHyst) {
+ class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
+ }
return(OutputsModel)
}
diff --git a/R/RunModel_CemaNeigeGR4J.R b/R/RunModel_CemaNeigeGR4J.R
index 732b52bb913bf890f9c9bb769b035adbdebc9f77..c426f90575c74440b9283d76f9a7cce58cf35e90 100644
--- a/R/RunModel_CemaNeigeGR4J.R
+++ b/R/RunModel_CemaNeigeGR4J.R
@@ -1,7 +1,15 @@
-RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
+RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE){
- NParam <- 6;
- FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp");
+
+ ## Arguments_check
+ if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
+ stop("'IsHyst' must be a 'logical' of length 1")
+ return(NULL)
+ }
+
+ NParam <- ifelse(IsHyst, 8L, 6L)
+ NStates <- 4L
+ FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp", "Gthreshold", "Glocalmax");
FortranOutputsMod <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
"AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
@@ -37,11 +45,11 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
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):length(Param)];
- NParamMod <- as.integer(length(Param)-2);
+ 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)-2*NLayers);
+ NStatesMod <- as.integer(length(RunOptions$IniStates)-NStates*NLayers);
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
@@ -52,6 +60,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
} else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
+
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
@@ -62,20 +71,20 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
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(2), ### number of model parameter = 2
- Param=ParamCemaNeige, ### parameter set
- NStates=as.integer(2), ### number of state variables used for model initialising = 2
+ 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(2)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
+ 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]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
@@ -85,7 +94,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
- names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep="");
+ names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ###ENDIF_RunSnowModule
if(inherits(RunOptions,"CemaNeige")==FALSE){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
@@ -121,12 +130,15 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
- if (ExportStateEnd) {
+ if (ExportStateEnd) {
+ idNStates <- seq_len(NStates*NLayers) %% NStates
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4J, 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 = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 1]],
- eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 0]],
+ 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)
}
@@ -161,6 +173,9 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige");
+ if(IsHyst) {
+ class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
+ }
return(OutputsModel);
}
diff --git a/R/RunModel_CemaNeigeGR5J.R b/R/RunModel_CemaNeigeGR5J.R
index 5457034a5448479d8497f05a45d59bce1be39e2d..fc58b9c9a32668f33d92ebd0350f2647a5b830a2 100644
--- a/R/RunModel_CemaNeigeGR5J.R
+++ b/R/RunModel_CemaNeigeGR5J.R
@@ -1,9 +1,17 @@
-RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
-
- NParam <- 7;
- FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp");
+RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE){
+
+
+ ## Arguments_check
+ if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
+ stop("'IsHyst' must be a 'logical' of length 1")
+ return(NULL)
+ }
+
+ NParam <- ifelse(IsHyst, 9L, 7L)
+ NStates <- 4L
+ FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp", "Gthreshold", "Glocalmax");
FortranOutputsMod <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
- "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
+ "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
@@ -37,11 +45,11 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
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):length(Param)];
- NParamMod <- as.integer(length(Param)-2);
+ 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)-2*NLayers);
+ NStatesMod <- as.integer(length(RunOptions$IniStates)-NStates*NLayers);
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
@@ -52,6 +60,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
} else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
+
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
@@ -62,15 +71,16 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
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(2), ### number of model parameter = 2
- Param=ParamCemaNeige, ### parameter set
- NStates=as.integer(2), ### number of state variables used for model initialising = 2
+ 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(2)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
+ 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;
@@ -84,7 +94,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
- names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep="");
+ names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ###ENDIF_RunSnowModule
if(inherits(RunOptions,"CemaNeige")==FALSE){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
@@ -121,11 +131,14 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
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,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL, UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
- GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 1]],
- eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 0]],
+ 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)
}
@@ -160,6 +173,9 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige");
+ if(IsHyst) {
+ class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
+ }
return(OutputsModel);
}
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index 0bbef76f9cbb9e1f9d9fea39e874ac6b329d43ff..862dad0eccb3a5e016a7714ffe3ba00470b3e6fc 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -1,7 +1,15 @@
-RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
-
- NParam <- 8;
- FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp");
+RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE){
+
+
+ ## Arguments_check
+ if (!is.logical(IsHyst) | length(IsHyst) != 1L) {
+ stop("'IsHyst' must be a 'logical' of length 1")
+ return(NULL)
+ }
+
+ NParam <- ifelse(IsHyst, 10L, 8L)
+ NStates <- 4L
+ FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp", "Gthreshold", "Glocalmax");
FortranOutputsMod <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1",
"Rout", "Exch", "AExch1", "AExch2", "AExch", "QR", "QRExp", "Exp", "QD", "Qsim");
@@ -41,22 +49,22 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
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):length(Param)];
- NParamMod <- as.integer(length(Param)-2);
+ 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)-2*NLayers);
+ 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(FortranOutputsCemaNeige));
} else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim); }
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
+
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
@@ -67,15 +75,16 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
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(2), ### number of model parameter = 2
- Param=ParamCemaNeige, ### parameter set
- NStates=as.integer(2), ### number of state variables used for model initialising = 2
+ 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(2)) ### state variables at the end of the model run (reservoir levels [mm] and HU)
+ 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;
@@ -89,7 +98,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); }
rm(RESULTS);
} ###ENDFOR_iLayer
- names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep="");
+ names(CemaNeigeLayers) <- sprintf("Layer%02i", seq_len(NLayers))
} ###ENDIF_RunSnowModule
if(inherits(RunOptions,"CemaNeige")==FALSE){
CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- NULL;
@@ -127,11 +136,14 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
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,
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(2*NLayers)[seq_len(2*NLayers) %% 2 == 1]],
- eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %% 2 == 0]],
+ 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)
}
@@ -166,6 +178,9 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR","CemaNeige");
+ if(IsHyst) {
+ class(OutputsModel) <- c(class(OutputsModel), "hysteresis")
+ }
return(OutputsModel);
}
diff --git a/R/TransfoParam_CemaNeigeHyst.R b/R/TransfoParam_CemaNeigeHyst.R
new file mode 100644
index 0000000000000000000000000000000000000000..09a889342026cd58205ac59a58375ce98e2366db
--- /dev/null
+++ b/R/TransfoParam_CemaNeigeHyst.R
@@ -0,0 +1,46 @@
+TransfoParam_CemaNeigeHyst <- function(ParamIn, Direction) {
+
+
+ NParam <- 4L
+
+
+ ## check_arguments
+ Bool <- is.vector(ParamIn)
+ if (Bool) {
+ ParamIn <- matrix(ParamIn, nrow = 1)
+ warning("'ParamIn' automatically convert to 'matrix'")
+ }
+ if (!inherits(ParamIn, "matrix")) {
+ stop("'ParamIn' must be of class 'matrix'")
+ }
+ if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
+ stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
+ }
+ if (ncol(ParamIn) != NParam) {
+ stop(sprintf( "the CemaNeige module with hysteresis requires %i parameters", NParam))
+ }
+
+
+ if (Direction == "TR") {
+ ParamOut <- ParamIn
+ ParamOut[, 1] <- (ParamIn[, 1] + 9.99) / 19.98 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
+ ParamOut[, 2] <- exp(ParamIn[, 2]) / 200 ### CemaNeige X2 (degree-day melt coefficient)
+ ParamOut[, 3] <- (ParamIn[, 3] * 5) + 50 ### Hyst Gaccum
+ ParamOut[, 4] <- (ParamIn[, 4] / 19.98) + 0.5 ### Hyst CV
+ }
+ if (Direction == "RT") {
+ ParamOut <- ParamIn
+ ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
+ ParamOut[, 2] <- log(ParamIn[, 2] * 200) ### CemaNeige X2 (degree-day melt coefficient)
+ ParamOut[, 3] <- (ParamIn[, 3] - 50) / 5 ### Hyst Gaccum
+ ParamOut[, 4] <- (ParamIn[, 4] - 0.5) * 19.98 ### Hyst CV
+ }
+
+ if (Bool) {
+ ParamOut <- as.vector(ParamOut)
+ }
+
+ return(ParamOut)
+
+}
+
diff --git a/man/CreateIniStates.Rd b/man/CreateIniStates.Rd
index aed2025f4623a2d460cd5a6a91c8d258d707a72c..685efb964677d71b3b64cdf2ede5dc85c80648f8 100644
--- a/man/CreateIniStates.Rd
+++ b/man/CreateIniStates.Rd
@@ -13,6 +13,7 @@ CreateIniStates(FUN_MOD, InputsModel,
ProdStore = 350, RoutStore = 90, ExpStore = NULL,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+ GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
verbose = TRUE)
}
@@ -36,6 +37,10 @@ CreateIniStates(FUN_MOD, InputsModel,
\item{eTGCemaNeigeLayers}{(optional) [numeric] snow pack thermal state [°C], possibly used to create the CemaNeige model initial state}
+\item{GthrCemaNeigeLayers}{}
+
+\item{GlocmaxCemaNeigeLayers}{}
+
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
}
diff --git a/man/RunModel_CemaNeige.Rd b/man/RunModel_CemaNeige.Rd
index dee53c8aa08c798951626754461db681df9f55e8..0da708fe53c5e74285610e77dda1464e5bb5d9b6 100644
--- a/man/RunModel_CemaNeige.Rd
+++ b/man/RunModel_CemaNeige.Rd
@@ -9,7 +9,7 @@
\usage{
-RunModel_CemaNeige(InputsModel, RunOptions, Param)
+RunModel_CemaNeige(InputsModel, RunOptions, Param, IsHyst = FALSE)
}
@@ -23,6 +23,8 @@ RunModel_CemaNeige(InputsModel, RunOptions, Param)
CemaNeige X1 \tab weighting coefficient for snow pack thermal state [-] \cr
CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d] \cr
}}
+
+\item{IsHyst}{}
}
diff --git a/man/RunModel_CemaNeigeGR4J.Rd b/man/RunModel_CemaNeigeGR4J.Rd
index eb646e31613134dd3204d9136b07d96f0e80337b..938767709b98c062f175865ca1ccb036de5dec8a 100644
--- a/man/RunModel_CemaNeigeGR4J.Rd
+++ b/man/RunModel_CemaNeigeGR4J.Rd
@@ -9,7 +9,7 @@
\usage{
-RunModel_CemaNeigeGR4J(InputsModel, RunOptions, Param)
+RunModel_CemaNeigeGR4J(InputsModel, RunOptions, Param, IsHyst = FALSE)
}
@@ -27,6 +27,8 @@ GR4J X4 \tab unit hydrograph time constant [d]
CemaNeige X1 \tab weighting coefficient for snow pack thermal state [-] \cr
CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d] \cr
}}
+
+\item{IsHyst}{}
}
diff --git a/man/RunModel_CemaNeigeGR5J.Rd b/man/RunModel_CemaNeigeGR5J.Rd
index c5d21ce347a1c23b42b22d7021e7f75f3b6b6a17..c6e0ae18e3f82ee103472946db845a79bf81dfe0 100644
--- a/man/RunModel_CemaNeigeGR5J.Rd
+++ b/man/RunModel_CemaNeigeGR5J.Rd
@@ -9,7 +9,7 @@
\usage{
-RunModel_CemaNeigeGR5J(InputsModel, RunOptions, Param)
+RunModel_CemaNeigeGR5J(InputsModel, RunOptions, Param, IsHyst = FALSE)
}
@@ -28,6 +28,7 @@ GR5J X5 \tab intercatchment exchange threshold [-] \cr
CemaNeige X1 \tab weighting coefficient for snow pack thermal state [-] \cr
CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d] \cr
}}
+\item{IsHyst}{}
}
diff --git a/man/RunModel_CemaNeigeGR6J.Rd b/man/RunModel_CemaNeigeGR6J.Rd
index 73df2b7c9362c6d4cdee9cffc03ea33e40f6a517..38fc40136a5372a9158e9888e0448e576f52d18b 100644
--- a/man/RunModel_CemaNeigeGR6J.Rd
+++ b/man/RunModel_CemaNeigeGR6J.Rd
@@ -9,7 +9,7 @@
\usage{
-RunModel_CemaNeigeGR6J(InputsModel, RunOptions, Param)
+RunModel_CemaNeigeGR6J(InputsModel, RunOptions, Param, IsHyst = FALSE)
}
@@ -29,6 +29,7 @@ GR6J X6 \tab coefficient for emptying exponential store [mm] \cr
CemaNeige X1 \tab weighting coefficient for snow pack thermal state [-] \cr
CemaNeige X2 \tab degree-day melt coefficient [mm/°C/d] \cr
}}
+\item{IsHyst}{}
}
diff --git a/src/airGR.c b/src/airGR.c
index 9c149457a5d36fe9bea49df47c352c49f82c59d9..4ccbc8789ca287a401b04a745ea070ea577e698d 100644
--- a/src/airGR.c
+++ b/src/airGR.c
@@ -7,7 +7,7 @@
*/
/* .Fortran calls */
-extern void F77_NAME(frun_cemaneige)(int *, double *, double *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
+extern void F77_NAME(frun_cemaneige)(int *, double *, double *, double *, double *, int *, double *, int *, double *, int *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr1a)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr2m)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr4h)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
@@ -16,7 +16,7 @@ extern void F77_NAME(frun_gr5j)(int *, double *, double *, int *, double *, int
extern void F77_NAME(frun_gr6j)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
static const R_FortranMethodDef FortranEntries[] = {
- {"frun_cemaneige", (DL_FUNC) &F77_NAME(frun_cemaneige), 13},
+ {"frun_cemaneige", (DL_FUNC) &F77_NAME(frun_cemaneige), 14},
{"frun_gr1a", (DL_FUNC) &F77_NAME(frun_gr1a), 11},
{"frun_gr2m", (DL_FUNC) &F77_NAME(frun_gr2m), 11},
{"frun_gr4h", (DL_FUNC) &F77_NAME(frun_gr4h), 11},
diff --git a/src/frun_CEMANEIGE.f b/src/frun_CEMANEIGE.f
index 3e66f254f3ba47c76e7ec27b469a0d78065cd39f..13dd51578177b31802da64b7e6e3e51557a63ba0 100644
--- a/src/frun_CEMANEIGE.f
+++ b/src/frun_CEMANEIGE.f
@@ -1,5 +1,7 @@
+
+
SUBROUTINE frun_CEMANEIGE(
!inputs
& LInputs , ! [integer] length of input and output series
@@ -11,6 +13,7 @@
& Param , ! [double] parameter set
& NStates , ! [integer] number of state variables used for model initialising = 2
& StateStart , ! [double] state variables used when the model run starts
+ & IsHyst , ! [logical] whether we should use the linear hysteresis or not
& NOutputs , ! [integer] number of output series
& IndOutputs , ! [integer] indices of output series
!outputs
@@ -25,39 +28,53 @@
!### input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, intent(in) :: MeanAnSolidPrecip
- doubleprecision, dimension(LInputs) :: InputsPrecip
- doubleprecision, dimension(LInputs) :: InputsFracSolidPrecip
- doubleprecision, dimension(LInputs) :: InputsTemp
- doubleprecision, dimension(NParam) :: Param
- doubleprecision, dimension(NStates) :: StateStart
- doubleprecision, dimension(NStates) :: StateEnd
- integer, dimension(NOutputs) :: IndOutputs
- doubleprecision, dimension(LInputs,NOutputs) :: Outputs
+ doubleprecision, intent(in), dimension(LInputs) :: InputsPrecip
+ doubleprecision, intent(in), dimension(LInputs) ::
+ & InputsFracSolidPrecip
+ doubleprecision, intent(in), dimension(LInputs) :: InputsTemp
+ doubleprecision, intent(in), dimension(NParam) :: Param
+ doubleprecision, intent(in), dimension(NStates) :: StateStart
+ logical, intent(in) :: IsHyst ! TRUE if linear hysteresis is used, FALSE otherwise
+ doubleprecision, intent(out), dimension(NStates) :: StateEnd
+ integer, intent(in), dimension(NOutputs) :: IndOutputs
+ doubleprecision, intent(out), dimension(LInputs,NOutputs) ::
+ & Outputs
!parameters, internal states and variables
doubleprecision CTG,Kf
- doubleprecision G,eTG,PliqAndMelt
- doubleprecision Tmelt,Gthreshold,MinSpeed
- doubleprecision Pliq,Psol,Gratio,PotMelt,Melt
+ doubleprecision G,eTG,PliqAndMelt,dG,prct
+ doubleprecision Tmelt,Gthreshold,MinSpeed,Gacc,Glocalmax
+ doubleprecision Pliq,Psol,Gratio,PotMelt,Melt,Ginit
integer I,K
!--------------------------------------------------------------
!Initialisations
!--------------------------------------------------------------
- !initilisation des constantes
+ !initialisation of constants
Tmelt=0.
- Gthreshold=0.9*MeanAnSolidPrecip
MinSpeed=0.1
- !initilisation of model states using StateStart
+ !initialisation of model states using StateStart
G=StateStart(1)
eTG=StateStart(2)
+ Gthreshold=StateStart(3)
+ Glocalmax=StateStart(4)
+ Gratio=0.
PliqAndMelt=0.
!setting parameter values
CTG=Param(1)
Kf=Param(2)
+ IF (IsHyst) THEN
+ Gacc=Param(3)
+ prct=Param(4)
+
+ Gthreshold=prct*MeanAnSolidPrecip
+ Glocalmax=Gthreshold
+ ELSE
+ Gthreshold=0.9*MeanAnSolidPrecip
+ ENDIF
!initialisation of model outputs
c StateEnd = -999.999 !initialisation made in R
@@ -71,14 +88,16 @@ c Outputs = -999.999 !initialisation made in R
DO k=1,LInputs
!SolidPrecip and LiquidPrecip
- Pliq=(1-InputsFracSolidPrecip(k))*InputsPrecip(k)
+ Pliq=(1.-InputsFracSolidPrecip(k))*InputsPrecip(k)
Psol=InputsFracSolidPrecip(k)*InputsPrecip(k)
!Snow pack volume before melt
+ IF (IsHyst) Ginit=G
G=G+Psol
+
!Snow pack thermal state before melt
- eTG=CTG*eTG + (1-CTG)*InputsTemp(k)
+ eTG=CTG*eTG + (1.-CTG)*InputsTemp(k)
IF(eTG.GT.0.) eTG=0.
!Potential melt
@@ -89,46 +108,69 @@ c Outputs = -999.999 !initialisation made in R
PotMelt=0.
ENDIF
- !Gratio
- IF(G.LT.Gthreshold) THEN
- Gratio=G/Gthreshold
+
+ IF (IsHyst) THEN
+ IF (Potmelt.GT.0.) THEN
+ IF (G.LT.Glocalmax.AND.Gratio.EQ.1.) Glocalmax=G !Update in case of potential melt and G lower than Gseuil
+ Gratio=MIN(G/Glocalmax,1.)
+ ENDIF
ELSE
- Gratio=1.
+ IF(G.LT.Gthreshold) THEN
+ Gratio=G/Gthreshold
+ ELSE
+ Gratio=1.
+ ENDIF
ENDIF
!Actual melt
- Melt=((1-MinSpeed)*Gratio+MinSpeed)*PotMelt
+ Melt=((1.-MinSpeed)*Gratio+MinSpeed)*PotMelt
!Update of snow pack volume
G=G-Melt
-
- !Update of Gratio
- IF(G.LT.Gthreshold) THEN
- Gratio=G/Gthreshold
+ IF (IsHyst) THEN
+ dG=G-Ginit !Melt in case of negative dG, accumulation otherwise
+
+
+ IF (dG.GT.0.) THEN
+ Gratio = MIN(Gratio+(Psol-Melt)/Gacc,1.) !Psol - Melt = dG
+ IF (Gratio.EQ.1.) Glocalmax = Gthreshold
+ ENDIF
+ IF (dG.LT.0.) THEN
+ Gratio=MIN(G/Glocalmax,1.)
+ ENDIF
ELSE
- Gratio=1.
+ IF(G.LT.Gthreshold) THEN
+ Gratio=G/Gthreshold
+ ELSE
+ Gratio=1.
+ ENDIF
ENDIF
+
!Water volume to pass to the hydrological model
PliqAndMelt=Pliq+Melt
!Storage of outputs
DO I=1,NOutputs
- IF(IndOutputs(I).EQ.1) Outputs(k,I)=Pliq ! Pliq ! observed liquid precipitation [mm/day]
- IF(IndOutputs(I).EQ.2) Outputs(k,I)=Psol ! Psol ! observed solid precipitation [mm/day]
- IF(IndOutputs(I).EQ.3) Outputs(k,I)=G ! SnowPack ! snow pack [mm]
- IF(IndOutputs(I).EQ.4) Outputs(k,I)=eTG ! ThermalState ! thermal state [°C]
- IF(IndOutputs(I).EQ.5) Outputs(k,I)=Gratio ! Gratio ! Gratio [-]
- IF(IndOutputs(I).EQ.6) Outputs(k,I)=PotMelt ! PotMelt ! potential snow melt [mm/day]
- IF(IndOutputs(I).EQ.7) Outputs(k,I)=Melt ! Melt ! melt [mm/day]
- IF(IndOutputs(I).EQ.8) Outputs(k,I)=PliqAndMelt ! PliqAndMelt ! liquid precipitation + melt [mm/day]
- IF(IndOutputs(I).EQ.9) Outputs(k,I)=InputsTemp(k) ! Temp ! air temperature [°C]
+ IF(IndOutputs(I).EQ.1) Outputs(k,I)=Pliq ! Pliq ! observed liquid precipitation [mm/day]
+ IF(IndOutputs(I).EQ.2) Outputs(k,I)=Psol ! Psol ! observed solid precipitation [mm/day]
+ IF(IndOutputs(I).EQ.3) Outputs(k,I)=G ! SnowPack ! snow pack [mm]
+ IF(IndOutputs(I).EQ.4) Outputs(k,I)=eTG ! ThermalState ! thermal state [°C]
+ IF(IndOutputs(I).EQ.5) Outputs(k,I)=Gratio ! Gratio ! Gratio [-]
+ IF(IndOutputs(I).EQ.6) Outputs(k,I)=PotMelt ! PotMelt ! potential snow melt [mm/day]
+ IF(IndOutputs(I).EQ.7) Outputs(k,I)=Melt ! Melt ! melt [mm/day]
+ IF(IndOutputs(I).EQ.8) Outputs(k,I)=PliqAndMelt ! PliqAndMelt ! liquid precipitation + melt [mm/day]
+ IF(IndOutputs(I).EQ.9) Outputs(k,I)=InputsTemp(k) ! Temp ! air temperature [°C]
+ IF(IndOutputs(I).EQ.10) Outputs(k,I)=Gthreshold ! Gthreshold ! melt threshold [mm]
+ IF(IndOutputs(I).EQ.11) Outputs(k,I)=Glocalmax ! Glocalmax ! local melt threshold for hysteresis [mm]
ENDDO
ENDDO
StateEnd(1)=G
StateEnd(2)=eTG
+ StateEnd(3)=Gthreshold
+ StateEnd(4)=Glocalmax
RETURN