RunModel.R

RunModel <- function (InputsModel, RunOptions, Param, FUN_MOD) {
  
  FUN_MOD <- match.fun(FUN_MOD)
  
  if (inherits(InputsModel, "SD")) {
    OutputsModelDown <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions, 
                                Param = Param[-length(Param)])
    OutputsModelDown$QsimDown <- OutputsModelDown$Qsim
    
    if (inherits(InputsModel, "daily")) {
      TimeStep <- 60 * 60 * 24
    }
    if (inherits(InputsModel, "hourly")) {
      TimeStep <- 60 * 60
    }
    
    # total area 
    AreaTot <- sum(InputsModel$BasinAreas)
    
    # propagation time from upstream meshes to outlet
    PT <- InputsModel$LengthHydro[1, ] / Param[length(Param)] / TimeStep
    HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
    
    NbUpBasins <- length(InputsModel$LengthHydro)
    LengthTs <- length(OutputsModelDown$QsimDown)
    OutputsModelDown$Qsim <- OutputsModelDown$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] / AreaTot
    
    for (upstream_basin in seq_len(NbUpBasins)) {
      QobsUpstr <- InputsModel$QobsUpstr[RunOptions$IndPeriod_Run, upstream_basin]
      OutputsModelDown$Qsim <- OutputsModelDown$Qsim + 
        c(rep(0, floor(PT[upstream_basin])), 
          QobsUpstr[(1 + floor(PT[upstream_basin])):LengthTs]) *
        HUTRANS[1, upstream_basin] * InputsModel$BasinAreas[upstream_basin] / AreaTot + 
        c(rep(0, floor(PT[upstream_basin] + 1)), 
          QobsUpstr[(2 + floor(PT[upstream_basin])):LengthTs]) *
        HUTRANS[2, upstream_basin] * InputsModel$BasinAreas[upstream_basin] / AreaTot
    }
    
  } else {
    
    OutputsModelDown <- FUN_MOD(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
    
  }
  return(OutputsModelDown)
}