From 5702ffd16d0efe32ba6d4ef13613c4c01b9b051f Mon Sep 17 00:00:00 2001 From: Delaigue Olivier <olivier.delaigue@irstea.priv> Date: Wed, 4 Dec 2019 15:37:40 +0100 Subject: [PATCH] v1.4.0.0 NEW: add a RunModel_GR5H fun #13 --- DESCRIPTION | 4 +- NEWS.md | 2 +- R/RunModel_GR5H.R | 113 ++++++++++++++ src/frun_GR5H.f | 372 ++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 488 insertions(+), 3 deletions(-) create mode 100644 R/RunModel_GR5H.R create mode 100644 src/frun_GR5H.f diff --git a/DESCRIPTION b/DESCRIPTION index 61be0e69..793c7f85 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.3.2.69 -Date: 2019-12-02 +Version: 1.4.0.0 +Date: 2019-12-04 Authors@R: c( person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")), person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@irstea.fr"), diff --git a/NEWS.md b/NEWS.md index 08bd8a4c..23a55f3e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,7 +2,7 @@ -### 1.3.2.69 Release Notes (2019-12-02) +### 1.4.0.0 Release Notes (2019-12-04) #### New features diff --git a/R/RunModel_GR5H.R b/R/RunModel_GR5H.R new file mode 100644 index 00000000..65d4cc4c --- /dev/null +++ b/R/RunModel_GR5H.R @@ -0,0 +1,113 @@ +RunModel_GR5H <- function(InputsModel,RunOptions,Param){ + + NParam <- 5; + FortranOutputs <- .FortranOutputs(GR = "GR5H")$GR + IsIntStore <- inherits(RunOptions, "interception") + + ##Arguments_check + if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") } + if(inherits(InputsModel,"hourly" )==FALSE){ stop("'InputsModel' must be of class 'hourly' ") } + if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") } + if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") } + if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") } + if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") } + if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) } + Param <- as.double(Param); + + Param_X1X3_threshold <- 1e-2 + Param_X4_threshold <- 0.5 + if (Param[1L] < Param_X1X3_threshold) { + warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold)) + Param[1L] <- Param_X1X3_threshold + } + if (Param[3L] < Param_X1X3_threshold) { + warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold)) + Param[3L] <- Param_X1X3_threshold + } + if (Param[4L] < Param_X4_threshold) { + warning(sprintf("Param[4] (X4: unit hydrograph time constant [h]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold)) + Param[4L] <- Param_X4_threshold + } + + ##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); } + + ##Output_data_preparation + IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries; + ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim; + ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim; + + ##Use_of_IniResLevels + if(!is.null(RunOptions$IniResLevels)){ + RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm) + RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm) + if(IsIntStore){ + RunOptions$IniStates[4] <- RunOptions$IniResLevels[4]*RunOptions$Imax; ### interception store level (mm) + } + } + + ##Call_fortan + RESULTS <- .Fortran("frun_gr5h",PACKAGE="airGR", + ##inputs + LInputs=LInputSeries, ### length of input and output series + InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/h] + InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/h] + 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 + IsIntStore=as.integer(IsIntStore), ### use of interception store + Imax=RunOptions$Imax, ### maximal capacity of interception store + 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 or mm/h] + 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; + if (ExportStateEnd) { + RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location + RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel, + ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL, + IntStore = RESULTS$StateEnd[4L], + UH1 = NULL, UH2 = RESULTS$StateEnd[(1:(40*24))+(7+20*24)], + GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL, + verbose = FALSE) + } + + ##Output_data_preparation + ##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_StateEnd_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_StateEnd + 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","hourly","GR"); + if(IsIntStore) { + class(OutputsModel) <- c(class(OutputsModel), "interception") + } + return(OutputsModel); + +} diff --git a/src/frun_GR5H.f b/src/frun_GR5H.f new file mode 100644 index 00000000..8728e49b --- /dev/null +++ b/src/frun_GR5H.f @@ -0,0 +1,372 @@ +!------------------------------------------------------------------------------ +! Subroutines relative to the annual GR5H model +!------------------------------------------------------------------------------ +! TITLE : airGR +! PROJECT : airGR +! FILE : frun_GR5H.f +!------------------------------------------------------------------------------ +! AUTHORS +! Original code: N. Le Moine, A. Ficchì +! Cleaning and formatting for airGR: L. Coron +! Further cleaning: G. Thirel +!------------------------------------------------------------------------------ +! Creation date: 2006 +! Last modified: 26/11/2019 +!------------------------------------------------------------------------------ +! REFERENCES +! Ficchì, A., Perrin, C., and Andréassian, V. (2019). Hydrological modelling at +! multiple sub-daily time steps: model improvement via flux-matching, Journal +! of Hydrology, 575, 1308-1327, https://doi.org/10.1016/j.jhydrol.2019.05.084. +! +! Le Moine, N. (2008). Le bassin versant de surface vu par le souterrain : une +! voie d'amélioration des performances et du réalisme des modèles +! pluie-débit ? PhD thesis (French), UPMC, Paris, France. +! +! Pushpalatha, R., C. Perrin, N. Le Moine, T. Mathevet, and V. Andréassian +! (2011). A downward structural sensitivity analysis of hydrological models to +! improve low-flow simulation. Journal of Hydrology, 411(1-2), 66-76. +! doi:10.1016/j.jhydrol.2011.09.034. +!------------------------------------------------------------------------------ +! Quick description of public procedures: +! 1. frun_gr5h +! 2. MOD_GR5H +!------------------------------------------------------------------------------ + + + SUBROUTINE frun_gr5h(LInputs,InputsPrecip,InputsPE,NParam,Param, + & NStates,StateStart,IsIntStore,Imax,NOutputs,IndOutputs, + & Outputs,StateEnd) +! Subroutine that initializes GR5H, get its parameters, performs the call +! to the MOD_GR5H subroutine at each time step, and stores the final states +! Inputs +! LInputs ! Integer, length of input and output series +! InputsPrecip ! Vector of real, input series of total precipitation [mm/hour] +! InputsPE ! Vector of real, input series of potential evapotranspiration (PE) [mm/hour] +! NParam ! Integer, number of model parameters +! Param ! Vector of real, parameter set +! NStates ! Integer, number of state variables +! StateStart ! Vector of real, state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm]) +! IsIntStore ! Integer, whether we should use the interception store [1] or not [0] +! Imax ! Real, fixed capacity of the interception store [mm] (used only if IsIntStore=1) +! NOutputs ! Integer, number of output series +! IndOutputs ! Vector of integer, indices of output series +! Outputs +! Outputs ! Vector of real, output series +! StateEnd ! Vector of real, state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm]) + + + !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr5h + + + Implicit None + + !! dummies + ! in + integer, intent(in) :: LInputs,NParam,NStates,NOutputs + doubleprecision, dimension(LInputs), intent(in) :: InputsPrecip + doubleprecision, dimension(LInputs), intent(in) :: InputsPE + doubleprecision, dimension(NParam), intent(in) :: Param + doubleprecision, dimension(NStates), intent(in) :: StateStart + doubleprecision, intent(in) :: Imax ! interception capacity (fixed parameter) used only if IsIntStore=1 + integer, intent(in) :: IsIntStore ! 1 if interception store is used, 0 otherwise + integer, dimension(NOutputs), intent(in) :: IndOutputs + ! out + doubleprecision, dimension(NStates), intent(out) :: StateEnd + doubleprecision, dimension(LInputs,NOutputs), + & intent(out) :: Outputs + + !! locals + logical :: IsIntStoreBool ! TRUE if interception store is used, FALSE otherwise + integer :: I,K + integer, parameter :: NH=480,NMISC=30 + doubleprecision, dimension(3) :: St + doubleprecision, dimension(2*NH) :: StUH2, OrdUH2 + doubleprecision, dimension(NMISC) :: MISC + doubleprecision :: D,P1,E,Q + + IF (IsIntStore .EQ. 1) IsIntStoreBool = .TRUE. + IF (IsIntStore .EQ. 0) IsIntStoreBool = .FALSE. + + !-------------------------------------------------------------- + ! Initializations + !-------------------------------------------------------------- + + ! initialization of model states to zero + St=0. + StUH2=0. + + ! initialization of model states using StateStart + St(1) = StateStart(1) + St(2) = StateStart(2) + IF (IsIntStoreBool .EQ. .TRUE.) St(3) = StateStart(4) + + DO I=1,2*NH + StUH2(I)=StateStart(7+NH+I) + ENDDO + + ! parameter values + ! Param(1) : production store capacity (X1 - PROD) [mm] + ! Param(2) : intercatchment exchange coefficient (X2 - CES1) [mm/h] + ! Param(3) : routing store capacity (X3 - ROUT) [mm] + ! Param(4) : time constant of unit hydrograph (X4 - TB) [hour] + ! Param(5) : intercatchment exchange threshold (X5 - CES2) [-] + + !computation of HU ordinates + OrdUH2 = 0. + + D=1.25 + CALL UH2_H(OrdUH2,Param(4),D) + + ! initialization of model outputs + Q = -999.999 + MISC = -999.999 +! StateEnd = -999.999 !initialization made in R +! Outputs = -999.999 !initialization made in R + + + + !-------------------------------------------------------------- + ! Time loop + !-------------------------------------------------------------- + DO k=1,LInputs + P1=InputsPrecip(k) + E =InputsPE(k) +! Q = -999.999 +! MISC = -999.999 + ! model run on one time step + CALL MOD_GR5H(St,StUH2,OrdUH2,Param,IsIntStoreBool,Imax,P1,E,Q,MISC) + ! storage of outputs + DO I=1,NOutputs + Outputs(k,I)=MISC(IndOutputs(I)) + ENDDO + ENDDO + ! model states at the end of the run + StateEnd(1) = St(1) + StateEnd(2) = St(2) + StateEnd(4) = St(3) + DO K=1,2*NH + StateEnd(7+NH+K)=StUH2(K) + ENDDO + + RETURN + + ENDSUBROUTINE + + + + + +!################################################################################################################################ + + + + +!********************************************************************** + SUBROUTINE MOD_GR5H(St,StUH2,OrdUH2,Param,IsIntStoreBool,Imax,P1,E,Q, + &MISC) +! Calculation of streamflow on a single time step (hour) with the GR5H model +! Inputs: +! St Vector of real, model states in stores at the beginning of the time step [mm] +! StUH2 Vector of real, model states in Unit Hydrograph 2 at the beginning of the time step [mm/h] +! OrdUH2 Vector of real, ordinates in UH2 [-] +! Param Vector of real, model parameters [various units] +! IsIntStoreBool Logical, whether interception store is used +! Imax Real, value of interception store capacity, fixed parameter [mm] +! P1 Real, value of rainfall during the time step [mm/h] +! E Real, value of potential evapotranspiration during the time step [mm/h] +! Outputs: +! St Vector of real, model states in stores at the end of the time step [mm/h] +! StUH2 Vector of real, model states in Unit Hydrograph 2 at the end of the time step [mm/h] +! Q Real, value of simulated flow at the catchment outlet for the time step [mm/h] +! MISC Vector of real, model outputs for the time step [mm or mm/h] +!********************************************************************** + Implicit None + + !! locals + integer, parameter :: NParam=5,NMISC=30,NH=480 + doubleprecision :: A,EN,ER,PN,PR,PS,WS,tanHyp,EI + doubleprecision :: PERC,EXCH,QR,QD,Q1,Q9 + doubleprecision :: AE,AEXCH1,AEXCH2 + integer :: K + doubleprecision, parameter :: B=0.9 + doubleprecision, parameter :: stored_val=759.69140625 + doubleprecision :: expWS, TWS, Sr, Rr ! speed-up + + !! dummies + ! in + doubleprecision, dimension(NParam), intent(in) :: Param + logical, intent(in) :: IsIntStoreBool + doubleprecision, intent(in) :: P1,E,Imax + doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2 + ! inout + doubleprecision, dimension(3), intent(inout) :: St + doubleprecision, dimension(2*NH), intent(inout) :: StUH2 + ! out + doubleprecision, intent(out) :: Q + doubleprecision, dimension(NMISC), intent(out) :: MISC + + A=Param(1) + + +! Interception and production store + IF (IsIntStore .EQ. .TRUE.) THEN + + ! MODIFIED - A. Ficchi + ! Calculation of interception fluxes [EI] and throughfall [PTH] + ! & update of the Interception store state, St(3) + + ! Interception store calculation, with evaporation prior to throughfall + EI=MIN(E, P1+St(3)) + PN=MAX(0., P1-(Imax-St(3))-EI) + St(3)=St(3)+P1-EI-PN + EN=MAX(0., E-EI) + + ! Production (SMA) store, saving the total actual evaporation including evaporation from interception store (EI) + IF(EN.GT.0) THEN + WS=EN/A + IF(WS.GT.13)WS=13. + expWS=exp(2.*WS) + TWS = (expWS - 1.)/(expWS + 1.) + Sr = St(1)/A + ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS) + St(1)=St(1)-ER + AE=ER+EI + ELSE + AE=EI + ENDIF + + IF(PN.GT.0.) THEN + WS=PN/A + IF(WS.GT.13) WS=13. + expWS=exp(2.*WS) + TWS = (expWS - 1.)/(expWS + 1.) + Sr = St(1)/A + PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS) + PR=PN-PS + St(1)=St(1)+PS + ELSE + PS=0. + PR=0. + ENDIF + ENDIF + + + IF (IsIntStore .EQ. .FALSE.) THEN + IF(P1.LE.E) THEN + EN=E-P1 + PN=0. + WS=EN/A + IF(WS.GT.13.) WS=13. + ! speed-up + expWS=exp(2.*WS) + TWS = (expWS - 1.)/(expWS + 1.) + Sr = St(1)/A + ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS) + ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS)) + ! end speed-up + AE=ER+P1 + EI = P1 + St(1)=St(1)-ER + PS=0. + PR=0. + ELSE + EN=0. + AE=E + EI = 0. + PN=P1-E + WS=PN/A + IF(WS.GT.13.) WS=13. + ! speed-up + expWS=exp(2.*WS) + TWS = (expWS - 1.)/(expWS + 1.) + Sr = St(1)/A + PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS) + ! PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS)) + ! end speed-up + + PR=PN-PS + St(1)=St(1)+PS + ENDIF + ENDIF + +! Percolation from production store + IF(St(1).LT.0.) St(1)=0. + + ! speed-up + ! (21/4)**4 = 759.69140625 = stored_val + Sr = St(1)/Param(1) + Sr = Sr * Sr + Sr = Sr * Sr + PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/stored_val))) + ! PERC=X(2)*(1.-(1.+(X(2)/(21./4.*Param(1)))**4.)**(-0.25)) + ! end speed-up + + St(1)=St(1)-PERC + + PR=PR+PERC + +! Convolution of unit hydrograph UH2 + DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1.))) + StUH2(K)=StUH2(K+1)+OrdUH2(K)*PR + ENDDO + StUH2(2*NH)=OrdUH2(2*NH)*PR + +! Split of unit hydrograph first component into the two routing components + Q9=StUH2(1)*B + Q1=StUH2(1)*(1.-B) + +! Potential intercatchment semi-exchange + EXCH=Param(2)*(St(2)/Param(3)-Param(5)) + +! Routing store + AEXCH1=EXCH + IF((St(2)+Q9+EXCH).LT.0.) AEXCH1=-St(2)-Q9 + St(2)=St(2)+Q9+EXCH + IF(St(2).LT.0.) St(2)=0. + + ! speed-up + Rr = St(2)/Param(3) + Rr = Rr * Rr + Rr = Rr * Rr + QR = St(2)*(1.-1./SQRT(SQRT(1.+Rr))) + ! QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.)) + ! end speed-up + + St(2)=St(2)-QR + +! Runoff from direct branch QD + AEXCH2=EXCH + IF((Q1+EXCH).LT.0.) AEXCH2=-Q1 + QD=MAX(0.d0,Q1+EXCH) + +! Total runoff + Q=QR+QD + IF(Q.LT.0.) Q=0. + +! Variables storage + MISC( 1)=E ! PE ! observed potential evapotranspiration [mm/h] + MISC( 2)=P1 ! Precip ! observed total precipitation [mm/h] + MISC( 3)=St(3) ! Interc ! interception store level (St(3)) [mm] + MISC( 4)=St(1) ! Prod ! production store level (St(1)) [mm] + MISC( 5)=PN ! Pn ! net rainfall [mm/h] + MISC( 6)=PS ! Ps ! part of Ps filling the production store [mm/h] + MISC( 7)=AE ! AE ! actual evapotranspiration [mm/h] + MISC( 8)=EI ! EI ! evapotranspiration from rainfall neutralisation or interception store [mm/h] + MISC( 9)=AE-EI ! AE-EI ! evapotranspiration from production store [mm/h] + MISC(10)=PERC ! Perc ! percolation (PERC) [mm/h] + MISC(11)=PR ! PR ! PR=PN-PS+PERC [mm/h] + MISC(12)=Q9 ! Q9 ! outflow from UH1 (Q9) [mm/h] + MISC(13)=Q1 ! Q1 ! outflow from UH2 (Q1) [mm/h] + MISC(14)=St(2) ! Rout ! routing store level (St(2)) [mm] + MISC(15)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/h] + MISC(16)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/h] + MISC(17)=AEXCH2 ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/h] + MISC(18)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/h] + MISC(19)=QR ! QR ! outflow from routing store (QR) [mm/h] + MISC(20)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/h] + MISC(21)=Q ! Qsim ! simulated outflow at catchment outlet [mm/h] + + + ENDSUBROUTINE + + -- GitLab