Commit baca0e67 authored by unknown's avatar unknown
Browse files

v1.0.8.2 new function CreateIniStates to help user to format IniStates for CreateRunOptions #4422

parent 8acfff47
^.*\.Rproj$
^\.Rproj\.user$
^\.Rprofile$
^packrat/
......@@ -2,3 +2,4 @@
.Rhistory
.RData
airGR.Rproj
packrat/lib*/
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.0.8.1
Date: 2017-06-14
Version: 1.0.8.2
Date: 2017-06-21
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl")),
person("Charles", "Perrin", role = c("aut", "ths")),
......
......@@ -18,6 +18,7 @@ S3method("plot", "OutputsModel")
export(Calibration)
export(Calibration_Michel)
export(CreateCalibOptions)
export(CreateIniStates)
export(CreateInputsCrit)
export(CreateInputsModel)
export(CreateRunOptions)
......
CreateIniStates <- function(FUN_MOD, InputsModel,
ProdStore = 0.3, RoutStore = 0.5, ExpStore = NULL,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = TRUE) {
ObjectClass <- NULL
UH1n <- 20L
UH2n <- UH1n * 2L
## check FUN_MOD
BOOL <- FALSE
if (identical(FUN_MOD, RunModel_GR4H)) {
ObjectClass <- c(ObjectClass, "GR", "hourly")
BOOL <- TRUE
}
if (identical(FUN_MOD, RunModel_GR4J) |
identical(FUN_MOD, RunModel_GR5J) |
identical(FUN_MOD, RunModel_GR6J)) {
ObjectClass <- c(ObjectClass, "GR", "daily")
BOOL <- TRUE
}
if (identical(FUN_MOD, RunModel_GR2M)) {
ObjectClass <- c(ObjectClass, "GR", "monthly")
BOOL <- TRUE
}
if (identical(FUN_MOD, RunModel_GR1A)) {
stop("'RunModel_GR1A' does not require 'IniStates' object")
}
if (identical(FUN_MOD, RunModel_CemaNeige)) {
ObjectClass <- c(ObjectClass, "CemaNeige", "daily")
BOOL <- TRUE
}
if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "daily")
BOOL <- TRUE
}
if (!BOOL) {
stop("Incorrect 'FUN_MOD' for use in 'CreateIniStates'")
return(NULL)
}
## check InputsModel
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
return(NULL)
}
if ("GR" %in% ObjectClass & !inherits(InputsModel, "GR")) {
stop("'InputsModel' must be of class 'GR'")
return(NULL)
}
if ("CemaNeige" %in% ObjectClass &
!inherits(InputsModel, "CemaNeige")) {
stop("'RunModel_CemaNeigeGR*' must be of class 'CemaNeige'")
return(NULL)
}
## check states
if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
if (is.null(ExpStore)) {
stop("'RunModel_*GR6J' need an 'ExpStore' value")
return(NULL)
}
} else if (!is.null(ExpStore) & verbose) {
warning("This 'FUN_MOD' does not require 'ExpStore'. Value set to NA")
ExpStore <- Inf
}
if (identical(FUN_MOD, RunModel_GR2M) & verbose) {
if (!is.null(UH1)) {
warning("This 'FUN_MOD' does not require 'UH1'. Values set to NA")
UH1 <- rep(Inf, UH1n)
}
if (!is.null(UH2)) {
warning("This 'FUN_MOD' does not require 'UH2'. Values set to NA")
UH2 <- rep(Inf, UH2n)
}
}
if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !is.null(UH1) & verbose) {
warning("This 'FUN_MOD' does not require 'UH1'. Values set to NA")
UH1 <- rep(Inf, UH1n)
}
if("CemaNeige" %in% ObjectClass &
(is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers))) {
stop("'RunModel_CemaNeigeGR*' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'")
return(NULL)
}
if(!"CemaNeige" %in% ObjectClass &
(!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers)) & verbose) {
warning("This 'FUN_MOD' does not require 'GCemaNeigeLayers' and 'GCemaNeigeLayers'. Values set to NA")
GCemaNeigeLayers <- Inf
eTGCemaNeigeLayers <- Inf
}
## set states
if("CemaNeige" %in% ObjectClass) {
NLayers <- length(InputsModel$LayerPrecip)
} else {
NLayers <- 1
}
## manage NULL values
if (is.null(ExpStore)) {
ExpStore <- Inf
}
if (is.null(UH1)) {
if ("hourly" %in% ObjectClass) {
k <- 24
} else {
k <- 1
}
UH1 <- rep(Inf, 20 * k)
}
if (is.null(UH2)) {
if ("hourly" %in% ObjectClass) {
k <- 24
} else {
k <- 1
}
UH2 <- rep(Inf, 40 * k)
}
if (is.null(GCemaNeigeLayers)) {
GCemaNeigeLayers <- rep(Inf, NLayers)
}
if (is.null(eTGCemaNeigeLayers)) {
eTGCemaNeigeLayers <- rep(Inf, NLayers)
}
## check length
if (!is.numeric(ProdStore) || length(ProdStore) != 1L) {
stop("'ProdStore' must be numeric of length one")
}
if (!is.numeric(RoutStore) || length(RoutStore) != 1L) {
stop("'RoutStore' must be numeric of length one")
}
if (!is.numeric(ExpStore) || length(ExpStore) != 1L) {
stop("'ExpStore' must be numeric of length one")
}
if ( "hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != 480L)) {
stop(sprintf("'UH1' must be numeric of length 480 (%i * 24)", UH1n))
}
if (!"hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != 20L)) {
stop(sprintf("'UH1' must be numeric of length %i", UH1n))
}
if ( "hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != 960L)) {
stop(sprintf("'UH2' must be numeric of length 960 (%i * 24)", UH2n))
}
if (!"hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != 40L)) {
stop(sprintf("'UH2' must be numeric of length %i (2 * %i)", UH2n, UH1n))
}
if (!is.numeric(GCemaNeigeLayers) || length(GCemaNeigeLayers) != NLayers) {
stop(sprintf("'GCemaNeigeLayers' must be numeric of length %i", NLayers))
}
if (!is.numeric(eTGCemaNeigeLayers) || length(eTGCemaNeigeLayers) != NLayers) {
stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
}
# if ("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass) {
# if ("hourly" %in% ObjectClass) {
# NState <- 3 * 24 * 20 + 7
# }
# if ("daily" %in% ObjectClass) {
# if (identical(FUN_MOD, RunModel_GR5J)) {
# NState <-
# 2 * 20 + 2 * NLayers + 7
# } else {
# NState <- 3 * 20 + 2 * NLayers + 7
# }
# }
# if ("monthly" %in% ObjectClass) {
# NState <- 2
# }
# if ("yearly" %in% ObjectClass) {
# NState <- 1
# }
# }
# if (!is.null(IniStates)) {
# if (!is.vector(IniStates) | !is.numeric(IniStates)) {
# stop("IniStates must be a vector of numeric values")
# return(NULL)
# }
# if (length(IniStates) != NState) {
# stop(paste0(
# "The length of IniStates must be ",
# NState,
# " for the chosen FUN_MOD"
# ))
# return(NULL)
# }
# } else {
# IniStates <- as.double(rep(0.0, NState))
# IniStates[1:3] <- NA
# }
# if ("yearly" %in% ObjectClass) {
# IniStates <- c(ProdStore)
# }
# else if ("monthly" %in% ObjectClass) {
# IniStates <- c(ProdStore, RoutStore)
# }
IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore),
UH = list(UH1 = UH1, UH2 = UH2),
CemaNeigeLayers = list(G = GCemaNeigeLayers, eTG = eTGCemaNeigeLayers))
IniStatesNA <- unlist(IniStates)
IniStatesNA[is.infinite(IniStatesNA)] <- NA
IniStatesNA <- relist(IniStatesNA, skeleton = IniStates)
class(IniStatesNA) <- c("IniStates", ObjectClass)
return(IniStatesNA)
}
CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod_Run,IniStates=NULL,IniResLevels=NULL,
Outputs_Cal=NULL,Outputs_Sim="all",RunSnowModule=TRUE,MeanAnSolidPrecip=NULL, verbose = TRUE){
CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, IniStates = NULL, IniResLevels = NULL,
Outputs_Cal = NULL, Outputs_Sim = "all", RunSnowModule = TRUE, MeanAnSolidPrecip = NULL, verbose = TRUE) {
ObjectClass <- NULL;
ObjectClass <- NULL
##check_FUN_MOD
BOOL <- FALSE;
......@@ -93,37 +93,158 @@ CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod
}
if(!is.null(WTxt) & verbose){ warning(WTxt); }
##check_IniStates_and_IniResLevels
if(is.null(IniStates) & is.null(IniResLevels) & verbose){
warning("\t Model states initialisation not defined -> default configuration used \n"); }
## check IniResLevels
if ("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)) {
if (!is.null(IniResLevels)) {
if (!is.vector(IniResLevels) | !is.numeric(IniResLevels) | any(is.na(IniResLevels))) {
stop("IniResLevels must be a vector of numeric values \n")
return(NULL)
}
if ((identical(FUN_MOD, RunModel_GR4H) |
identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
identical(FUN_MOD, RunModel_GR2M)) &
length(IniResLevels) != 2) {
stop("The length of IniStates must be 2 for the chosen FUN_MOD \n")
return(NULL)
}
if ((identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) &
length(IniResLevels) != 3) {
stop("The length of IniStates must be 3 for the chosen FUN_MOD \n")
return(NULL)
}
} else if (is.null(IniStates)) {
if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
IniResLevels <- as.double(c(0.3, 0.5, 0))
} else {
IniResLevels <- as.double(c(0.3, 0.5, NA))
}
}
} else {
if (!is.null(IniResLevels)) {
stop("IniResLevels can only be used with monthly or daily or hourly GR models \n")
}
}
## check IniStates
if (is.null(IniStates) & is.null(IniResLevels) & verbose) {
warning("\t Model states initialisation not defined -> default configuration used \n")
}
if (!is.null(IniStates) & !is.null(IniResLevels) & verbose) {
warning("\t IniStates and IniResLevels are both defined -> Store levels are taken from IniResLevels \n")
}
if("CemaNeige" %in% ObjectClass){ NLayers <- length(InputsModel$LayerPrecip); } else { NLayers <- 0; }
NState <- NULL;
if("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass){
if("hourly" %in% ObjectClass){ NState <- 3*24*20 + 7; }
if("daily" %in% ObjectClass){ if (identical(FUN_MOD,RunModel_GR5J)){NState <- 2*20 + 2*NLayers + 7; } else {NState <- 3*20 + 2*NLayers + 7; }}
if("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass){
if("hourly" %in% ObjectClass){ NState <- 7 + 3*24*20 }
if("daily" %in% ObjectClass){ NState <- 7 + 3*20 + 2*NLayers }
if("monthly" %in% ObjectClass){ NState <- 2; }
if("yearly" %in% ObjectClass){ NState <- 1; }
}
if(!is.null(IniStates)){
if(!is.vector( IniStates) ){ stop("IniStates must be a vector of numeric values \n"); return(NULL); }
if(!is.numeric(IniStates) ){ stop("IniStates must be a vector of numeric values \n"); return(NULL); }
if(length(IniStates)!=NState){ stop(paste("The length of IniStates must be ",NState," for the chosen FUN_MOD \n",sep="")); return(NULL); }
} else {
IniStates <- as.double(rep(0.0,NState));
}
if("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)){
if(!is.null(IniResLevels)){
if(!is.vector(IniResLevels) ){ stop("IniResLevels must be a vector of numeric values \n"); return(NULL); }
if(!is.numeric(IniResLevels)){ stop("IniResLevels must be a vector of numeric values \n"); return(NULL); }
if(length(IniResLevels)!=2 ) { stop("The length of IniStates must be 2 for the chosen FUN_MOD \n"); return(NULL); }
} else {
IniResLevels <- as.double(c(0.3,0.5));
if (!is.null(IniStates)) {
if (!inherits(IniStates, "IniStates")) {
stop("IniStates must be an object of class IniStates\n")
return(NULL)
}
# print(str(IniStates))
# print(ObjectClass)
# print(class(IniStates))
# print(sum(ObjectClass %in% class(IniStates)))
if (sum(ObjectClass %in% class(IniStates)) < 2) {
stop(paste0("Non convenient IniStates for this FUN_MOD\n"))
return(NULL)
}
if (identical(FUN_MOD, RunModel_GR1A) & !is.null(IniStates)) { ## GR1A
stop(paste0("IniStates is not available for this FUN_MOD\n"))
return(NULL)
}
if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !all(is.na(IniStates$UH$UH1))) { ## GR5J
stop(paste0("Non convenient IniStates for this FUN_MOD. In IniStates, UH1 have to be a vector of NA for GR5J \n"))
return(NULL)
}
if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & is.na(IniStates$Store$Exp)) { ## GR6J
stop(paste0("Non convenient IniStates for this FUN_MOD. GR6J need an exponential store value in IniStates\n"))
return(NULL)
}
if (!(identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & !is.na(IniStates$Store$Exp)) { ## except GR6J
stop(paste0("Non convenient IniStates for this FUN_MOD\n"))
return(NULL)
}
# if (length(na.omit(unlist(IniStates))) != NState) {
# stop(paste0("The length of IniStates must be ", NState, " for the chosen FUN_MOD \n"))
# return(NULL)
# }
if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$G ))) {
IniStates$CemaNeigeLayers$G <- NULL
}
if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$eTG))) {
IniStates$CemaNeigeLayers$eTG <- NULL
}
IniStates$Store$Rest <- rep(NA, 4)
IniStates <- unlist(IniStates)
IniStates[is.na(IniStates)] <- 0
# print(IniStates)
if ("monthly" %in% ObjectClass) {
IniStates <- IniStates[seq_len(NState)]
# print(NState)
# print(IniStates)
}
} else {
if(!is.null(IniResLevels)){ stop("IniResLevels can only be used with monthly or daily or hourly GR models \n") }
IniStates <- as.double(rep(0.0, NState))
}
# if (is.null(IniStates) & is.null(IniResLevels) & verbose) {
# warning("\t Model states initialisation not defined -> default configuration used \n")
# }
# if("CemaNeige" %in% ObjectClass){ NLayers <- length(InputsModel$LayerPrecip); } else { NLayers <- 0; }
# NState <- NULL;
# if("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass){
# if("hourly" %in% ObjectClass){ NState <- 3*24*20 + 7; }
# if("daily" %in% ObjectClass){ if (identical(FUN_MOD,RunModel_GR5J)){NState <- 2*20 + 2*NLayers + 7; } else {NState <- 3*20 + 2*NLayers + 7; }}
# if("monthly" %in% ObjectClass){ NState <- 2; }
# if("yearly" %in% ObjectClass){ NState <- 1; }
# }
# if (!is.null(IniStates)) {
# if (!is.vector( IniStates) | !is.numeric(IniStates)) {
# stop("IniStates must be a vector of numeric values \n")
# return(NULL)
# }
# if (length(IniStates) != NState) {
# stop(paste0("The length of IniStates must be ", NState, " for the chosen FUN_MOD \n"))
# return(NULL)
# }
# } else {
# IniStates <- as.double(rep(0.0, NState))
# IniStates[1:3] <- NA
# }
# if ("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)) {
# if (!is.null(IniResLevels)) {
# if (!is.vector(IniResLevels) | !is.numeric(IniResLevels)) {
# stop("IniResLevels must be a vector of numeric values \n")
# return(NULL)
# }
# if ((identical(FUN_MOD, RunModel_GR4H) |
# identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
# identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
# identical(FUN_MOD, RunModel_GR2M)) &
# length(IniResLevels) != 2) {
# stop("The length of IniStates must be 2 for the chosen FUN_MOD \n")
# return(NULL)
# }
# if ((identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) &
# length(IniResLevels) != 3) {
# stop("The length of IniStates must be 3 for the chosen FUN_MOD \n")
# return(NULL)
# }
# } else {
# IniResLevels <- as.double(c(0.3, 0.5, 0))
# }
# } else {
# if (!is.null(IniResLevels)) {
# stop("IniResLevels can only be used with monthly or daily or hourly GR models \n")
# }
# }
##check_Outputs_Cal_and_Sim
......
......@@ -54,7 +54,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
......@@ -74,6 +74,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
)
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]);
......@@ -97,7 +98,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
} else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if("IniResLevels" %in% names(RunOptions)){
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)
}
......@@ -120,6 +121,13 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
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]],
verbose = FALSE)
if(RunOptions$RunSnowModule & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
##Output_data_preparation
......@@ -138,14 +146,14 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers),
list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputsMod[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(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); }
##End
......
......@@ -54,7 +54,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
......@@ -97,7 +97,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
} else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if("IniResLevels" %in% names(RunOptions)){
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)
}
......@@ -120,6 +120,13 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
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]],
verbose = FALSE)
if(RunOptions$RunSnowModule & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
##Output_data_preparation
......@@ -138,14 +145,14 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(CemaNeigeLayers),
list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputsMod[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(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd"); }
##End
......
......@@ -59,7 +59,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
##Call_DLL_CemaNeige_________________________
for(iLayer in 1:NLayers){
StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
......@@ -102,9 +102,10 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
} else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim); }
##Use_of_IniResLevels
if("IniResLevels" %in% names(RunOptions)){
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(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)