From da5bb03ca4b8e4849430ec88dbe48f66a91e8a7a Mon Sep 17 00:00:00 2001 From: Delaigue Olivier <olivier.delaigue@irstea.priv> Date: Wed, 30 Jan 2019 13:48:52 +0100 Subject: [PATCH] v1.1.2.38 CLEAN: RunModel_CemaNeige syntax cleaned --- DESCRIPTION | 2 +- NEWS.rmd | 2 +- R/RunModel_CemaNeige.R | 243 +++++++++++++++++++++++++---------------- 3 files changed, 150 insertions(+), 97 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 35f8adad..7921b101 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: airGR Type: Package Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling -Version: 1.1.2.37 +Version: 1.1.2.38 Date: 2019-01-29 Authors@R: c( person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")), diff --git a/NEWS.rmd b/NEWS.rmd index a0f4f31a..94c0ac08 100644 --- a/NEWS.rmd +++ b/NEWS.rmd @@ -13,7 +13,7 @@ output: -### 1.1.2.37 Release Notes (2019-01-29) +### 1.1.2.38 Release Notes (2019-01-29) diff --git a/R/RunModel_CemaNeige.R b/R/RunModel_CemaNeige.R index 775377fb..52e3789e 100644 --- a/R/RunModel_CemaNeige.R +++ b/R/RunModel_CemaNeige.R @@ -1,97 +1,150 @@ -RunModel_CemaNeige <- function(InputsModel,RunOptions,Param){ - - NParam <- 2 - FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp"); - - ##Arguments_check - if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); } - if(inherits(InputsModel,"daily" )==FALSE){ stop("InputsModel must be of class 'daily' \n"); return(NULL); } - if(inherits(InputsModel,"CemaNeige" )==FALSE){ stop("InputsModel must be of class 'CemaNeige' \n"); return(NULL); } - if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("RunOptions must be of class 'RunOptions' \n"); return(NULL); } - if(inherits(RunOptions,"CemaNeige" )==FALSE){ stop("RunOptions must be of class 'CemaNeige' \n"); return(NULL); } - if(!is.vector(Param) | !is.numeric(Param)){ stop("Param must be a numeric vector \n"); return(NULL); } - if(sum(!is.na(Param))!=NParam){ stop(paste("Param must be a vector of length ",NParam," and contain no NA \n",sep="")); return(NULL); } - Param <- as.double(Param); - - ##Input_data_preparation - if(identical(RunOptions$IndPeriod_WarmUp,0)){ RunOptions$IndPeriod_WarmUp <- NULL; } - IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run); - IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):length(IndPeriod1); - ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; - ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; - - - - ##SNOW_MODULE________________________________________________________________________________## - ParamCemaNeige <- Param; - NLayers <- length(InputsModel$LayerPrecip); - if(sum(is.na(ParamCemaNeige))!=0){ stop("Param contains missing values \n"); return(NULL); } - 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)] - 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] - 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 - StateStart=StateStartCemaNeige, ### state variables used when the model run starts - NOutputs=as.integer(length(IndOutputsCemaNeige)), ### number of output series - IndOutputs=IndOutputsCemaNeige, ### indices of output series - ##outputs - Outputs=matrix(as.double(-999.999),nrow=length(IndPeriod1),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) - ) - 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]; - if(ExportStateEnd){ CemaNeigeStateEnd <- c(CemaNeigeStateEnd,RESULTS$StateEnd); } - rm(RESULTS); - } ###ENDFOR_iLayer - names(CemaNeigeLayers) <- paste("Layer",formatC(1:NLayers,width=2,flag="0"),sep=""); - - if (ExportStateEnd) { - 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]], - verbose = FALSE) - } - - ##Output_data_preparation - if(ExportDatesR==FALSE & ExportStateEnd==FALSE){ - OutputsModel <- list(CemaNeigeLayers); - names(OutputsModel) <- NameCemaNeigeLayers; } - if(ExportDatesR==TRUE & ExportStateEnd==FALSE){ - OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), - list(CemaNeigeLayers) ); - names(OutputsModel) <- c("DatesR",NameCemaNeigeLayers); } - if(ExportDatesR==FALSE & ExportStateEnd==TRUE){ - OutputsModel <- c( list(CemaNeigeLayers), - list(CemaNeigeStateEnd)); - names(OutputsModel) <- c(NameCemaNeigeLayers,"StateEnd"); } - if(ExportDatesR==TRUE & ExportStateEnd==TRUE){ - OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), - list(CemaNeigeLayers), - list(CemaNeigeStateEnd)); - names(OutputsModel) <- c("DatesR",NameCemaNeigeLayers,"StateEnd"); } - - ##End - class(OutputsModel) <- c("OutputsModel","daily","CemaNeige"); - return(OutputsModel); +RunModel_CemaNeige <- function(InputsModel, RunOptions, Param) { + + + NParam <- 2L + FortranOutputsCemaNeige <- c("Pliq", "Psol", + "SnowPack", "ThermalState", "Gratio", + "PotMelt", "Melt", "PliqAndMelt", "Temp") + + + + ## Arguments_check + if (!inherits(InputsModel, "InputsModel")) { + stop("'InputsModel' must be of class 'InputsModel'") + return(NULL) + } + if (!inherits(InputsModel, "daily")) { + stop("'InputsModel' must be of class 'daily'") + return(NULL) + } + if (!inherits(InputsModel, "CemaNeige")) { + stop("'InputsModel' must be of class 'CemaNeige'") + return(NULL) + } + if (!inherits(RunOptions, "RunOptions")) { + stop("'RunOptions' must be of class 'RunOptions'") + return(NULL) + } + if (!inherits(RunOptions, "CemaNeige")) { + stop("'RunOptions' must be of class 'CemaNeige'") + return(NULL) + } + if (!is.vector(Param) | !is.numeric(Param)) { + stop("'Param' must be a numeric vector") + return(NULL) + } + if (sum(!is.na(Param)) != NParam) { + stop(sprintf("'Param' must be a vector of length %i and contain no NA", NParam)) + return(NULL) + } + + ## Input_data_preparation + if (identical(RunOptions$IndPeriod_WarmUp, 0)) { + RunOptions$IndPeriod_WarmUp <- NULL + } + IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run) + IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp) + 1):length(IndPeriod1) + ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim + ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim + -} + + + ## SNOW_MODULE________________________________________________________________________________ + ParamCemaNeige <- Param + NLayers <- length(InputsModel$LayerPrecip) + + if (sum(is.na(ParamCemaNeige)) != 0) { + stop("Param contains missing values") + return(NULL) + } + if ("all" %in% RunOptions$Outputs_Sim) { + IndOutputsCemaNeige <- 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)] + 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] + 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 = as.double(ParamCemaNeige), ### parameter set + NStates = as.integer(2), ### number of state variables used for model initialising = 2 + StateStart = StateStartCemaNeige, ### state variables used when the model run starts + 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) + ) + 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] + if (ExportStateEnd) { + CemaNeigeStateEnd <- c(CemaNeigeStateEnd, RESULTS$StateEnd) + } + rm(RESULTS) + + } ### ENDFOR_iLayer + + names(CemaNeigeLayers) <- paste("Layer", formatC(1:NLayers, width = 2, flag = "0"), sep = "") + + if (ExportStateEnd) { + 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]], + verbose = FALSE) + } + + ## Output_data_preparation + if (!ExportDatesR & !ExportStateEnd) { + OutputsModel <- list(CemaNeigeLayers) + names(OutputsModel) <- NameCemaNeigeLayers + } + if (ExportDatesR & !ExportStateEnd) { + OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), + list(CemaNeigeLayers)) + names(OutputsModel) <- c("DatesR", NameCemaNeigeLayers) + } + if (!ExportDatesR & ExportStateEnd) { + OutputsModel <- c(list(CemaNeigeLayers), + list(CemaNeigeStateEnd)) + names(OutputsModel) <- c(NameCemaNeigeLayers, "StateEnd") + } + if (ExportDatesR & ExportStateEnd) { + OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]), + list(CemaNeigeLayers), + list(CemaNeigeStateEnd)) + names(OutputsModel) <- c("DatesR", NameCemaNeigeLayers, "StateEnd") + } + + + ## End + class(OutputsModel) <- c("OutputsModel", "daily", "CemaNeige") + return(OutputsModel) +} + -- GitLab