RunModel_GR4J.R 9.86 KB
Newer Older
Delaigue Olivier's avatar
Delaigue Olivier committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#*****************************************************************************************************************
#' Function which performs a single run for the GR4J daily lumped model.
#'
#' For further details on the model, see the references section.
#' For further details on the argument structures and initialisation options, see \code{\link{CreateRunOptions}}.
#*****************************************************************************************************************
#' @title Run with the GR4J hydrological model
#' @author Laurent Coron (December 2013)
#' @references
#'   Perrin, C., C. Michel and V. Andréassian (2003), 
#'       Improvement of a parsimonious model for streamflow simulation, 
#'       Journal of Hydrology, 279(1-4), 275-289, doi:10.1016/S0022-1694(03)00225-7.
#' @seealso \code{\link{RunModel_GR5J}}, \code{\link{RunModel_GR6J}}, \code{\link{RunModel_CemaNeigeGR4J}},
#'          \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}.
#' @example tests/example_RunModel_GR4J.R
#' @useDynLib airgr
#' @encoding UTF-8
#' @export
#_FunctionInputs__________________________________________________________________________________________________
#' @param  InputsModel         [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details
#' @param  RunOptions          [object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details
#' @param  Param               [numeric] vector of 4 parameters                                                             
#'                             \tabular{ll}{                                                                      
#'                             GR4J X1      \tab production store capacity [mm]                                \cr
#'                             GR4J X2      \tab intercatchment exchange coefficient [mm/d]                    \cr
#'                             GR4J X3      \tab routing store capacity [mm]                                   \cr
#'                             GR4J X4      \tab unit hydrograph time constant [d]                             \cr
#'                             }                                                                                  
#_FunctionOutputs_________________________________________________________________________________________________
#' @return  [list] list containing the function outputs organised as follows:                                         
#'          \tabular{ll}{                                                                                         
#'          \emph{$DatesR  }          \tab [POSIXlt] series of dates                                                     \cr
#'          \emph{$PotEvap }          \tab [numeric] series of input potential evapotranspiration [mm/d]                 \cr
#'          \emph{$Precip  }          \tab [numeric] series of input total precipitation [mm/d]                          \cr
#'          \emph{$Prod    }          \tab [numeric] series of production store level (X(2)) [mm]                        \cr
#'          \emph{$AE      }          \tab [numeric] series of actual evapotranspiration [mm/d]                          \cr
#'          \emph{$Perc    }          \tab [numeric] series of percolation (PERC) [mm/d]                                 \cr
#'          \emph{$PR      }          \tab [numeric] series of PR=PN-PS+PERC [mm/d]                                      \cr
#'          \emph{$Q9      }          \tab [numeric] series of HU1 outflow (Q9) [mm/d]                                   \cr
#'          \emph{$Q1      }          \tab [numeric] series of HU2 outflow (Q1) [mm/d]                                   \cr
#'          \emph{$Rout    }          \tab [numeric] series of routing store level (X(1)) [mm]                           \cr
#'          \emph{$Exch    }          \tab [numeric] series of potential semi-exchange between catchments [mm/d]         \cr
#'          \emph{$AExch   }          \tab [numeric] series of actual exchange between catchments (1+2) [mm/d]           \cr
#'          \emph{$QR      }          \tab [numeric] series of routing store outflow (QR) [mm/d]                         \cr
#'          \emph{$QD      }          \tab [numeric] series of direct flow from HU2 after exchange (QD) [mm/d]           \cr
#'          \emph{$Qsim    }          \tab [numeric] series of Qsim [mm/d]                                               \cr
#'          \emph{$StateEnd}          \tab [numeric] states at the end of the run (res. levels, HU1 levels, HU2 levels) [mm] \cr
#'          }                                                                                                     
#'          (refer to the provided references or to the package source code for further details on these model outputs)
#*****************************************************************************************************************'
RunModel_GR4J <- function(InputsModel,RunOptions,Param){

    NParam <- 4;
    FortranOutputs <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");

    ##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,"GR"         )==FALSE){ stop("InputsModel must be of class 'GR'          \n"); return(NULL); }  
      if(inherits(RunOptions,"RunOptions"  )==FALSE){ stop("RunOptions must be of class 'RunOptions'   \n"); return(NULL); }  
      if(inherits(RunOptions,"GR"          )==FALSE){ stop("RunOptions must be of class 'GR'           \n"); return(NULL); }  
      if(!is.vector(Param)){ stop("Param must be a 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,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);  }

    ##Use_of_IniResLevels
      if("IniResLevels" %in% names(RunOptions)){
        RunOptions$IniStates[1] <- RunOptions$IniResLevels[2]*Param[3];  ### routing store level (mm)
        RunOptions$IniStates[2] <- RunOptions$IniResLevels[1]*Param[1];  ### production store level (mm)
      }

    ##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;

    ##Output_data_preparation
      IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
      ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim;
      ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
      ##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","daily","GR");
      return(OutputsModel);

}