Commit 2d9f5f0a authored by Delaigue Olivier's avatar Delaigue Olivier

Merge branch 'cleanFortran' into 'master'

Clean fortran

See merge request !2
parents e0ad6416 bc912126
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"),
......
......@@ -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_&#42;GR&#42;()</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_&#42;GR&#42;()</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&#42;()</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
......
......@@ -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)
}
......@@ -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
......
man/figures/diagramGR4J-EN.png

44.4 KB | W: | H:

man/figures/diagramGR4J-EN.png

44.6 KB | W: | H:

man/figures/diagramGR4J-EN.png
man/figures/diagramGR4J-EN.png
man/figures/diagramGR4J-EN.png
man/figures/diagramGR4J-EN.png
  • 2-up
  • Swipe
  • Onion skin
man/figures/diagramGR5J-EN.png

39.5 KB | W: | H:

man/figures/diagramGR5J-EN.png

39.8 KB | W: | H:

man/figures/diagramGR5J-EN.png
man/figures/diagramGR5J-EN.png
man/figures/diagramGR5J-EN.png
man/figures/diagramGR5J-EN.png
  • 2-up
  • Swipe
  • Onion skin
man/figures/diagramGR6J-EN.png

55.3 KB | W: | H:

man/figures/diagramGR6J-EN.png

55.7 KB | W: | H:

man/figures/diagramGR6J-EN.png
man/figures/diagramGR6J-EN.png
man/figures/diagramGR6J-EN.png
man/figures/diagramGR6J-EN.png
  • 2-up
  • Swipe
  • Onion skin
!------------------------------------------------------------------------------
! 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
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)
!------------------------------------------------------------------------------
! 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
! 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 NMISC
parameter (NMISC=3)
doubleprecision MISC(NMISC)
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
integer I,K
!--------------------------------------------------------------
! Initializations
......@@ -49,8 +79,8 @@
MISC = -999.999
! Outputs = -999.999 ! initialization made in R
StateStart = -999.999 ! CRAN-compatibility updates
StateEnd = -999.999 ! CRAN-compatibility updates
! StateStart = -999.999 ! CRAN-compatibility updates
! StateEnd = -999.999 ! CRAN-compatibility updates
!--------------------------------------------------------------
......@@ -66,7 +96,7 @@
CALL MOD_GR1A(Param,P0,P1,E1,Q,MISC)
! storage of outputs
DO I=1,NOutputs
Outputs(k,I)=MISC(IndOutputs(I))
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
......@@ -85,25 +115,31 @@
!**********************************************************************
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]
! 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
INTEGER NMISC,NParam
PARAMETER (NMISC=3)
PARAMETER (NParam=1)
DOUBLEPRECISION Param(NParam)
DOUBLEPRECISION MISC(NMISC)
DOUBLEPRECISION P0,P1,E1,Q
DOUBLEPRECISION tt ! speed-up
!! 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
......@@ -118,6 +154,6 @@
MISC( 3)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/year]
ENDSUBROUTINE
END SUBROUTINE
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])
!------------------------------------------------------------------------------
! 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
! 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 NMISC
parameter (NMISC=30)
doubleprecision St(2)
doubleprecision MISC(NMISC)
doubleprecision P,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 :: NMISC=30
doubleprecision, dimension(2) :: St
doubleprecision, dimension(NMISC) :: MISC
doubleprecision :: P,E,Q
!--------------------------------------------------------------
! Initializations
......@@ -97,44 +131,51 @@
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]
! 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 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]
! 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
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
!! 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.
IF(WS.GT.13.) WS=13.
! speed-up
TWS = tanHyp(WS)
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.
IF(WS.GT.13.) WS=13.
! speed-up
TWS = tanHyp(WS)
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
......@@ -178,6 +219,6 @@
MISC(10)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/month]
ENDSUBROUTINE
END SUBROUTINE
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
!**********************************************************************