diff --git a/DESCRIPTION b/DESCRIPTION
index a136852130f107196cf0d7051ab96a8bebeb4075..61be0e691b72dc962135d329d25a13983c8510f7 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.59
-Date: 2019-11-20
+Version: 1.3.2.69
+Date: 2019-12-02
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 d6a5478552d99f2f2c8bcaae8bfc7504a5017f33..08bd8a4cc37eb86e976a3e45a55c2faf809e3ee6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
-### 1.3.2.59 Release Notes (2019-11-19)
+### 1.3.2.69 Release Notes (2019-12-02)
#### New features
@@ -15,10 +15,14 @@
- Fixed bug in <code>TransfoParam_GR1A()</code>. The number of model parameters was wrong (2 instead of 1) which caused an error when the GR1 model during the calibration.
- Fixed bug in <code>plot.OutputsModel()</code>. The function does not return any error message when <code>log_scale = TRUE</code>, <code>Qobs = NULL</code> and user want to draw flows time series.
- Fixed bug in <code>RunModel_*GR*()</code>. The functions do not return any error message anymore due to slightly negative values returned by GR4H, GR4J, GR5J or GR6J Fortran codes (the message was returned by <code>CreateIniStates()</code> when the final states were created). The <code>RunModel_*GR*()</code> functions now returns zero instead of these slightly negative values, except for the ExpStore where negatives values are allowed.
+- Fixed bug in the <code>.ErrorCrit()</code> function. The Box-Cox transformation formula is now corrected when the <code>ErrorCrit*()</code> are used.
+
#### Minor user-visible changes
-- Spurious flows set to <code>NA</code> into the <code>BasinObs</code> time series of the <code>L0123001</code> dataset.
+- Added the diagram of GR2M in the <code>RunModel_GR2M ()</code> documentation.
+- Fotran codes cleaned and translated from F77 to F90.
+
#### CRAN-compatibility updates
diff --git a/R/Utils.R b/R/Utils.R
index bce343222338fbbaf8ebdb9cfb77d5692665876f..435bced0e5c0d218bb70c102da3bb4ef53699c17 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -148,7 +148,7 @@
warning("zeroes detected in 'Qsim': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
}
}
- if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) {
+ if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon) & !(InputsCrit$transfo == "boxcox")) {
VarObs <- VarObs + InputsCrit$epsilon
VarSim <- VarSim + InputsCrit$epsilon
}
@@ -173,8 +173,9 @@
InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
}
if (InputsCrit$transfo == "boxcox") {
- VarSim <- (VarSim^0.25 - 0.01 * mean(VarSim, na.rm = TRUE)) / 0.25
- VarObs <- (VarObs^0.25 - 0.01 * mean(VarObs, na.rm = TRUE)) / 0.25
+ muTransfoVarObs <- (0.01 * mean(VarObs, na.rm = TRUE))^0.25
+ VarSim <- (VarSim^0.25 - muTransfoVarObs) / 0.25
+ VarObs <- (VarObs^0.25 - muTransfoVarObs) / 0.25
}
if (grepl("\\^", InputsCrit$transfo)) {
VarObs <- VarObs^transfoPow
@@ -221,7 +222,6 @@
CritCompute = CritCompute,
TS_ignore = TS_ignore,
Ind_TS_ignore = Ind_TS_ignore)
-
}
diff --git a/man/CreateInputsCrit.Rd b/man/CreateInputsCrit.Rd
index cbfb7e092c0fbda81bea45fbef948fad5f128c97..1ebdf5b7f108cc210a5bb675d8ba3f7be7fce63b 100644
--- a/man/CreateInputsCrit.Rd
+++ b/man/CreateInputsCrit.Rd
@@ -32,7 +32,7 @@ CreateInputsCrit(FUN_CRIT, InputsModel, RunOptions,
\item{BoolCrit}{(optional) [boolean (atomic or list)] boolean (the same length as \code{Obs}) giving the time steps to consider in the computation (all time steps are considered by default)}
-\item{transfo}{(optional) [character (atomic or list)] name of the transformation (e.g. \code{""}, \code{"sqrt"}, \code{"log"}, \code{"inv"}, \code{"sort"}, \code{"boxcox"} or a numeric value for power transformation (see details))}
+\item{transfo}{(optional) [character (atomic or list)] name of the transformation applied to the variables (e.g. \code{""}, \code{"sqrt"}, \code{"log"}, \code{"inv"}, \code{"sort"}, \code{"boxcox"} or a numeric value for power transformation (see details))}
\item{Weights}{(optional) [numeric (atomic or list)] vector of weights necessary to calculate a composite criterion (the same length as \code{FUN_CRIT}) giving the weights to use for elements of \code{FUN_CRIT} [-]. See details}
@@ -72,12 +72,29 @@ Creation of the \code{InputsCrit} object required to the \code{ErrorCrit_*} func
\details{
Users wanting to use \code{FUN_CRIT} functions that are not included in the package must create their own InputsCrit object accordingly. \cr \cr
-The syntax of the power transformation allows a numeric or a string of characters. For example for a squared transformation, the following can be used: \code{transfo = 2}, \code{transfo = "2"} or \code{transfo = "^2"}. Negative values are allowed. Fraction values are not allowed (e.g., \code{"-1/2"} must instead be written \code{"-0.5"}).
-\cr \cr
-In order to make sure that KGE and KGE2 remain dimensionless and not impacted by zero values, the Box-Cox transformation (\code{transfo = "boxcox"}) uses the formulation given in Equation 10 of Santos et al. (2018). Lambda is set to 0.25 accordingly.
-\cr \cr
-The epsilon value is useful when \code{"log"} or \code{"inv"} transformations are used (to avoid calculation of the inverse or of the logarithm of zero). The impact of this value and a recommendation about the epsilon value to use (usually one hundredth of average observation) are discussed in Pushpalatha et al. (2012) for NSE and in Santos et al. (2018) for KGE and KGE'. \cr \cr
+
+## ---- Transformations
+
+Transformations are simple functions applied to the observed and simulated variables used in order to change their distribution. Transformations are often used in hydrology for modifying the weight put on errors made for high flows or low flows. The following transformations are available: \cr \cr
+\itemize{
+ \item \code{""}: no transformation is used (default case)
+ \item \code{"sqrt"}: squared root transformation
+ \item \code{"log"}: logarithmic transformation (see below regarding the specific case of KGE or KGE2)
+ \item \code{"inv"}: inverse transformation
+ \item \code{"sort"}: sort transformation (the simulated and observed variables are sorted from lowest to highest)
+ \item \code{"boxcox"}: Box-Cox transformation (see below for details)
+ \item numeric: power transformation (see below for details)
+}
We do not advise computing KGE or KGE' with log-transformation as it might be wrongly influenced by discharge values close to 0 or 1 and the criterion value is dependent on the discharge unit. See Santos et al. (2018) for more details and alternative solutions (see the references list below). \cr \cr
+In order to make sure that KGE and KGE2 remain dimensionless and are not impacted by zero values, the Box-Cox transformation (\code{transfo = "boxcox"}) uses the formulation given in Equation 10 of Santos et al. (2018). Lambda is set to 0.25 accordingly. \cr \cr
+The syntax of the power transformation allows a numeric or a string of characters. For example for a squared transformation, the following can be used: \code{transfo = 2}, \code{transfo = "2"} or \code{transfo = "^2"}. Negative values are allowed. Fraction values are not allowed (e.g., \code{"-1/2"} must instead be written \code{"-0.5"}).\cr \cr
+
+## ---- The epsilon value
+
+The epsilon value is useful when \code{"log"} or \code{"inv"} transformations are used (to avoid calculation of the inverse or of the logarithm of zero). If an epsilon value is provided, then it is added to the observed and simulated variable time series at each time step and before the application of a transformation. The epsilon value has no effect when the \code{"boxcox"} transformation is used. The impact of this value and a recommendation about the epsilon value to use (usually one hundredth of average observation) are discussed in Pushpalatha et al. (2012) for NSE and in Santos et al. (2018) for KGE and KGE'. \cr \cr
+
+## ---- Single, multiple or composite criteria calculation
+
Users can set the following arguments as atomic or list: \code{FUN_CRIT}, \code{Obs}, \code{VarObs}, \code{BoolCrit}, \code{transfo}, \code{Weights}. If the list format is chosen, all the lists must have the same length. \cr
Calculation of a single criterion (e.g. NSE computed on discharge) is prepared by providing to \code{CreateInputsCrit} arguments atomics only. \cr
Calculation of multiple criteria (e.g. NSE computed on discharge and RMSE computed on discharge) is prepared by providing to \code{CreateInputsCrit} arguments lists except for \code{Weights} that must be set as \code{NULL}. \cr
diff --git a/man/figures/diagramGR4J-EN.pdf b/man/figures/diagramGR4J-EN.pdf
index 47e7cb57b50f0143e306aa5649525d6a44f2158e..f442f181dcbef255f8f4740282f08cccc518b34a 100644
Binary files a/man/figures/diagramGR4J-EN.pdf and b/man/figures/diagramGR4J-EN.pdf differ
diff --git a/man/figures/diagramGR4J-EN.png b/man/figures/diagramGR4J-EN.png
index 7aaf6abe4071ad9c28af52dc10f9ef63b18da849..61ea032594d5b70b08466be9f08186becc8dbdb8 100644
Binary files a/man/figures/diagramGR4J-EN.png and b/man/figures/diagramGR4J-EN.png differ
diff --git a/man/figures/diagramGR5J-EN.pdf b/man/figures/diagramGR5J-EN.pdf
index 0db91d2a60b93c786d5705a9152f6c9e450a6563..7b300d046ed22e9d568cfe7c0e2bc8a04576dabe 100644
Binary files a/man/figures/diagramGR5J-EN.pdf and b/man/figures/diagramGR5J-EN.pdf differ
diff --git a/man/figures/diagramGR5J-EN.png b/man/figures/diagramGR5J-EN.png
index 184a01ee07374144cd1a61b471a58df074da8e24..3b9b0c30538424b0f445ff0d5dfaaf62e9ed25ae 100644
Binary files a/man/figures/diagramGR5J-EN.png and b/man/figures/diagramGR5J-EN.png differ
diff --git a/man/figures/diagramGR6J-EN.pdf b/man/figures/diagramGR6J-EN.pdf
index c6d08b57e6cdfd415e45fe75fefab498d6a07691..f122ff4044bd5862352320b1e6a2644e1adbe4c9 100644
Binary files a/man/figures/diagramGR6J-EN.pdf and b/man/figures/diagramGR6J-EN.pdf differ
diff --git a/man/figures/diagramGR6J-EN.png b/man/figures/diagramGR6J-EN.png
index 90d0e3651653e84f4b26887e866dff680005e918..f9a115fad1e78b40a0e22c5c0c30c4f4a802ba26 100644
Binary files a/man/figures/diagramGR6J-EN.png and b/man/figures/diagramGR6J-EN.png differ
diff --git a/src/frun_CEMANEIGE.f b/src/frun_CEMANEIGE.f90
similarity index 55%
rename from src/frun_CEMANEIGE.f
rename to src/frun_CEMANEIGE.f90
index c4ca0aeaefb272a33635313f20d783d046c19ff5..b19e1aae4c9598601a25e012454d500bacd9b58b 100644
--- a/src/frun_CEMANEIGE.f
+++ b/src/frun_CEMANEIGE.f90
@@ -1,52 +1,91 @@
+!------------------------------------------------------------------------------
+! Subroutines relative to the annual GR4J model
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : frun_CEMANEIGE.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: A. Valéry, P. Riboust
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2011
+! Last modified: 22/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! Riboust, P., G. Thirel, N. Le Moine and P. Ribstein (2019), Revisiting a
+! simple degree-day model for integrating satellite data: implementation of
+! SWE-SCA hystereses. Journal of Hydrology and Hydromechanics,
+! doi:10.2478/johh-2018-0004, 67, 1, 70–81.
+!
+! Valéry, A., V. Andréassian and C. Perrin (2014), "As simple as possible but
+! not simpler": what is useful in a temperature-based snow-accounting routine?
+! Part 1 - Comparison of six snow accounting routines on 380 catchments,
+! Journal of Hydrology, doi:10.1016/j.jhydrol.2014.04.059.
+!
+! Valéry, A., V. Andréassian and C. Perrin (2014), "As simple as possible but
+! not simpler": What is useful in a temperature-based snow-accounting routine?
+! Part 2 - Sensitivity analysis of the Cemaneige snow accounting routine on
+! 380 catchments, Journal of Hydrology, doi:10.1016/j.jhydrol.2014.04.058.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+! 1. frun_cemaneige
+!------------------------------------------------------------------------------
+
+
+ SUBROUTINE frun_cemaneige(LInputs,InputsPrecip, &
+ InputsFracSolidPrecip,InputsTemp, &
+ MeanAnSolidPrecip,NParam,Param,NStates, &
+ StateStart,IsHyst,NOutputs,IndOutputs, &
+ Outputs,StateEnd)
+! Subroutine that runs the CemaNeige model 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/time step]
+! InputsFracSolidPrecip ! Vector of real, input series of fraction of solid precipitation [0-1]
+! InputsTemp ! Vector of real, input series of air mean temperature [degC]
+! MeanAnSolidPrecip ! Real, value of annual mean solid precip [mm/y]
+! 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 [-] and thresholds [mm])
+! IsHyst ! integer, whether we should use the linear hysteresis [1] or not [0]
+! 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 [-] and thresholds [mm])
-
- SUBROUTINE frun_cemaneige(
- ! inputs
- & LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm/time step]
- & InputsFracSolidPrecip, ! [double] input series of fraction of solid precipitation [0-1]
- & InputsTemp , ! [double] input series of air mean temperature [degC]
- & MeanAnSolidPrecip , ! [double] value of annual mean solid precip [mm/y]
- & NParam , ! [integer] number of model parameter
- & Param , ! [double] parameter set
- & NStates , ! [integer] number of state variables used for model initialization = 4
- & StateStart , ! [double] state variables used when the model run starts
- & IsHyst , ! [integer] whether we should use the linear hysteresis or not
- & NOutputs , ! [integer] number of output series
- & IndOutputs , ! [integer] indices of output series
- ! outputs
- & Outputs , ! [double] output series
- & StateEnd ) ! [double] state variables at the end of the model run
-
-
!DEC$ ATTRIBUTES DLLEXPORT :: frun_cemaneige
Implicit None
- ! input and output variables
+
+ !! dummies
+ ! in
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, intent(in) :: MeanAnSolidPrecip
doubleprecision, intent(in), dimension(LInputs) :: InputsPrecip
- doubleprecision, intent(in), dimension(LInputs) ::
- & InputsFracSolidPrecip
+ doubleprecision, intent(in), dimension(LInputs) :: InputsFracSolidPrecip
doubleprecision, intent(in), dimension(LInputs) :: InputsTemp
doubleprecision, intent(in), dimension(NParam) :: Param
doubleprecision, intent(in), dimension(NStates) :: StateStart
- integer, intent(in) :: IsHyst ! 1 if linear hysteresis is used, 0 otherwise
- logical :: IsHystBool ! TRUE if linear hysteresis is used, FALSE otherwise
- doubleprecision, intent(out), dimension(NStates) :: StateEnd
+ integer, intent(in) :: IsHyst ! 1 if linear hysteresis is used, 0 otherwise
integer, intent(in), dimension(NOutputs) :: IndOutputs
- doubleprecision, intent(out), dimension(LInputs,NOutputs) ::
- & Outputs
+ ! out
+ doubleprecision, intent(out), dimension(NStates) :: StateEnd
+ doubleprecision, intent(out), dimension(LInputs,NOutputs) :: Outputs
- ! parameters, internal states and variables
- doubleprecision CTG,Kf
- doubleprecision G,eTG,PliqAndMelt,dG,prct
- doubleprecision Tmelt,Gthreshold,MinSpeed,Gacc,Glocalmax
- doubleprecision Pliq,Psol,Gratio,PotMelt,Melt,Ginit
- integer I,K
+ !! locals
+ logical :: IsHystBool ! TRUE if linear hysteresis is used, FALSE otherwise
+ doubleprecision :: CTG,Kf
+ doubleprecision :: G,eTG,PliqAndMelt,dG,prct
+ doubleprecision :: Tmelt,Gthreshold,MinSpeed,Gacc,Glocalmax
+ doubleprecision :: Pliq,Psol,Gratio,PotMelt,Melt,Ginit
+ integer :: I,K
IF (IsHyst .EQ. 1) IsHystBool = .TRUE.
IF (IsHyst .EQ. 0) IsHystBool = .FALSE.
@@ -180,5 +219,5 @@
RETURN
- ENDSUBROUTINE
+ END SUBROUTINE
diff --git a/src/frun_GR1A.f b/src/frun_GR1A.f
deleted file mode 100644
index 00df0efe7400e8132b73eeab261a7c814a418475..0000000000000000000000000000000000000000
--- a/src/frun_GR1A.f
+++ /dev/null
@@ -1,123 +0,0 @@
-
-
- SUBROUTINE frun_gr1a(
- ! inputs
- & LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm/year]
- & InputsPE , ! [double] input series of potential evapotranspiration (PE) [mm/year]
- & NParam , ! [integer] number of model parameters
- & Param , ! [double] parameter set
- & NStates , ! [integer] number of state variables
- & StateStart , ! [double] state variables used when the model run starts (none here)
- & NOutputs , ! [integer] number of output series
- & IndOutputs , ! [integer] indices of output series
- ! outputs
- & Outputs , ! [double] output series
- & StateEnd ) ! [double] state variables at the end of the model run (none here)
-
-
- !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr1a
-
-
- Implicit None
- ! input and output variables
- integer, intent(in) :: LInputs,NParam,NStates,NOutputs
- doubleprecision, dimension(LInputs) :: InputsPrecip
- doubleprecision, dimension(LInputs) :: InputsPE
- doubleprecision, dimension(NParam) :: Param
- doubleprecision, dimension(NStates) :: StateStart
- doubleprecision, dimension(NStates) :: StateEnd
- integer, dimension(NOutputs) :: IndOutputs
- doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
- ! parameters, internal states and variables
- integer NMISC
- parameter (NMISC=3)
- doubleprecision MISC(NMISC)
- doubleprecision P0,P1,E1,Q
- integer I,K
-
- !--------------------------------------------------------------
- ! Initializations
- !--------------------------------------------------------------
-
- ! parameter values
- ! Param(1) : PE adjustment factor [-]
-
- ! initialization of model outputs
- Q = -999.999
- MISC = -999.999
-! Outputs = -999.999 ! initialization made in R
-
- StateStart = -999.999 ! CRAN-compatibility updates
- StateEnd = -999.999 ! CRAN-compatibility updates
-
-
- !--------------------------------------------------------------
- ! Time loop
- !--------------------------------------------------------------
- DO k=2,LInputs
- P0=InputsPrecip(k-1)
- P1=InputsPrecip(k)
- E1=InputsPE(k)
-! Q = -999.999
-! MISC = -999.999
- ! model run on one time step
- CALL MOD_GR1A(Param,P0,P1,E1,Q,MISC)
- ! storage of outputs
- DO I=1,NOutputs
- Outputs(k,I)=MISC(IndOutputs(I))
- ENDDO
- ENDDO
-
- RETURN
-
- ENDSUBROUTINE
-
-
-
-
-
-!################################################################################################################################
-
-
-
-
-!**********************************************************************
- SUBROUTINE MOD_GR1A(Param,P0,P1,E1,Q,MISC)
-! Calculation of streamflow on a single time step with the GR1A model
-! Inputs:
-! Param Vector of model parameters (Param(1) [-])
-! P0 Value of rainfall during the previous time step [mm/year]
-! P1 Value of rainfall during the current time step [mm/year]
-! E1 Value of potential evapotranspiration during the current time step [mm/year]
-! Outputs:
-! Q Value of simulated flow at the catchment outlet for the current time step [mm/year]
-! MISC Vector of model outputs for the time step [mm/year]
-!**********************************************************************
- Implicit None
- INTEGER NMISC,NParam
- PARAMETER (NMISC=3)
- PARAMETER (NParam=1)
- DOUBLEPRECISION Param(NParam)
- DOUBLEPRECISION MISC(NMISC)
- DOUBLEPRECISION P0,P1,E1,Q
-
- DOUBLEPRECISION tt ! speed-up
-
-! Runoff
- ! speed-up
- tt = (0.7*P1+0.3*P0)/Param(1)/E1
- Q=P1*(1.-1./SQRT(1.+tt*tt))
- ! Q=P1*(1.-1./(1.+((0.7*P1+0.3*P0)/Param(1)/E1)**2.)**0.5)
- ! end speed-up
-
-! Variables storage
- MISC( 1)=E1 ! PE ! [numeric] observed potential evapotranspiration [mm/year]
- MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/year]
- MISC( 3)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/year]
-
-
- ENDSUBROUTINE
-
-
diff --git a/src/frun_GR1A.f90 b/src/frun_GR1A.f90
new file mode 100644
index 0000000000000000000000000000000000000000..a6ff7cf89d0b13b7bbb77597a225351394ac68fc
--- /dev/null
+++ b/src/frun_GR1A.f90
@@ -0,0 +1,159 @@
+!------------------------------------------------------------------------------
+! Subroutines relative to the annual GR1A model
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : frun_GR1A.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: S. Mouelhi
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2003
+! Last modified: 25/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! Mouelhi S. (2003). Vers une chaîne cohérente de modèles pluie-débit
+! conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et
+! journalier. PhD thesis (in French), ENGREF, Cemagref Antony, France.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+! 1. frun_gr1a
+! 2. MOD_GR1A
+!------------------------------------------------------------------------------
+
+
+ SUBROUTINE frun_gr1a(LInputs,InputsPrecip,InputsPE,NParam,Param, &
+ NStates,StateStart,NOutputs,IndOutputs, &
+ Outputs,StateEnd)
+! Subroutine that initializes GR1A, get its parameters, performs the call
+! to the MOD_GR1A 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/year]
+! InputsPE ! Vector of real, input series of potential evapotranspiration (PE) [mm/year]
+! 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 (none here)
+! 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 (none here)
+
+
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr1a
+
+
+ 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
+ integer, dimension(NOutputs), intent(in) :: IndOutputs
+ ! out
+ doubleprecision, dimension(NStates), intent(out) :: StateEnd
+ doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
+
+ !! locals
+ integer :: I,K
+ integer, parameter :: NMISC=3
+ doubleprecision, dimension(NMISC) :: MISC
+ doubleprecision P0,P1,E1,Q
+
+ !--------------------------------------------------------------
+ ! Initializations
+ !--------------------------------------------------------------
+
+ ! parameter values
+ ! Param(1) : PE adjustment factor [-]
+
+ ! initialization of model outputs
+ Q = -999.999
+ MISC = -999.999
+! Outputs = -999.999 ! initialization made in R
+
+! StateStart = -999.999 ! CRAN-compatibility updates
+! StateEnd = -999.999 ! CRAN-compatibility updates
+
+
+ !--------------------------------------------------------------
+ ! Time loop
+ !--------------------------------------------------------------
+ DO k=2,LInputs
+ P0=InputsPrecip(k-1)
+ P1=InputsPrecip(k)
+ E1=InputsPE(k)
+! Q = -999.999
+! MISC = -999.999
+ ! model run on one time step
+ CALL MOD_GR1A(Param,P0,P1,E1,Q,MISC)
+ ! storage of outputs
+ DO I=1,NOutputs
+ Outputs(k,I)=MISC(IndOutputs(I))
+ ENDDO
+ ENDDO
+
+ RETURN
+
+ ENDSUBROUTINE
+
+
+
+
+
+!################################################################################################################################
+
+
+
+
+!**********************************************************************
+ SUBROUTINE MOD_GR1A(Param,P0,P1,E1,Q,MISC)
+! Calculation of streamflow on a single time step (year) with the GR1A model
+! Inputs
+! Param ! Vector of real, model parameters (Param(1) [-])
+! P0 ! Real, value of rainfall during the previous time step [mm/year]
+! P1 ! Real, value of rainfall during the current time step [mm/year]
+! E1 ! Real, value of potential evapotranspiration during the current time step [mm/year]
+! Outputs
+! Q ! Real, value of simulated flow at the catchment outlet for the current time step [mm/year]
+! MISC ! Vector of real, model outputs for the time step [mm/year]
+!**********************************************************************
+ Implicit None
+
+ !! locals
+ integer, parameter :: NMISC=3
+ integer, parameter :: NParam=1
+ doubleprecision :: tt ! speed-up
+
+ !! dummies
+ ! in
+ doubleprecision, dimension(NParam), intent(in) :: Param
+ doubleprecision, intent(in) :: P0,P1,E1
+ ! out
+ doubleprecision, dimension(NMISC), intent(out) :: MISC
+ doubleprecision, intent(out) :: Q
+
+
+! Runoff
+ ! speed-up
+ tt = (0.7*P1+0.3*P0)/Param(1)/E1
+ Q=P1*(1.-1./SQRT(1.+tt*tt))
+ ! Q=P1*(1.-1./(1.+((0.7*P1+0.3*P0)/Param(1)/E1)**2.)**0.5)
+ ! end speed-up
+
+! Variables storage
+ MISC( 1)=E1 ! PE ! [numeric] observed potential evapotranspiration [mm/year]
+ MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/year]
+ MISC( 3)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/year]
+
+
+ END SUBROUTINE
+
+
diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f
deleted file mode 100644
index 904d9befe5003a71a3fa6535d2b4ee83470390c3..0000000000000000000000000000000000000000
--- a/src/frun_GR2M.f
+++ /dev/null
@@ -1,183 +0,0 @@
-
-
- SUBROUTINE frun_gr2m(
- ! inputs
- & LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm/month]
- & InputsPE , ! [double] input series of potential evapotranspiration (PE) [mm/month]
- & NParam , ! [integer] number of model parameters
- & Param , ! [double] parameter set
- & NStates , ! [integer] number of state variables
- & StateStart , ! [double] state variables used when the model run starts (store levels [mm])
- & NOutputs , ! [integer] number of output series
- & IndOutputs , ! [integer] indices of output series
- ! outputs
- & Outputs , ! [double] output series
- & StateEnd ) ! [double] state variables at the end of the model run (store levels [mm])
-
-
- !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr2m
-
-
- Implicit None
- ! input and output variables
- integer, intent(in) :: LInputs,NParam,NStates,NOutputs
- doubleprecision, dimension(LInputs) :: InputsPrecip
- doubleprecision, dimension(LInputs) :: InputsPE
- doubleprecision, dimension(NParam) :: Param
- doubleprecision, dimension(NStates) :: StateStart
- doubleprecision, dimension(NStates) :: StateEnd
- integer, dimension(NOutputs) :: IndOutputs
- doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
- ! parameters, internal states and variables
- integer NMISC
- parameter (NMISC=30)
- doubleprecision St(2)
- doubleprecision MISC(NMISC)
- doubleprecision P,E,Q
- integer I,K
-
- !--------------------------------------------------------------
- ! Initializations
- !--------------------------------------------------------------
-
- ! initialization of model states to zero
- St=0.
-
- ! initialization of model states using StateStart
- St(1)=StateStart(1)
- St(2)=StateStart(2)
-
- ! parameter values
- ! Param(1) : production store capacity [mm]
- ! Param(2) : groundwater exchange coefficient [-]
-
- ! 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
- P =InputsPrecip(k)
- E =InputsPE(k)
-! Q = -999.999
-! MISC = -999.999
- !model run on one time step
- CALL MOD_GR2M(St,Param,P,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)
-
- RETURN
-
- ENDSUBROUTINE
-
-
-
-
-
-!################################################################################################################################
-
-
-
-
-!**********************************************************************
- SUBROUTINE MOD_GR2M(St,Param,P,E,Q,MISC)
-! Calculation of streamflow on a single time step (month) with the GR2M model
-! Inputs:
-! St Vector of model states at the beginning of the time step [mm]
-! Param Vector of model parameters (Param(1) [mm]; Param(2) [-])
-! P Value of rainfall during the time step [mm/month]
-! E Value of potential evapotranspiration during the time step [mm/month]
-! Outputs:
-! St Vector of model states at the end of the time step [mm]
-! Q Value of simulated flow at the catchment outlet for the time step [mm/month]
-! MISC Vector of model outputs for the time step [mm/month]
-!**********************************************************************
- Implicit None
- INTEGER NMISC,NParam
- PARAMETER (NMISC=30)
- PARAMETER (NParam=2)
- DOUBLEPRECISION St(2)
- DOUBLEPRECISION Param(NParam)
- DOUBLEPRECISION MISC(NMISC)
- DOUBLEPRECISION P,E,Q
- DOUBLEPRECISION WS,tanHyp,S1,S2
- DOUBLEPRECISION P1,P2,P3,R1,R2,AE,EXCH
-
- DOUBLEPRECISION TWS, Sr ! speed-up
-
-! Production store
- WS=P/Param(1)
- IF(WS.GT.13.)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- S1=(St(1)+Param(1)*TWS)/(1.+St(1)/Param(1)*TWS)
- ! S1=(X(1)+Param(1)*tanHyp(WS))/(1.+X(1)/Param(1)*tanHyp(WS))
- ! end speed-up
-
- P1=P+St(1)-S1
- WS=E/Param(1)
- IF(WS.GT.13.)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- S2=S1*(1.-TWS)/(1.+(1.-S1/Param(1))*TWS)
- ! S2=S1*(1.-tanHyp(WS))/(1.+(1.-S1/Param(1))*tanHyp(WS))
- ! end speed-up
- AE = S1 - S2
-
-! Percolation
- ! speed-up
- Sr = S2/Param(1)
- Sr = Sr * Sr * Sr + 1.
- St(1)=S2/Sr**(1./3.)
- ! X(1)=S2/(1+(S2/Param(1))**3.)**(1./3.)
- ! end speed-up
-
- P2=S2-St(1)
- P3=P1+P2
-
-! QR calculation (routing store)
- R1=St(2)+P3
-
-! Water exchange
- R2=Param(2)*R1
- EXCH = R2 - R1
-
-! Total runoff
- Q=R2*R2/(R2+60.)
-
-! Updating store level
- St(2)=R2-Q
-
-
-! Variables storage
- MISC( 1)=E ! PE ! [numeric] observed potential evapotranspiration [mm/month]
- MISC( 2)=P ! Precip ! [numeric] observed total precipitation [mm/month]
- MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
- MISC( 4)=P1 ! Pn ! [numeric] net rainfall (P1) [mm/month]
- MISC( 5)=AE ! AE ! [numeric] actual evapotranspiration [mm/month]
- MISC( 6)=P2 ! Perc ! [numeric] percolation (P2) [mm/month]
- MISC( 7)=P3 ! PR ! [numeric] P3=P1+P2 [mm/month]
- MISC( 8)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
- MISC( 9)=EXCH ! EXCH ! [numeric] groundwater exchange (EXCH) [mm/month]
- MISC(10)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/month]
-
-
- ENDSUBROUTINE
-
-
diff --git a/src/frun_GR2M.f90 b/src/frun_GR2M.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9be6de91b8e8522b920a7d9a7f1158ac06dc0baa
--- /dev/null
+++ b/src/frun_GR2M.f90
@@ -0,0 +1,224 @@
+!------------------------------------------------------------------------------
+! Subroutines relative to the annual GR2M model
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : frun_GR2M.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: S. Mouelhi
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2003
+! Last modified: 25/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! Mouelhi S. (2003). Vers une chaîne cohérente de modèles pluie-débit
+! conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et
+! journalier. PhD thesis (in French), ENGREF, Cemagref Antony, France.
+!
+! Mouelhi, S., C. Michel, C. Perrin and V. Andréassian (2006). Stepwise
+! development of a two-parameter monthly water balance model. Journal of
+! Hydrology, 318(1-4), 200-214. doi:10.1016/j.jhydrol.2005.06.014.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+! 1. frun_gr2m
+! 2. MOD_GR2M
+!------------------------------------------------------------------------------
+
+
+ SUBROUTINE frun_gr2m(LInputs,InputsPrecip,InputsPE,NParam,Param, &
+ NStates,StateStart,NOutputs,IndOutputs, &
+ Outputs,StateEnd)
+! Subroutine that initializes GR2M, get its parameters, performs the call
+! to the MOD_GR2M 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/month]
+! InputsPE ! Vector of real, input series of potential evapotranspiration (PE) [mm/month]
+! 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])
+! 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])
+
+
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr2m
+
+
+ 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
+ integer, dimension(NOutputs), intent(in) :: IndOutputs
+ ! out
+ doubleprecision, dimension(NStates), intent(out) :: StateEnd
+ doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
+
+ !! locals
+ integer :: I,K
+ integer, parameter :: NMISC=30
+ doubleprecision, dimension(2) :: St
+ doubleprecision, dimension(NMISC) :: MISC
+ doubleprecision :: P,E,Q
+
+ !--------------------------------------------------------------
+ ! Initializations
+ !--------------------------------------------------------------
+
+ ! initialization of model states to zero
+ St=0.
+
+ ! initialization of model states using StateStart
+ St(1)=StateStart(1)
+ St(2)=StateStart(2)
+
+ ! parameter values
+ ! Param(1) : production store capacity [mm]
+ ! Param(2) : groundwater exchange coefficient [-]
+
+ ! 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
+ P =InputsPrecip(k)
+ E =InputsPE(k)
+! Q = -999.999
+! MISC = -999.999
+ !model run on one time step
+ CALL MOD_GR2M(St,Param,P,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)
+
+ RETURN
+
+ ENDSUBROUTINE
+
+
+
+
+
+!################################################################################################################################
+
+
+
+
+!**********************************************************************
+ SUBROUTINE MOD_GR2M(St,Param,P,E,Q,MISC)
+! Calculation of streamflow on a single time step (month) with the GR2M model
+! Inputs:
+! St Vector of real, model states at the beginning of the time step [mm]
+! Param Vector of real, model parameters (Param(1) [mm]; Param(2) [-])
+! P Real, value of rainfall during the time step [mm/month]
+! E Real, value of potential evapotranspiration during the time step [mm/month]
+! Outputs:
+! St Vector of real, model states at the end of the time step [mm]
+! Q Real, value of simulated flow at the catchment outlet for the time step [mm/month]
+! MISC Vector of real, model outputs for the time step [mm/month]
+!**********************************************************************
+ Implicit None
+
+ !! locals
+ integer, parameter :: NParam=2,NMISC=30
+ doubleprecision :: WS,tanHyp,S1,S2
+ doubleprecision :: P1,P2,P3,R1,R2,AE,EXCH
+ doubleprecision :: expWS, TWS, Sr ! speed-up
+
+ !! dummies
+ ! in
+ doubleprecision, dimension(NParam), intent(in) :: Param
+ doubleprecision, intent(in) :: P,E
+ ! inout
+ doubleprecision, dimension(2), intent(inout) :: St
+ ! out
+ doubleprecision, intent(out) :: Q
+ doubleprecision, dimension(NMISC), intent(out) :: MISC
+
+! Production store
+ WS=P/Param(1)
+ IF(WS.GT.13.) WS=13.
+
+ ! speed-up
+ expWS = exp(2.*WS)
+ TWS = (expWS - 1.)/(expWS + 1.)
+ S1=(St(1)+Param(1)*TWS)/(1.+St(1)/Param(1)*TWS)
+ ! S1=(X(1)+Param(1)*tanHyp(WS))/(1.+X(1)/Param(1)*tanHyp(WS))
+ ! end speed-up
+
+ P1=P+St(1)-S1
+ WS=E/Param(1)
+ IF(WS.GT.13.) WS=13.
+
+ ! speed-up
+ expWS = exp(2.*WS)
+ TWS = (expWS - 1.)/(expWS + 1.)
+ S2=S1*(1.-TWS)/(1.+(1.-S1/Param(1))*TWS)
+ ! S2=S1*(1.-tanHyp(WS))/(1.+(1.-S1/Param(1))*tanHyp(WS))
+ ! end speed-up
+ AE = S1 - S2
+
+! Percolation
+ ! speed-up
+ Sr = S2/Param(1)
+ Sr = Sr * Sr * Sr + 1.
+ St(1)=S2/Sr**(1./3.)
+ ! X(1)=S2/(1+(S2/Param(1))**3.)**(1./3.)
+ ! end speed-up
+
+ P2=S2-St(1)
+ P3=P1+P2
+
+! QR calculation (routing store)
+ R1=St(2)+P3
+
+! Water exchange
+ R2=Param(2)*R1
+ EXCH = R2 - R1
+
+! Total runoff
+ Q=R2*R2/(R2+60.)
+
+! Updating store level
+ St(2)=R2-Q
+
+
+! Variables storage
+ MISC( 1)=E ! PE ! [numeric] observed potential evapotranspiration [mm/month]
+ MISC( 2)=P ! Precip ! [numeric] observed total precipitation [mm/month]
+ MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
+ MISC( 4)=P1 ! Pn ! [numeric] net rainfall (P1) [mm/month]
+ MISC( 5)=AE ! AE ! [numeric] actual evapotranspiration [mm/month]
+ MISC( 6)=P2 ! Perc ! [numeric] percolation (P2) [mm/month]
+ MISC( 7)=P3 ! PR ! [numeric] P3=P1+P2 [mm/month]
+ MISC( 8)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
+ MISC( 9)=EXCH ! EXCH ! [numeric] groundwater exchange (EXCH) [mm/month]
+ MISC(10)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/month]
+
+
+ END SUBROUTINE
+
+
diff --git a/src/frun_GR4H.f b/src/frun_GR4H.f
deleted file mode 100644
index 5ffe19e4cf3e35a267c70b8b73ec38bb2d77fd59..0000000000000000000000000000000000000000
--- a/src/frun_GR4H.f
+++ /dev/null
@@ -1,283 +0,0 @@
-
-
- SUBROUTINE frun_gr4h(
- ! inputs
- & LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm/hour]
- & InputsPE , ! [double] input series of potential evapotranspiration (PE) [mm/hour]
- & NParam , ! [integer] number of model parameters
- & Param , ! [double] parameter set
- & NStates , ! [integer] number of state variables
- & StateStart , ! [double] state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
- & NOutputs , ! [integer] number of output series
- & IndOutputs , ! [integer] indices of output series
- ! outputs
- & Outputs , ! [double] output series
- & StateEnd ) ! [double] state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
-
-
- !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr4h
-
-
- Implicit None
- ! input and output variables
- integer, intent(in) :: LInputs,NParam,NStates,NOutputs
- doubleprecision, dimension(LInputs) :: InputsPrecip
- doubleprecision, dimension(LInputs) :: InputsPE
- doubleprecision, dimension(NParam) :: Param
- doubleprecision, dimension(NStates) :: StateStart
- doubleprecision, dimension(NStates) :: StateEnd
- integer, dimension(NOutputs) :: IndOutputs
- doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
- ! parameters, internal states and variables
- integer NH,NMISC
- parameter (NH=480,NMISC=30)
- doubleprecision St(2), StUH1(NH), StUH2(2*NH)
- doubleprecision OrdUH1(NH), OrdUH2(2*NH)
- doubleprecision MISC(NMISC)
- doubleprecision D
- doubleprecision P1,E,Q
- integer I,K
-
- !--------------------------------------------------------------
- ! Initializations
- !--------------------------------------------------------------
-
- ! initialization of model states using StateStart
- St=0.
- StUH1=0.
- StUH2=0.
-
- ! initialization of model states using StateStart
- St(1) = StateStart(1)
- St(2) = StateStart(2)
- DO I=1,NH
- StUH1(I)=StateStart(7+I)
- ENDDO
- DO I=1,2*NH
- StUH2(I)=StateStart(7+I+NH)
- ENDDO
-
- ! parameter values
- ! Param(1) : production store capacity (X1 - PROD) [mm]
- ! Param(2) : intercatchment exchange coefficient (X2 - CES) [mm/hour]
- ! Param(3) : routing store capacity (X3 - ROUT) [mm]
- ! Param(4) : time constant of unit hydrograph (X4 - TB) [hour]
-
- !computation of UH ordinates
- OrdUH1 = 0.
- OrdUH2 = 0.
-
- D=1.25
- CALL UH1_H(OrdUH1,Param(4),D)
- 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_GR4H(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,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)
- DO K=1,NH
- StateEnd(7+K)=StUH1(K)
- ENDDO
- DO K=1,2*NH
- StateEnd(7+NH+K)=StUH2(K)
- ENDDO
-
- RETURN
-
- ENDSUBROUTINE
-
-
-
-
-
-!################################################################################################################################
-
-
-
-
-!**********************************************************************
- SUBROUTINE MOD_GR4H(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
- &MISC)
-! Run on a single time step with the GR4H model
-! Inputs:
-! St Vector of model states in stores at the beginning of the time step [mm]
-! StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-! OrdUH1 Vector of ordinates in UH1 [-]
-! OrdUH2 Vector of ordinates in UH2 [-]
-! Param Vector of model parameters [various units]
-! P1 Value of rainfall during the time step [mm]
-! E Value of potential evapotranspiration during the time step [mm]
-! Outputs:
-! St Vector of model states in stores at the beginning of the time step [mm]
-! StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-! Q Value of simulated flow at the catchment outlet for the time step [mm/hour]
-! MISC Vector of model outputs for the time step [mm/hour]
-!**********************************************************************
- Implicit None
- INTEGER NH,NMISC,NParam
- PARAMETER (NH=480,NMISC=30)
- PARAMETER (NParam=4)
- DOUBLEPRECISION St(2),StUH1(NH)
- DOUBLEPRECISION StUH2(2*NH),OrdUH1(NH),OrdUH2(2*NH)
- DOUBLEPRECISION Param(NParam)
- DOUBLEPRECISION MISC(NMISC)
- DOUBLEPRECISION P1,E,Q
- DOUBLEPRECISION A,B,EN,ER,PN,PR,PS,WS,tanHyp
- DOUBLEPRECISION PERC,PRHU1,PRHU2,EXCH,QR,QD
- DOUBLEPRECISION AE,AEXCH1,AEXCH2
- INTEGER K
-
- DATA B/0.9/
-
- DOUBLEPRECISION TWS, Sr, Rr ! speed-up
-
- A=Param(1)
-
-
-! Interception and production store
- IF(P1.LE.E) THEN
- EN=E-P1
- PN=0.
- WS=EN/A
- IF(WS.GT.13.)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- 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
- St(1)=St(1)-ER
- PR=0.
- ELSE
- EN=0.
- AE=E
- PN=P1-E
- WS=PN/A
- IF(WS.GT.13.)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- 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
-
-! Percolation from production store
- IF(St(1).LT.0.)St(1)=0.
-
- ! speed-up
- ! (21/4)**4 = 759.69140625
- Sr = St(1)/Param(1)
- Sr = Sr * Sr
- Sr = Sr * Sr
- PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/759.69140625)))
- ! 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
-
-! Split of effective rainfall into the two routing components
- PRHU1=PR*B
- PRHU2=PR*(1.-B)
-
-! Convolution of unit hydrograph UH1
- DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
- StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
- ENDDO
- StUH1(NH)=OrdUH1(NH)*PRHU1
-
-! 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)*PRHU2
- ENDDO
- StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
-
-! Potential intercatchment semi-exchange
- ! speed-up
- Rr = St(2)/Param(3)
- EXCH=Param(2)*Rr*Rr*Rr*SQRT(Rr)
- ! EXCH=Param(2)*(X(1)/Param(3))**3.5
- ! end speed-up
-
-! Routing store
- AEXCH1=EXCH
- IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
- St(2)=St(2)+StUH1(1)+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((StUH2(1)+EXCH).LT.0.) AEXCH2=-StUH2(1)
- QD=MAX(0.d0,StUH2(1)+EXCH)
-
-! Total runoff
- Q=QR+QD
- IF(Q.LT.0.) Q=0.
-
-! Variables storage
- MISC( 1)=E ! PE ! [numeric] observed potential evapotranspiration [mm/hour]
- MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/hour]
- MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
- MISC( 4)=AE ! AE ! [numeric] actual evapotranspiration [mm/hour]
- MISC( 5)=PERC ! Perc ! [numeric] percolation (PERC) [mm/hour]
- MISC( 6)=PR ! PR ! [numeric] PR=PN-PS+PERC [mm/hour]
- MISC( 7)=StUH1(1) ! Q9 ! [numeric] outflow from UH1 (Q9) [mm/hour]
- MISC( 8)=StUH2(1) ! Q1 ! [numeric] outflow from UH2 (Q1) [mm/hour]
- MISC( 9)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
- MISC(10)=EXCH ! Exch ! [numeric] potential semi-exchange between catchments (EXCH) [mm/hour]
- MISC(11)=AEXCH1+AEXCH2 ! AExch ! [numeric] actual total exchange between catchments (AEXCH1+AEXCH2) [mm/hour]
- MISC(12)=QR ! QR ! [numeric] outflow from routing store (QR) [mm/hour]
- MISC(13)=QD ! QD ! [numeric] outflow from UH2 branch after exchange (QD) [mm/hour]
- MISC(14)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/hour]
-
-
-
-
- ENDSUBROUTINE
-
-
diff --git a/src/frun_GR4H.f90 b/src/frun_GR4H.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ec12cc787151ab74822469c44407ca9e945a16cd
--- /dev/null
+++ b/src/frun_GR4H.f90
@@ -0,0 +1,322 @@
+!------------------------------------------------------------------------------
+! Subroutines relative to the annual GR4H model
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : frun_GR4H.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: C. Perrin
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2003
+! Last modified: 25/11/2019
+!------------------------------------------------------------------------------
+! 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.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+! 1. frun_gr4h
+! 2. MOD_GR4H
+!------------------------------------------------------------------------------
+
+
+ SUBROUTINE frun_gr4h(LInputs,InputsPrecip,InputsPE,NParam,Param, &
+ NStates,StateStart,NOutputs,IndOutputs, &
+ Outputs,StateEnd)
+! Subroutine that initializes GR4H, get its parameters, performs the call
+! to the MOD_GR4H 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])
+! 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_gr4h
+
+
+ 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
+ integer, dimension(NOutputs), intent(in) :: IndOutputs
+ ! out
+ doubleprecision, dimension(NStates), intent(out) :: StateEnd
+ doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
+
+ !! locals
+ integer :: I,K
+ integer, parameter :: NH=480,NMISC=30
+ doubleprecision, dimension(2) :: St
+ doubleprecision, dimension(NH) :: StUH1, OrdUH1
+ doubleprecision, dimension(2*NH) :: StUH2, OrdUH2
+ doubleprecision, dimension(NMISC) :: MISC
+ doubleprecision :: D,P1,E,Q
+
+ !--------------------------------------------------------------
+ ! Initializations
+ !--------------------------------------------------------------
+
+ ! initialization of model states using StateStart
+ St=0.
+ StUH1=0.
+ StUH2=0.
+
+ ! initialization of model states using StateStart
+ St(1) = StateStart(1)
+ St(2) = StateStart(2)
+ DO I=1,NH
+ StUH1(I)=StateStart(7+I)
+ ENDDO
+ DO I=1,2*NH
+ StUH2(I)=StateStart(7+I+NH)
+ ENDDO
+
+ ! parameter values
+ ! Param(1) : production store capacity (X1 - PROD) [mm]
+ ! Param(2) : intercatchment exchange coefficient (X2 - CES) [mm/hour]
+ ! Param(3) : routing store capacity (X3 - ROUT) [mm]
+ ! Param(4) : time constant of unit hydrograph (X4 - TB) [hour]
+
+ !computation of UH ordinates
+ OrdUH1 = 0.
+ OrdUH2 = 0.
+
+ D=1.25
+ CALL UH1_H(OrdUH1,Param(4),D)
+ 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_GR4H(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,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)
+ DO K=1,NH
+ StateEnd(7+K)=StUH1(K)
+ ENDDO
+ DO K=1,2*NH
+ StateEnd(7+NH+K)=StUH2(K)
+ ENDDO
+
+ RETURN
+
+ ENDSUBROUTINE
+
+
+
+
+
+!################################################################################################################################
+
+
+
+
+!**********************************************************************
+ SUBROUTINE MOD_GR4H(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,MISC)
+! Calculation of streamflow on a single time step (hour) with the GR4H model
+! Inputs:
+! St Vector of real, model states in stores at the beginning of the time step [mm]
+! StUH1 Vector of real, model states in Unit Hydrograph 1 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]
+! OrdUH1 Vector of real, ordinates in UH1 [-]
+! OrdUH2 Vector of real, ordinates in UH2 [-]
+! Param Vector of real, model parameters [various units]
+! P1 Real, value of rainfall during the time step [mm]
+! E Real, value of potential evapotranspiration during the time step [mm]
+! Outputs:
+! St Vector of real, model states in stores at the beginning of the time step [mm]
+! StUH1 Vector of real, model states in Unit Hydrograph 1 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]
+! Q Real, value of simulated flow at the catchment outlet for the time step [mm/hour]
+! MISC Vector of real, model outputs for the time step [mm/hour]
+!**********************************************************************
+ Implicit None
+
+ !! locals
+ integer, parameter :: NParam=4,NMISC=30,NH=480
+ doubleprecision :: A,EN,ER,PN,PR,PS,WS,tanHyp
+ doubleprecision :: PERC,PRHU1,PRHU2,EXCH,QR,QD
+ 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
+ doubleprecision, intent(in) :: P1,E
+ doubleprecision, dimension(NH), intent(inout) :: OrdUH1
+ doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
+ ! inout
+ doubleprecision, dimension(2), intent(inout) :: St
+ doubleprecision, dimension(NH), intent(inout) :: StUH1
+ 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(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
+ St(1)=St(1)-ER
+ PR=0.
+ ELSE
+ EN=0.
+ AE=E
+ 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
+
+! 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
+
+! Split of effective rainfall into the two routing components
+ PRHU1=PR*B
+ PRHU2=PR*(1.-B)
+
+! Convolution of unit hydrograph UH1
+ DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
+ StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
+ ENDDO
+ StUH1(NH)=OrdUH1(NH)*PRHU1
+
+! 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)*PRHU2
+ ENDDO
+ StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
+
+! Potential intercatchment semi-exchange
+ ! speed-up
+ Rr = St(2)/Param(3)
+ EXCH=Param(2)*Rr*Rr*Rr*SQRT(Rr)
+ ! EXCH=Param(2)*(X(1)/Param(3))**3.5
+ ! end speed-up
+
+! Routing store
+ AEXCH1=EXCH
+ IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
+ St(2)=St(2)+StUH1(1)+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((StUH2(1)+EXCH).LT.0.) AEXCH2=-StUH2(1)
+ QD=MAX(0.d0,StUH2(1)+EXCH)
+
+! Total runoff
+ Q=QR+QD
+ IF(Q.LT.0.) Q=0.
+
+! Variables storage
+ MISC( 1)=E ! PE ! [numeric] observed potential evapotranspiration [mm/hour]
+ MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/hour]
+ MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
+ MISC( 4)=AE ! AE ! [numeric] actual evapotranspiration [mm/hour]
+ MISC( 5)=PERC ! Perc ! [numeric] percolation (PERC) [mm/hour]
+ MISC( 6)=PR ! PR ! [numeric] PR=PN-PS+PERC [mm/hour]
+ MISC( 7)=StUH1(1) ! Q9 ! [numeric] outflow from UH1 (Q9) [mm/hour]
+ MISC( 8)=StUH2(1) ! Q1 ! [numeric] outflow from UH2 (Q1) [mm/hour]
+ MISC( 9)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
+ MISC(10)=EXCH ! Exch ! [numeric] potential semi-exchange between catchments (EXCH) [mm/hour]
+ MISC(11)=AEXCH1+AEXCH2 ! AExch ! [numeric] actual total exchange between catchments (AEXCH1+AEXCH2) [mm/hour]
+ MISC(12)=QR ! QR ! [numeric] outflow from routing store (QR) [mm/hour]
+ MISC(13)=QD ! QD ! [numeric] outflow from UH2 branch after exchange (QD) [mm/hour]
+ MISC(14)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/hour]
+
+
+
+
+ END SUBROUTINE
+
+
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f
deleted file mode 100644
index 14fb0c1d41405e61880b8f6f38329dc015e32ff3..0000000000000000000000000000000000000000
--- a/src/frun_GR4J.f
+++ /dev/null
@@ -1,279 +0,0 @@
-
-
- SUBROUTINE frun_gr4j(
- ! inputs
- & LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm/day]
- & InputsPE , ! [double] input series of potential evapotranspiration (PE) [mm/day]
- & NParam , ! [integer] number of model parameters
- & Param , ! [double] parameter set
- & NStates , ! [integer] number of state variables
- & StateStart , ! [double] state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
- & NOutputs , ! [integer] number of output series
- & IndOutputs , ! [integer] indices of output series
- ! outputs
- & Outputs , ! [double] output series
- & StateEnd ) ! [double] state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
-
-
- !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr4j
-
-
- Implicit None
- ! input and output variables
- integer, intent(in) :: LInputs,NParam,NStates,NOutputs
- doubleprecision, dimension(LInputs) :: InputsPrecip
- doubleprecision, dimension(LInputs) :: InputsPE
- doubleprecision, dimension(NParam) :: Param
- doubleprecision, dimension(NStates) :: StateStart
- doubleprecision, dimension(NStates) :: StateEnd
- integer, dimension(NOutputs) :: IndOutputs
- doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
- ! parameters, internal states and variables
- integer NH,NMISC
- parameter (NH=20,NMISC=30)
- doubleprecision St(2), StUH1(NH), StUH2(2*NH)
- doubleprecision OrdUH1(NH), OrdUH2(2*NH)
- doubleprecision MISC(NMISC)
- doubleprecision D
- doubleprecision P1,E,Q
- integer I,K
-
- !--------------------------------------------------------------
- ! Initializations
- !--------------------------------------------------------------
-
- ! initialization of model states to zero
- St=0.
- StUH1=0.
- StUH2=0.
-
- ! initialization of model states using StateStart
- St(1) = StateStart(1)
- St(2) = StateStart(2)
- DO I=1,NH
- StUH1(I)=StateStart(7+I)
- ENDDO
- DO I=1,2*NH
- StUH2(I)=StateStart(7+I+NH)
- ENDDO
-
- ! parameter values
- ! Param(1) : production store capacity (X1 - PROD) [mm]
- ! Param(2) : intercatchment exchange coefficient (X2 - CES) [mm/day]
- ! Param(3) : routing store capacity (X3 - ROUT) [mm]
- ! Param(4) : time constant of unit hydrograph (X4 - TB) [day]
-
- ! computation of UH ordinates
- OrdUH1 = 0.
- OrdUH2 = 0.
-
- D=2.5
- CALL UH1(OrdUH1,Param(4),D)
- CALL UH2(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_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,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)
- DO K=1,NH
- StateEnd(7+K)=StUH1(K)
- ENDDO
- DO K=1,2*NH
- StateEnd(7+NH+K)=StUH2(K)
- ENDDO
-
- RETURN
-
- ENDSUBROUTINE
-
-
-
-
-
-!################################################################################################################################
-
-
-
-
-!**********************************************************************
- SUBROUTINE MOD_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
- &MISC)
-! Run on a single time step with the GR4J model
-! Inputs:
-! St Vector of model states in stores at the beginning of the time step [mm]
-! StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-! OrdUH1 Vector of ordinates in UH1 [-]
-! OrdUH2 Vector of ordinates in UH2 [-]
-! Param Vector of model parameters [various units]
-! P1 Value of rainfall during the time step [mm/day]
-! E Value of potential evapotranspiration during the time step [mm/day]
-! Outputs:
-! St Vector of model states in stores at the end of the time step [mm]
-! StUH1 Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
-! StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
-! Q Value of simulated flow at the catchment outlet for the time step [mm/day]
-! MISC Vector of model outputs for the time step [mm/day]
-!**********************************************************************
- Implicit None
- INTEGER NH,NMISC,NParam
- PARAMETER (NH=20,NMISC=30)
- PARAMETER (NParam=4)
- DOUBLEPRECISION St(2),StUH1(NH),StUH2(2*NH)
- DOUBLEPRECISION OrdUH1(NH),OrdUH2(2*NH)
- DOUBLEPRECISION Param(NParam)
- DOUBLEPRECISION MISC(NMISC)
- DOUBLEPRECISION P1,E,Q
- DOUBLEPRECISION A,B,EN,ER,PN,PR,PS,WS,tanHyp
- DOUBLEPRECISION PERC,PRHU1,PRHU2,EXCH,QR,QD
- DOUBLEPRECISION AE,AEXCH1,AEXCH2
- INTEGER K
-
- DATA B/0.9/
- DOUBLEPRECISION TWS, Sr, Rr ! speed-up
-
- A=Param(1)
-
-
-! Interception and production store
- IF(P1.LE.E) THEN
- EN=E-P1
- PN=0.
- WS=EN/A
- IF(WS.GT.13.)WS=13.
- ! speed-up
- TWS = tanHyp(WS)
- 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
- St(1)=St(1)-ER
- PS=0.
- PR=0.
- ELSE
- EN=0.
- AE=E
- PN=P1-E
- WS=PN/A
- IF(WS.GT.13.)WS=13.
- ! speed-up
- TWS = tanHyp(WS)
- 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
-
-! Percolation from production store
- IF(St(1).LT.0.)St(1)=0.
- ! speed-up
- ! (9/4)**4 = 25.62891
- Sr = St(1)/Param(1)
- Sr = Sr * Sr
- Sr = Sr * Sr
- PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62891)))
- ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
- ! end speed-up
- St(1)=St(1)-PERC
-
- PR=PR+PERC
-
-! Split of effective rainfall into the two routing components
- PRHU1=PR*B
- PRHU2=PR*(1.-B)
-
-! Convolution of unit hydrograph UH1
- DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
- StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
- ENDDO
- StUH1(NH)=OrdUH1(NH)*PRHU1
-
-! 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)*PRHU2
- ENDDO
- StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
-
-! Potential intercatchment semi-exchange
- ! speed-up
- Rr = St(2)/Param(3)
- EXCH=Param(2)*Rr*Rr*Rr*SQRT(Rr)
- ! EXCH=Param(2)*(X(1)/Param(3))**3.5
- ! end speed-up
-
-! Routing store
- AEXCH1=EXCH
- IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
- St(2)=St(2)+StUH1(1)+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((StUH2(1)+EXCH).LT.0.) AEXCH2=-StUH2(1)
- QD=MAX(0.d0,StUH2(1)+EXCH)
-
-! Total runoff
- Q=QR+QD
- IF(Q.LT.0.) Q=0.
-
-! Variables storage
- MISC( 1)=E ! PE ! observed potential evapotranspiration [mm/day]
- MISC( 2)=P1 ! Precip ! observed total precipitation [mm/day]
- MISC( 3)=St(1) ! Prod ! production store level (St(1)) [mm]
- MISC( 4)=PN ! Pn ! net rainfall [mm/day]
- MISC( 5)=PS ! Ps ! part of Ps filling the production store [mm/day]
- MISC( 6)=AE ! AE ! actual evapotranspiration [mm/day]
- MISC( 7)=PERC ! Perc ! percolation (PERC) [mm/day]
- MISC( 8)=PR ! PR ! PR=PN-PS+PERC [mm/day]
- MISC( 9)=StUH1(1) ! Q9 ! outflow from UH1 (Q9) [mm/day]
- MISC(10)=StUH2(1) ! Q1 ! outflow from UH2 (Q1) [mm/day]
- MISC(11)=St(2) ! Rout ! routing store level (St(2)) [mm]
- MISC(12)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
- MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
- MISC(14)=AEXCH2 ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
- MISC(15)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
- MISC(16)=QR ! QR ! outflow from routing store (QR) [mm/day]
- MISC(17)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/day]
- MISC(18)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
-
-
-
-
- ENDSUBROUTINE
-
-
diff --git a/src/frun_GR4J.f90 b/src/frun_GR4J.f90
new file mode 100644
index 0000000000000000000000000000000000000000..fa486f31dd290c4501ca8f55f5251c88d2eca985
--- /dev/null
+++ b/src/frun_GR4J.f90
@@ -0,0 +1,319 @@
+!------------------------------------------------------------------------------
+! Subroutines relative to the annual GR4J model
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : frun_GR4J.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: C. Perrin
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2000
+! Last modified: 25/11/2019
+!------------------------------------------------------------------------------
+! 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.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+! 1. frun_gr4j
+! 2. MOD_GR4J
+!------------------------------------------------------------------------------
+
+
+ SUBROUTINE frun_gr4j(LInputs,InputsPrecip,InputsPE,NParam,Param, &
+ NStates,StateStart,NOutputs,IndOutputs, &
+ Outputs,StateEnd)
+! Subroutine that initializes GR4J, get its parameters, performs the call
+! to the MOD_GR4J 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/day]
+! InputsPE ! Vector of real, input series of potential evapotranspiration (PE) [mm/day]
+! 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])
+! 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_gr4j
+
+
+ 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
+ integer, dimension(NOutputs), intent(in) :: IndOutputs
+ ! out
+ doubleprecision, dimension(NStates), intent(out) :: StateEnd
+ doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
+
+ !! locals
+ integer :: I,K
+ integer, parameter :: NH=20,NMISC=30
+ doubleprecision, dimension(2) :: St
+ doubleprecision, dimension(NH) :: StUH1, OrdUH1
+ doubleprecision, dimension(2*NH) :: StUH2, OrdUH2
+ doubleprecision, dimension(NMISC) :: MISC
+ doubleprecision :: D,P1,E,Q
+
+ !--------------------------------------------------------------
+ ! Initializations
+ !--------------------------------------------------------------
+
+ ! initialization of model states to zero
+ St=0.
+ StUH1=0.
+ StUH2=0.
+
+ ! initialization of model states using StateStart
+ St(1) = StateStart(1)
+ St(2) = StateStart(2)
+ DO I=1,NH
+ StUH1(I)=StateStart(7+I)
+ ENDDO
+ DO I=1,2*NH
+ StUH2(I)=StateStart(7+I+NH)
+ ENDDO
+
+ ! parameter values
+ ! Param(1) : production store capacity (X1 - PROD) [mm]
+ ! Param(2) : intercatchment exchange coefficient (X2 - CES) [mm/day]
+ ! Param(3) : routing store capacity (X3 - ROUT) [mm]
+ ! Param(4) : time constant of unit hydrograph (X4 - TB) [day]
+
+ ! computation of UH ordinates
+ OrdUH1 = 0.
+ OrdUH2 = 0.
+
+ D=2.5
+ CALL UH1(OrdUH1,Param(4),D)
+ CALL UH2(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_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,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)
+ DO K=1,NH
+ StateEnd(7+K)=StUH1(K)
+ ENDDO
+ DO K=1,2*NH
+ StateEnd(7+NH+K)=StUH2(K)
+ ENDDO
+
+ RETURN
+
+ ENDSUBROUTINE
+
+
+
+
+
+!################################################################################################################################
+
+
+
+
+!**********************************************************************
+ SUBROUTINE MOD_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,MISC)
+! Calculation of streamflow on a single time step (day) with the GR4J model
+! Inputs:
+! St Vector of real, model states in stores at the beginning of the time step [mm]
+! StUH1 Vector of real, model states in Unit Hydrograph 1 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]
+! OrdUH1 Vector of real, ordinates in UH1 [-]
+! OrdUH2 Vector of real, ordinates in UH2 [-]
+! Param Vector of real, model parameters [various units]
+! P1 Real, value of rainfall during the time step [mm/day]
+! E Real, value of potential evapotranspiration during the time step [mm/day]
+! Outputs:
+! St Vector of real, model states in stores at the end of the time step [mm]
+! StUH1 Vector of real, model states in Unit Hydrograph 1 at the end of the time step [mm]
+! StUH2 Vector of real, model states in Unit Hydrograph 2 at the end of the time step [mm]
+! Q Real, value of simulated flow at the catchment outlet for the time step [mm/day]
+! MISC Vector of real, model outputs for the time step [mm/day]
+!**********************************************************************
+ Implicit None
+
+ !! locals
+ integer, parameter :: NParam=4,NMISC=30,NH=20
+ doubleprecision :: A,EN,ER,PN,PR,PS,WS,tanHyp
+ doubleprecision :: PERC,PRHU1,PRHU2,EXCH,QR,QD
+ doubleprecision :: AE,AEXCH1,AEXCH2
+ integer :: K
+ doubleprecision, parameter :: B=0.9
+ doubleprecision, parameter :: stored_val=25.62890625
+ doubleprecision :: expWS, TWS, Sr, Rr ! speed-up
+
+ !! dummies
+ ! in
+ doubleprecision, dimension(NParam), intent(in) :: Param
+ doubleprecision, intent(in) :: P1,E
+ doubleprecision, dimension(NH), intent(inout) :: OrdUH1
+ doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
+ ! inout
+ doubleprecision, dimension(2), intent(inout) :: St
+ doubleprecision, dimension(NH), intent(inout) :: StUH1
+ 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(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
+ St(1)=St(1)-ER
+ PS=0.
+ PR=0.
+ ELSE
+ EN=0.
+ AE=E
+ 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
+
+! Percolation from production store
+ IF(St(1).LT.0.) St(1)=0.
+ ! speed-up
+ ! (9/4)**4 = 25.62890625 = 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)/(9./4.*Param(1)))**4.)**(-0.25))
+ ! end speed-up
+ St(1)=St(1)-PERC
+
+ PR=PR+PERC
+
+! Split of effective rainfall into the two routing components
+ PRHU1=PR*B
+ PRHU2=PR*(1.-B)
+
+! Convolution of unit hydrograph UH1
+ DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
+ StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
+ ENDDO
+ StUH1(NH)=OrdUH1(NH)*PRHU1
+
+! 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)*PRHU2
+ ENDDO
+ StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
+
+! Potential intercatchment semi-exchange
+ ! speed-up
+ Rr = St(2)/Param(3)
+ EXCH=Param(2)*Rr*Rr*Rr*SQRT(Rr)
+ ! EXCH=Param(2)*(X(1)/Param(3))**3.5
+ ! end speed-up
+
+! Routing store
+ AEXCH1=EXCH
+ IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
+ St(2)=St(2)+StUH1(1)+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((StUH2(1)+EXCH).LT.0.) AEXCH2=-StUH2(1)
+ QD=MAX(0.d0,StUH2(1)+EXCH)
+
+! Total runoff
+ Q=QR+QD
+ IF(Q.LT.0.) Q=0.
+
+! Variables storage
+ MISC( 1)=E ! PE ! observed potential evapotranspiration [mm/day]
+ MISC( 2)=P1 ! Precip ! observed total precipitation [mm/day]
+ MISC( 3)=St(1) ! Prod ! production store level (St(1)) [mm]
+ MISC( 4)=PN ! Pn ! net rainfall [mm/day]
+ MISC( 5)=PS ! Ps ! part of Ps filling the production store [mm/day]
+ MISC( 6)=AE ! AE ! actual evapotranspiration [mm/day]
+ MISC( 7)=PERC ! Perc ! percolation (PERC) [mm/day]
+ MISC( 8)=PR ! PR ! PR=PN-PS+PERC [mm/day]
+ MISC( 9)=StUH1(1) ! Q9 ! outflow from UH1 (Q9) [mm/day]
+ MISC(10)=StUH2(1) ! Q1 ! outflow from UH2 (Q1) [mm/day]
+ MISC(11)=St(2) ! Rout ! routing store level (St(2)) [mm]
+ MISC(12)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
+ MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
+ MISC(14)=AEXCH2 ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
+ MISC(15)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
+ MISC(16)=QR ! QR ! outflow from routing store (QR) [mm/day]
+ MISC(17)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/day]
+ MISC(18)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
+
+
+
+
+ END SUBROUTINE
+
+
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f
deleted file mode 100644
index c9be386bafb934b81342e8c5b427a3708c37adc8..0000000000000000000000000000000000000000
--- a/src/frun_GR5J.f
+++ /dev/null
@@ -1,262 +0,0 @@
-
-
- SUBROUTINE frun_gr5j(
- ! inputs
- & LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm/day]
- & InputsPE , ! [double] input series of potential evapotranspiration (PE) [mm/day]
- & NParam , ! [integer] number of model parameters
- & Param , ! [double] parameter set
- & NStates , ! [integer] number of state variables
- & StateStart , ! [double] state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
- & NOutputs , ! [integer] number of output series
- & IndOutputs , ! [integer] indices of output series
- ! outputs
- & Outputs , ! [double] output series
- & StateEnd ) ! [double] state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
-
-
- !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr5j
-
-
- Implicit None
- ! input and output variables
- integer, intent(in) :: LInputs,NParam,NStates,NOutputs
- doubleprecision, dimension(LInputs) :: InputsPrecip
- doubleprecision, dimension(LInputs) :: InputsPE
- doubleprecision, dimension(NParam) :: Param
- doubleprecision, dimension(NStates) :: StateStart
- doubleprecision, dimension(NStates) :: StateEnd
- integer, dimension(NOutputs) :: IndOutputs
- doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
- ! parameters, internal states and variables
- integer NH,NMISC
- parameter (NH=20,NMISC=30)
- doubleprecision St(2), StUH2(2*NH)
- doubleprecision OrdUH2(2*NH)
- doubleprecision MISC(NMISC)
- doubleprecision D
- doubleprecision P1,E,Q
- integer I,K
-
- !--------------------------------------------------------------
- ! 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)
-
- 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/d]
- ! Param(3) : routing store capacity (X3 - ROUT) [mm]
- ! Param(4) : time constant of unit hydrograph (X4 - TB) [day]
- ! Param(5) : intercatchment exchange threshold (X5 - CES2) [-]
-
- !computation of HU ordinates
- OrdUH2 = 0.
-
- D=2.5
- CALL UH2(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_GR5J(St,StUH2,OrdUH2,Param,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)
- DO K=1,2*NH
- StateEnd(7+NH+K)=StUH2(K)
- ENDDO
-
- RETURN
-
- ENDSUBROUTINE
-
-
-
-
-
-!################################################################################################################################
-
-
-
-
-!**********************************************************************
- SUBROUTINE MOD_GR5J(St,StUH2,OrdUH2,Param,P1,E,Q,
- &MISC)
-! Run on a single time step with the GR5J model
-! Inputs:
-! St Vector of model states in stores at the beginning of the time step [mm]
-! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-! OrdUH2 Vector of ordinates in UH2 [-]
-! Param Vector of model parameters [various units]
-! P1 Value of rainfall during the time step [mm]
-! E Value of potential evapotranspiration during the time step [mm]
-! Outputs:
-! St Vector of model states in stores at the end of the time step [mm]
-! StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
-! Q Value of simulated flow at the catchment outlet for the time step [mm/day]
-! MISC Vector of model outputs for the time step [mm or mm/day]
-!**********************************************************************
- Implicit None
- INTEGER NH,NMISC,NParam
- PARAMETER (NH=20,NMISC=30)
- PARAMETER (NParam=5)
- DOUBLEPRECISION St(2),StUH2(2*NH)
- DOUBLEPRECISION OrdUH2(2*NH)
- DOUBLEPRECISION Param(NParam)
- DOUBLEPRECISION MISC(NMISC)
- DOUBLEPRECISION P1,E,Q
- DOUBLEPRECISION A,B,EN,ER,PN,PR,PS,WS,tanHyp
- DOUBLEPRECISION PERC,EXCH,QR,QD,Q1,Q9
- DOUBLEPRECISION AE,AEXCH1,AEXCH2
- INTEGER K
-
- DATA B/0.9/
- DOUBLEPRECISION TWS, Sr, Rr ! speed-up
-
- A=Param(1)
-
-
-! Interception and production store
- IF(P1.LE.E) THEN
- EN=E-P1
- PN=0.
- WS=EN/A
- IF(WS.GT.13.)WS=13.
- ! speed-up
- TWS = tanHyp(WS)
- 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
- St(1)=St(1)-ER
- PS=0.
- PR=0.
- ELSE
- EN=0.
- AE=E
- PN=P1-E
- WS=PN/A
- IF(WS.GT.13.)WS=13.
- ! speed-up
- TWS = tanHyp(WS)
- 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
-
-! Percolation from production store
- IF(St(1).LT.0.)St(1)=0.
-
- ! speed-up
- ! (9/4)**4 = 25.62890625
- Sr = St(1)/Param(1)
- Sr = Sr * Sr
- Sr = Sr * Sr
- PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
- ! PERC=X(2)*(1.-(1.+(X(2)/(9./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/day]
- MISC( 2)=P1 ! Precip ! observed total precipitation [mm/day]
- MISC( 3)=St(1) ! Prod ! production store level (St(1)) [mm]
- MISC( 4)=PN ! Pn ! net rainfall [mm/day]
- MISC( 5)=PS ! Ps ! part of Ps filling the production store [mm/day]
- MISC( 6)=AE ! AE ! actual evapotranspiration [mm/day]
- MISC( 7)=PERC ! Perc ! percolation (PERC) [mm/day]
- MISC( 8)=PR ! PR ! PR=PN-PS+PERC [mm/day]
- MISC( 9)=Q9 ! Q9 ! outflow from UH1 (Q9) [mm/day]
- MISC(10)=Q1 ! Q1 ! outflow from UH2 (Q1) [mm/day]
- MISC(11)=St(2) ! Rout ! routing store level (St(2)) [mm]
- MISC(12)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
- MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
- MISC(14)=AEXCH2 ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
- MISC(15)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
- MISC(16)=QR ! QR ! outflow from routing store (QR) [mm/day]
- MISC(17)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/day]
- MISC(18)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
-
-
- ENDSUBROUTINE
-
-
diff --git a/src/frun_GR5J.f90 b/src/frun_GR5J.f90
new file mode 100644
index 0000000000000000000000000000000000000000..d7b12ba03a16c445711b34fbaed6bb3261bcd266
--- /dev/null
+++ b/src/frun_GR5J.f90
@@ -0,0 +1,304 @@
+!------------------------------------------------------------------------------
+! Subroutines relative to the annual GR5J model
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : frun_GR5J.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: N. Le Moine
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2006
+! Last modified: 25/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! 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_gr5j
+! 2. MOD_GR5J
+!------------------------------------------------------------------------------
+
+
+ SUBROUTINE frun_gr5j(LInputs,InputsPrecip,InputsPE,NParam,Param, &
+ NStates,StateStart,NOutputs,IndOutputs, &
+ Outputs,StateEnd)
+! Subroutine that initializes GR5J, get its parameters, performs the call
+! to the MOD_GR5J 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/day]
+! InputsPE ! Vector of real, input series of potential evapotranspiration (PE) [mm/day]
+! 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])
+! 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_gr5j
+
+
+ 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
+ integer, dimension(NOutputs), intent(in) :: IndOutputs
+ ! out
+ doubleprecision, dimension(NStates), intent(out) :: StateEnd
+ doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
+
+ !! locals
+ integer :: I,K
+ integer, parameter :: NH=20,NMISC=30
+ doubleprecision, dimension(2) :: St
+ doubleprecision, dimension(2*NH) :: StUH2, OrdUH2
+ doubleprecision, dimension(NMISC) :: MISC
+ doubleprecision :: D,P1,E,Q
+
+ !--------------------------------------------------------------
+ ! 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)
+
+ 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/d]
+ ! Param(3) : routing store capacity (X3 - ROUT) [mm]
+ ! Param(4) : time constant of unit hydrograph (X4 - TB) [day]
+ ! Param(5) : intercatchment exchange threshold (X5 - CES2) [-]
+
+ !computation of HU ordinates
+ OrdUH2 = 0.
+
+ D=2.5
+ CALL UH2(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_GR5J(St,StUH2,OrdUH2,Param,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)
+ DO K=1,2*NH
+ StateEnd(7+NH+K)=StUH2(K)
+ ENDDO
+
+ RETURN
+
+ ENDSUBROUTINE
+
+
+
+
+
+!################################################################################################################################
+
+
+
+
+!**********************************************************************
+ SUBROUTINE MOD_GR5J(St,StUH2,OrdUH2,Param,P1,E,Q,MISC)
+! Calculation of streamflow on a single time step (day) with the GR5J 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]
+! OrdUH2 Vector of real, ordinates in UH2 [-]
+! Param Vector of real, model parameters [various units]
+! P1 Real, value of rainfall during the time step [mm]
+! E Real, value of potential evapotranspiration during the time step [mm]
+! Outputs:
+! St Vector of real, model states in stores at the end of the time step [mm]
+! StUH2 Vector of real, model states in Unit Hydrograph 2 at the end of the time step [mm]
+! Q Real, value of simulated flow at the catchment outlet for the time step [mm/day]
+! MISC Vector of real, model outputs for the time step [mm or mm/day]
+!**********************************************************************
+ Implicit None
+
+ !! locals
+ integer, parameter :: NParam=5,NMISC=30,NH=20
+ doubleprecision :: A,EN,ER,PN,PR,PS,WS,tanHyp
+ doubleprecision :: PERC,EXCH,QR,QD,Q1,Q9
+ doubleprecision :: AE,AEXCH1,AEXCH2
+ integer :: K
+ doubleprecision, parameter :: B=0.9
+ doubleprecision, parameter :: stored_val=25.62890625
+ doubleprecision :: expWS, TWS, Sr, Rr ! speed-up
+
+ !! dummies
+ ! in
+ doubleprecision, dimension(NParam), intent(in) :: Param
+ doubleprecision, intent(in) :: P1,E
+ doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
+ ! inout
+ doubleprecision, dimension(2), 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(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
+ St(1)=St(1)-ER
+ PS=0.
+ PR=0.
+ ELSE
+ EN=0.
+ AE=E
+ 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
+
+! Percolation from production store
+ IF(St(1).LT.0.) St(1)=0.
+
+ ! speed-up
+ ! (9/4)**4 = 25.62890625 = 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)/(9./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/day]
+ MISC( 2)=P1 ! Precip ! observed total precipitation [mm/day]
+ MISC( 3)=St(1) ! Prod ! production store level (St(1)) [mm]
+ MISC( 4)=PN ! Pn ! net rainfall [mm/day]
+ MISC( 5)=PS ! Ps ! part of Ps filling the production store [mm/day]
+ MISC( 6)=AE ! AE ! actual evapotranspiration [mm/day]
+ MISC( 7)=PERC ! Perc ! percolation (PERC) [mm/day]
+ MISC( 8)=PR ! PR ! PR=PN-PS+PERC [mm/day]
+ MISC( 9)=Q9 ! Q9 ! outflow from UH1 (Q9) [mm/day]
+ MISC(10)=Q1 ! Q1 ! outflow from UH2 (Q1) [mm/day]
+ MISC(11)=St(2) ! Rout ! routing store level (St(2)) [mm]
+ MISC(12)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
+ MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
+ MISC(14)=AEXCH2 ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
+ MISC(15)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
+ MISC(16)=QR ! QR ! outflow from routing store (QR) [mm/day]
+ MISC(17)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/day]
+ MISC(18)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
+
+
+ END SUBROUTINE
+
+
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f90
similarity index 52%
rename from src/frun_GR6J.f
rename to src/frun_GR6J.f90
index 17927fd4dd436a89c847d1dc55311b7a0f40d935..751d9114b19b7d2342bb84754327db83bbb53808 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f90
@@ -1,44 +1,75 @@
-
-
- SUBROUTINE frun_gr6j(
- ! inputs
- & LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm/day]
- & InputsPE , ! [double] input series of potential evapotranspiration (PE) [mm/day]
- & NParam , ! [integer] number of model parameters
- & Param , ! [double] parameter set
- & NStates , ! [integer] number of state variables
- & StateStart , ! [double] state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
- & NOutputs , ! [integer] number of output series
- & IndOutputs , ! [integer] indices of output series
- ! outputs
- & Outputs , ! [double] output series
- & StateEnd ) ! [double] state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
+!------------------------------------------------------------------------------
+! Subroutines relative to the annual GR6J model
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : frun_GR6J.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: R. Pushpalatha
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2010
+! Last modified: 25/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! 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_gr6j
+! 2. MOD_GR6J
+!------------------------------------------------------------------------------
+
+
+ SUBROUTINE frun_gr6j(LInputs,InputsPrecip,InputsPE,NParam,Param, &
+ NStates,StateStart,NOutputs,IndOutputs, &
+ Outputs,StateEnd)
+! Subroutine that initializes GR6J, get its parameters, performs the call
+! to the MOD_GR6J 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/day]
+! InputsPE ! Vector of real, input series of potential evapotranspiration (PE) [mm/day]
+! 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])
+! 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_gr6j
Implicit None
- ! input and output variables
+
+ !! dummies
+ ! in
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
- doubleprecision, dimension(LInputs) :: InputsPrecip
- doubleprecision, dimension(LInputs) :: InputsPE
- doubleprecision, dimension(NParam) :: Param
- doubleprecision, dimension(NStates) :: StateStart
- doubleprecision, dimension(NStates) :: StateEnd
- integer, dimension(NOutputs) :: IndOutputs
- doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
- ! parameters, internal states and variables
- integer NH,NMISC
- parameter (NH=20,NMISC=30)
- doubleprecision St(3), StUH1(NH), StUH2(2*NH)
- doubleprecision OrdUH1(NH), OrdUH2(2*NH)
- doubleprecision MISC(NMISC)
- doubleprecision D
- doubleprecision P1,E,Q
- integer I,K
+ doubleprecision, dimension(LInputs), intent(in) :: InputsPrecip
+ doubleprecision, dimension(LInputs), intent(in) :: InputsPE
+ doubleprecision, dimension(NParam), intent(in) :: Param
+ doubleprecision, dimension(NStates), intent(in) :: StateStart
+ integer, dimension(NOutputs), intent(in) :: IndOutputs
+ ! out
+ doubleprecision, dimension(NStates), intent(out) :: StateEnd
+ doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
+
+ !! locals
+ integer :: I,K
+ integer, parameter :: NH=20,NMISC=30
+ doubleprecision, dimension(3) :: St
+ doubleprecision, dimension(NH) :: StUH1, OrdUH1
+ doubleprecision, dimension(2*NH) :: StUH2, OrdUH2
+ doubleprecision, dimension(NMISC) :: MISC
+ doubleprecision :: D,P1,E,Q
!--------------------------------------------------------------
! Initializations
@@ -124,42 +155,49 @@
!**********************************************************************
- SUBROUTINE MOD_GR6J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
- &MISC)
-! Run on a single time step with the GR6J model
+ SUBROUTINE MOD_GR6J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,MISC)
+! Calculation of streamflow on a single time step (day) with the GR6J model
! Inputs:
-! St Vector of model states in stores at the beginning of the time step [mm]
-! StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-! OrdUH1 Vector of ordinates in UH1 [-]
-! OrdUH2 Vector of ordinates in UH2 [-]
-! Param Vector of model parameters [various units]
-! P1 Value of rainfall during the time step [mm]
-! E Value of potential evapotranspiration during the time step [mm]
+! St Vector of real, model states in stores at the beginning of the time step [mm]
+! StUH1 Vector of real, model states in Unit Hydrograph 1 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]
+! OrdUH1 Vector of real, ordinates in UH1 [-]
+! OrdUH2 Vector of real, ordinates in UH2 [-]
+! Param Vector of real, model parameters [various units]
+! P1 Real, value of rainfall during the time step [mm]
+! E Real, value of potential evapotranspiration during the time step [mm]
! Outputs:
-! St Vector of model states in stores at the end of the time step [mm]
-! StUH1 Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
-! StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
-! Q Value of simulated flow at the catchment outlet for the time step [mm/day]
-! MISC Vector of model outputs for the time step [mm or mm/day]
+! St Vector of real, model states in stores at the end of the time step [mm]
+! StUH1 Vector of real, model states in Unit Hydrograph 1 at the end of the time step [mm]
+! StUH2 Vector of real, model states in Unit Hydrograph 2 at the end of the time step [mm]
+! Q Real, value of simulated flow at the catchment outlet for the time step [mm/day]
+! MISC Vector of real, model outputs for the time step [mm or mm/day]
!**********************************************************************
Implicit None
- INTEGER NH,NMISC,NParam
- PARAMETER (NH=20,NMISC=30)
- PARAMETER (NParam=6)
- DOUBLEPRECISION St(3)
- DOUBLEPRECISION StUH1(NH),StUH2(2*NH),OrdUH1(NH),OrdUH2(2*NH)
- DOUBLEPRECISION Param(NParam)
- DOUBLEPRECISION MISC(NMISC)
- DOUBLEPRECISION P1,E,Q
- DOUBLEPRECISION A,B,C,EN,ER,PN,PR,PS,WS,tanHyp,AR
- DOUBLEPRECISION PERC,PRUH1,PRUH2,EXCH,QR,QD,QRExp
- DOUBLEPRECISION AE,AEXCH1,AEXCH2
- INTEGER K
-
- DATA B/0.9/
- DATA C/0.4/
- DOUBLEPRECISION TWS, Sr, Rr ! speed-up
+
+ !! locals
+ integer, parameter :: NParam=6,NMISC=30,NH=20
+ doubleprecision :: A,EN,ER,PN,PR,PS,WS,tanHyp,AR
+ doubleprecision :: PERC,PRUH1,PRUH2,EXCH,QR,QD,QRExp
+ doubleprecision :: AE,AEXCH1,AEXCH2
+ integer :: K
+ doubleprecision, parameter :: B=0.9, C=0.4
+ doubleprecision, parameter :: stored_val=25.62890625
+ doubleprecision :: expWS, TWS, Sr, Rr ! speed-up
+
+ !! dummies
+ ! in
+ doubleprecision, dimension(NParam), intent(in) :: Param
+ doubleprecision, intent(in) :: P1,E
+ doubleprecision, dimension(NH), intent(inout) :: OrdUH1
+ doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
+ ! inout
+ doubleprecision, dimension(3), intent(inout) :: St
+ doubleprecision, dimension(NH), intent(inout) :: StUH1
+ doubleprecision, dimension(2*NH), intent(inout) :: StUH2
+ ! out
+ doubleprecision, intent(out) :: Q
+ doubleprecision, dimension(NMISC), intent(out) :: MISC
A=Param(1)
@@ -169,10 +207,11 @@
EN=E-P1
PN=0.
WS=EN/A
- IF(WS.GT.13)WS=13.
+ IF(WS.GT.13) WS=13.
! speed-up
- TWS = tanHyp(WS)
+ 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))
@@ -187,10 +226,11 @@
AE=E
PN=P1-E
WS=PN/A
- IF(WS.GT.13)WS=13.
+ IF(WS.GT.13) WS=13.
! speed-up
- TWS = tanHyp(WS)
+ 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))
@@ -201,13 +241,13 @@
ENDIF
! Percolation from production store
- IF(St(1).LT.0.)St(1)=0.
+ IF(St(1).LT.0.) St(1)=0.
! speed-up
- ! (9/4)**4 = 25.62890625
+ ! (9/4)**4 = 25.62890625 = stored_val
Sr = St(1)/Param(1)
Sr = Sr * Sr
Sr = Sr * Sr
- PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
+ PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/stored_val)))
! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
! end speed-up
@@ -221,13 +261,13 @@
! Convolution of unit hydrograph UH1
DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
- StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRUH1
+ StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRUH1
ENDDO
StUH1(NH)=OrdUH1(NH)*PRUH1
! 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)*PRUH2
+ StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRUH2
ENDDO
StUH2(2*NH)=OrdUH2(2*NH)*PRUH2
@@ -238,7 +278,7 @@
AEXCH1=EXCH
IF((St(2)+(1-C)*StUH1(1)+EXCH).LT.0) AEXCH1=-St(2)-(1-C)*StUH1(1)
St(2)=St(2)+(1-C)*StUH1(1)+EXCH
- IF(St(2).LT.0.)St(2)=0.
+ IF(St(2).LT.0.) St(2)=0.
! speed-up
Rr = St(2)/Param(3)
@@ -253,16 +293,16 @@
! Update of exponential store
St(3)=St(3)+C*StUH1(1)+EXCH
AR=St(3)/Param(6)
- IF(AR.GT.33.)AR=33.
- IF(AR.LT.-33.)AR=-33.
+ IF(AR.GT.33.) AR=33.
+ IF(AR.LT.-33.) AR=-33.
IF(AR.GT.7.)THEN
- QRExp=St(3)+Param(6)/EXP(AR)
+ QRExp=St(3)+Param(6)/EXP(AR)
GOTO 3
ENDIF
IF(AR.LT.-7.)THEN
- QRExp=Param(6)*EXP(AR)
+ QRExp=Param(6)*EXP(AR)
GOTO 3
ENDIF
@@ -303,6 +343,6 @@
MISC(20)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
- ENDSUBROUTINE
+ END SUBROUTINE
diff --git a/src/utils_D.f b/src/utils_D.f
deleted file mode 100644
index eaf0d6ff9cd49091dfc13bd34b4820721fa29d67..0000000000000000000000000000000000000000
--- a/src/utils_D.f
+++ /dev/null
@@ -1,114 +0,0 @@
-!**********************************************************************
- SUBROUTINE UH1(OrdUH1,C,D)
-! Computation of ordinates of GR unit hydrograph UH1 using successive differences on the S curve SS1
-! Inputs:
-! C: time constant
-! D: exponent
-! Outputs:
-! OrdUH1: NH ordinates of discrete hydrograph
-!**********************************************************************
- Implicit None
- INTEGER NH
- PARAMETER (NH=20)
- DOUBLEPRECISION OrdUH1(NH)
- DOUBLEPRECISION C,D,SS1
- INTEGER I
-
- DO I=1,NH
- OrdUH1(I)=SS1(I,C,D)-SS1(I-1,C,D)
- ENDDO
- ENDSUBROUTINE
-
-
-!**********************************************************************
- SUBROUTINE UH2(OrdUH2,C,D)
-! Computation of ordinates of GR unit hydrograph HU2 using successive differences on the S curve SS2
-! Inputs:
-! C: time constant
-! D: exponent
-! Outputs:
-! OrdUH2: 2*NH ordinates of discrete hydrograph
-!**********************************************************************
- Implicit None
- INTEGER NH
- PARAMETER (NH=20)
- DOUBLEPRECISION OrdUH2(2*NH)
- DOUBLEPRECISION C,D,SS2
- INTEGER I
-
- DO I =1,2*NH
- OrdUH2(I)=SS2(I,C,D)-SS2(I-1,C,D)
- ENDDO
- ENDSUBROUTINE
-
-
-!**********************************************************************
- FUNCTION SS1(I,C,D)
-! Values of the S curve (cumulative HU curve) of GR unit hydrograph UH1
-! Inputs:
-! C: time constant
-! D: exponent
-! I: time-step
-! Outputs:
-! SS1: Values of the S curve for I
-!**********************************************************************
- Implicit None
- DOUBLEPRECISION C,D,SS1
- INTEGER I,FI
-
- FI=I
- IF(FI.LE.0) THEN
- SS1=0.
- RETURN
- ENDIF
- IF(FI.LT.C) THEN
- SS1=(FI/C)**D
- RETURN
- ENDIF
- SS1=1.
- ENDFUNCTION
-
-
-!**********************************************************************
- FUNCTION SS2(I,C,D)
-! Values of the S curve (cumulative HU curve) of GR unit hydrograph UH2
-! Inputs:
-! C: time constant
-! D: exponent
-! I: time-step
-! Outputs:
-! SS2: Values of the S curve for I
-!**********************************************************************
- Implicit None
- DOUBLEPRECISION C,D,SS2
- INTEGER I,FI
-
- FI=I
- IF(FI.LE.0) THEN
- SS2=0.
- RETURN
- ENDIF
- IF(FI.LE.C) THEN
- SS2=0.5*(FI/C)**D
- RETURN
- ENDIF
- IF(FI.LT.2.*C) THEN
- SS2=1.-0.5*(2.-FI/C)**D
- RETURN
- ENDIF
- SS2=1.
- ENDFUNCTION
-
-
-!**********************************************************************
- FUNCTION tanHyp(Val)
-! Computation of hyperbolic tangent
-!**********************************************************************
- Implicit None
- DOUBLEPRECISION Val,ValExp,tanHyp
-
- ValExp=EXP(Val)
- tanHyp=(ValExp - 1./ValExp)/(ValExp + 1./ValExp)
- RETURN
- ENDFUNCTION
-
diff --git a/src/utils_D.f90 b/src/utils_D.f90
new file mode 100644
index 0000000000000000000000000000000000000000..d12b96d3cc119e36bad19395595f86e68544a03e
--- /dev/null
+++ b/src/utils_D.f90
@@ -0,0 +1,191 @@
+!------------------------------------------------------------------------------
+! Subroutines relative to the daily unit hydrographs and hyperbolic
+! tangent calculation
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : utils_D.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: Unknown soldier
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2000
+! Last modified: 21/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+!
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+! 1. UH1
+! 2. UH2
+! 3. SS1
+! 4. SS2
+! 5. tanHyp
+!------------------------------------------------------------------------------
+
+
+!**********************************************************************
+ SUBROUTINE UH1(OrdUH1,C,D)
+! Subroutine that computes the ordinates of the daily GR unit hydrograph
+! UH1 using successive differences on the S curve SS1
+! Inputs
+! C ! Real, time constant
+! D ! Real, exponent
+! Outputs
+! OrdUH1: Vector of real, NH ordinates of the discrete hydrograph
+!**********************************************************************
+
+ Implicit None
+
+ !! locals
+ integer, parameter :: NH=20
+ doubleprecision :: SS1
+ integer :: I
+
+ !! dummies
+ ! in
+ double precision, intent(in) :: C,D
+ ! out
+ double precision, dimension (NH), intent(out) :: OrdUH1
+
+ DO I=1,NH
+ OrdUH1(I)=SS1(I,C,D)-SS1(I-1,C,D)
+ ENDDO
+ ENDSUBROUTINE
+
+
+!**********************************************************************
+ SUBROUTINE UH2(OrdUH2,C,D)
+! Subroutine that computes the ordinates of the daily GR unit hydrograph
+! UH2 using successive differences on the S curve SS2
+! Inputs
+! C ! Real, time constant
+! D ! Real, exponent
+! Outputs
+! OrdUH2: Vector of real, 2*NH ordinates of the discrete hydrograph
+!**********************************************************************
+
+ Implicit None
+
+ !! locals
+ integer, parameter :: NH=20
+ doubleprecision :: SS2
+ integer :: I
+
+ !! dummies
+ ! in
+ double precision, intent(in) :: C,D
+ ! out
+ double precision, dimension (2*NH), intent(out) :: OrdUH2
+
+ DO I =1,2*NH
+ OrdUH2(I)=SS2(I,C,D)-SS2(I-1,C,D)
+ ENDDO
+ ENDSUBROUTINE
+
+
+!**********************************************************************
+ FUNCTION SS1(I,C,D)
+! Function that computes the values of the S curve (cumulative UH
+! curve) of GR unit hydrograph UH1
+! Inputs
+! C ! Real, time constant
+! D ! Real, exponent
+! I ! Integer, time-step
+! Outputs
+! SS1 ! Real, value of the S curve for I
+!**********************************************************************
+
+ Implicit None
+
+ !! dummies
+ ! in
+ doubleprecision, intent(in) :: C,D
+ integer, intent(in) :: I
+ ! out
+ doubleprecision ::SS1
+
+ !! locals
+ integer :: FI
+
+ FI=I
+ IF(FI.LE.0) THEN
+ SS1=0.
+ RETURN
+ ENDIF
+ IF(FI.LT.C) THEN
+ SS1=(FI/C)**D
+ RETURN
+ ENDIF
+ SS1=1.
+ ENDFUNCTION
+
+
+!**********************************************************************
+ FUNCTION SS2(I,C,D)
+! Function that computes the values of the S curve (cumulative UH
+! curve) of GR unit hydrograph UH2
+! Inputs
+! C ! Real, time constant
+! D ! Real, exponent
+! I ! Integer, time-step
+! Outputs
+! SS2 ! Real, value of the S curve for I
+!**********************************************************************
+
+ Implicit None
+
+ !! dummies
+ ! in
+ doubleprecision, intent(in) :: C,D
+ integer, intent(in) :: I
+ ! out
+ doubleprecision ::SS2
+
+ !! locals
+ integer :: FI
+
+ FI=I
+ IF(FI.LE.0) THEN
+ SS2=0.
+ RETURN
+ ENDIF
+ IF(FI.LE.C) THEN
+ SS2=0.5*(FI/C)**D
+ RETURN
+ ENDIF
+ IF(FI.LT.2.*C) THEN
+ SS2=1.-0.5*(2.-FI/C)**D
+ RETURN
+ ENDIF
+ SS2=1.
+ ENDFUNCTION
+
+
+!**********************************************************************
+ FUNCTION tanHyp(Val)
+! Function that calculates the hyperbolic tangent
+!**********************************************************************
+! Inputs
+! Val ! Real, value
+! Outputs
+! tanHyp ! Real, value of the hyperbolic tangent
+
+ Implicit None
+
+ !! dummies
+ ! in
+ doubleprecision, intent(in) :: Val
+ ! out
+ doubleprecision :: tanHyp
+
+ !! locals
+ doubleprecision :: ValExp
+
+ ValExp=EXP(Val)
+ tanHyp=(ValExp - 1./ValExp)/(ValExp + 1./ValExp)
+ RETURN
+ ENDFUNCTION
+
diff --git a/src/utils_H.f b/src/utils_H.f
deleted file mode 100644
index e3c2ada84d03cb62a29e943aa2b5b94bbf00e256..0000000000000000000000000000000000000000
--- a/src/utils_H.f
+++ /dev/null
@@ -1,103 +0,0 @@
-
-
-!**********************************************************************
- SUBROUTINE UH1_H(OrdUH1,C,D)
-! Computation of ordinates of GR unit hydrograph UH2 using successive differences on the S curve SS1
-! Inputs:
-! C: time constant
-! D: exponent
-! Outputs:
-! OrdUH1: NH ordinates of discrete hydrograph
-!**********************************************************************
- Implicit None
- INTEGER NH
- PARAMETER (NH=480)
- DOUBLEPRECISION OrdUH1(NH)
- DOUBLEPRECISION C,D,SS1_H
- INTEGER I
-
- DO I=1,NH
- OrdUH1(I)=SS1_H(I,C,D)-SS1_H(I-1,C,D)
- ENDDO
- ENDSUBROUTINE
-
-
-!**********************************************************************
- SUBROUTINE UH2_H(OrdUH2,C,D)
-! Computation of ordinates of GR unit hydrograph UH2 using successive differences on the S curve SS2
-! Inputs:
-! C: time constant
-! D: exponent
-! Outputs:
-! OrdUH2: 2*NH ordinates of discrete hydrograph
-!**********************************************************************
- Implicit None
- INTEGER NH
- PARAMETER (NH=480)
- DOUBLEPRECISION OrdUH2(2*NH)
- DOUBLEPRECISION C,D,SS2_H
- INTEGER I
-
- DO I =1,2*NH
- OrdUH2(I)=SS2_H(I,C,D)-SS2_H(I-1,C,D)
- ENDDO
- ENDSUBROUTINE
-
-
-!**********************************************************************
- FUNCTION SS1_H(I,C,D)
-! Values of the S curve (cumulative HU curve) of GR unit hydrograph HU1
-! Inputs:
-! C: time constant
-! D: exponent
-! I: time-step
-! Outputs:
-! SS1: Values of the S curve for I
-!**********************************************************************
- Implicit None
- DOUBLEPRECISION C,D,SS1_H
- INTEGER I,FI
-
- FI=I
- IF(FI.LE.0) THEN
- SS1_H=0.
- RETURN
- ENDIF
- IF(FI.LT.C) THEN
- SS1_H=(FI/C)**D
- RETURN
- ENDIF
- SS1_H=1.
- ENDFUNCTION
-
-
-!**********************************************************************
- FUNCTION SS2_H(I,C,D)
-! Values of the S curve (cumulative HU curve) of GR unit hydrograph HU2
-! Inputs:
-! C: time constant
-! D: exponent
-! I: time-step
-! Outputs:
-! SS2: Values of the S curve for I
-!**********************************************************************
- Implicit None
- DOUBLEPRECISION C,D,SS2_H
- INTEGER I,FI
-
- FI=I
- IF(FI.LE.0) THEN
- SS2_H=0.
- RETURN
- ENDIF
- IF(FI.LE.C) THEN
- SS2_H=0.5*(FI/C)**D
- RETURN
- ENDIF
- IF(FI.LT.2.*C) THEN
- SS2_H=1.-0.5*(2.-FI/C)**D
- RETURN
- ENDIF
- SS2_H=1.
- ENDFUNCTION
-
diff --git a/src/utils_H.f90 b/src/utils_H.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6d9021392dc6c6fbd8119a0532539036d75ab111
--- /dev/null
+++ b/src/utils_H.f90
@@ -0,0 +1,163 @@
+!------------------------------------------------------------------------------
+! Subroutines relative to the hourly unit hydrographs
+!------------------------------------------------------------------------------
+! TITLE : airGR
+! PROJECT : airGR
+! FILE : utils_H.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: Unknown soldier
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2000
+! Last modified: 22/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+!
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+! 1. UH1_H
+! 2. UH2_H
+! 3. SS1_H
+! 4. SS2_H
+!------------------------------------------------------------------------------
+
+
+!**********************************************************************
+ SUBROUTINE UH1_H(OrdUH1,C,D)
+! Subroutine that computes the ordinates of the hourly GR unit hydrograph
+! UH1 using successive differences on the S curve SS1
+! Inputs
+! C ! Real, time constant
+! D ! Real, exponent
+! Outputs
+! OrdUH1: Vector of real, NH ordinates of the discrete hydrograph
+!**********************************************************************
+
+ Implicit None
+
+ !! locals
+ integer, parameter :: NH=480
+ doubleprecision :: SS1_H
+ integer :: I
+
+ !! dummies
+ ! in
+ double precision, intent(in) :: C,D
+ ! out
+ double precision, dimension (NH), intent(out) :: OrdUH1
+
+ DO I=1,NH
+ OrdUH1(I)=SS1_H(I,C,D)-SS1_H(I-1,C,D)
+ ENDDO
+ ENDSUBROUTINE
+
+
+!**********************************************************************
+ SUBROUTINE UH2_H(OrdUH2,C,D)
+! Subroutine that computes the ordinates of the hourly GR unit hydrograph
+! UH2 using successive differences on the S curve SS2
+! Inputs
+! C ! Real, time constant
+! D ! Real, exponent
+! Outputs
+! OrdUH2: Vector of real, 2*NH ordinates of the discrete hydrograph
+!**********************************************************************
+
+ Implicit None
+
+ !! locals
+ integer, parameter :: NH=480
+ doubleprecision :: SS2_H
+ integer :: I
+
+ !! dummies
+ ! in
+ double precision, intent(in) :: C,D
+ ! out
+ double precision, dimension (2*NH), intent(out) :: OrdUH2
+
+ DO I =1,2*NH
+ OrdUH2(I)=SS2_H(I,C,D)-SS2_H(I-1,C,D)
+ ENDDO
+ ENDSUBROUTINE
+
+
+!**********************************************************************
+ FUNCTION SS1_H(I,C,D)
+! Function that computes the values of the S curve (cumulative UH
+! curve) of GR unit hydrograph UH1
+! Inputs
+! C ! Real, time constant
+! D ! Real, exponent
+! I ! Integer, time-step
+! Outputs
+! SS1_H ! Real, value of the S curve for I
+!**********************************************************************
+
+ Implicit None
+
+ !! dummies
+ ! in
+ doubleprecision, intent(in) :: C,D
+ integer, intent(in) :: I
+ ! out
+ doubleprecision ::SS1_H
+
+ !! locals
+ integer :: FI
+
+ FI=I
+ IF(FI.LE.0) THEN
+ SS1_H=0.
+ RETURN
+ ENDIF
+ IF(FI.LT.C) THEN
+ SS1_H=(FI/C)**D
+ RETURN
+ ENDIF
+ SS1_H=1.
+ ENDFUNCTION
+
+
+!**********************************************************************
+ FUNCTION SS2_H(I,C,D)
+! Function that computes the values of the S curve (cumulative UH
+! curve) of GR unit hydrograph UH2
+! Inputs
+! C ! Real, time constant
+! D ! Real, exponent
+! I ! Integer, time-step
+! Outputs
+! SS2_H ! Real, value of the S curve for I
+!**********************************************************************
+
+ Implicit None
+
+ !! dummies
+ ! in
+ doubleprecision, intent(in) :: C,D
+ integer, intent(in) :: I
+ ! out
+ doubleprecision ::SS2_H
+
+ !! locals
+ integer :: FI
+
+ FI=I
+ IF(FI.LE.0) THEN
+ SS2_H=0.
+ RETURN
+ ENDIF
+ IF(FI.LE.C) THEN
+ SS2_H=0.5*(FI/C)**D
+ RETURN
+ ENDIF
+ IF(FI.LT.2.*C) THEN
+ SS2_H=1.-0.5*(2.-FI/C)**D
+ RETURN
+ ENDIF
+ SS2_H=1.
+ ENDFUNCTION
+