Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 2558 additions and 1459 deletions
+2558 -1459
man/figures/diagramGR5J-EN.png
View replaced file @ 0967defb View file @ 4bccfee2
man/figures/diagramGR5J-EN.png

39 KB | W: 0px | H: 0px

man/figures/diagramGR5J-EN.png

39.8 KB | W: 0px | H: 0px

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.pdf
View file @ 4bccfee2
No preview for this file type
man/figures/diagramGR6J-EN.png
View replaced file @ 0967defb View file @ 4bccfee2
man/figures/diagramGR6J-EN.png

54.3 KB | W: 0px | H: 0px

man/figures/diagramGR6J-EN.png

55.7 KB | W: 0px | H: 0px

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
man/plot.OutputsModel.Rd
View file @ 4bccfee2
\encoding{UTF-8}
\name{plot.OutputsModel}
\name{plot}
\alias{plot.OutputsModel}
\alias{plot}
\alias{exampleSimPlot}
\alias{simGR4J}
\alias{simCNGR4J}
\title{Default preview of model outputs}
\description{
Function which creates a screen plot giving an overview of the model outputs.
}
\usage{
\method{plot}{OutputsModel}(x, Qobs = NULL, IndPeriod_Plot = NULL,
BasinArea = NULL, which = "all", log_scale = FALSE,
cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1, verbose = TRUE, ...)
BasinArea = NULL, which = "synth", log_scale = FALSE,
cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1,
AxisTS = function(x) axis.POSIXct(side = 1, x = x$DatesR, ...),
LayoutMat = NULL,
LayoutWidths = rep.int(1, ncol(LayoutMat)),
LayoutHeights = rep.int(1, nrow(LayoutMat)),
verbose = TRUE, ...)
}
\arguments{
\item{x}{[object of class \emph{OutputsModel}] list of model outputs (which must at least include DatesR, Precip and Qsim) [POSIXlt, mm, mm]}
\item{x}{[object of class \emph{OutputsModel}] list of model outputs (which must at least include DatesR, Precip and Qsim) [POSIXlt, mm/time step, mm/time step]}
\item{Qobs}{(optional) [numeric] time series of observed flow (for the same time steps than simulated) [mm/time step]}
......@@ -25,9 +38,9 @@
\item{BasinArea}{(optional) [numeric] basin area [km2], used to plot flow axes in m3/s}
\item{which}{(optional) [character] choice of plots \cr (e.g. c(\code{"Precip"}, \code{"Temp"}, \code{"SnowPack"}, \code{"Flows"}, \code{"Regime"}, \code{"CumFreq"}, \code{"CorQQ"})), default = \code{"all"}}
\item{which}{(optional) [character] choice of plots \cr (e.g. c(\code{"Precip"}, \code{"Temp"}, \code{"SnowPack"}, \code{"Flows"}, \code{"Error"}, \code{"Regime"}, \code{"CumFreq"}, \code{"CorQQ"})), default = \code{"synth"}, see details below}
\item{log_scale}{(optional) [boolean] indicating if the flow axis is to be logarithmic, default = \code{FALSE}}
\item{log_scale}{(optional) [boolean] indicating if the flow and the error time series axis and the flow error time series axis are to be logarithmic, default = \code{FALSE}}
\item{cex.axis}{(optional) [numeric] the magnification to be used for axis annotation relative to the current setting of \code{cex}}
......@@ -37,9 +50,17 @@
\item{lwd}{(optional) [numeric] the line width (a positive number)}
\item{AxisTS}{(optional) [function] to manage x-axis representing calendar dates and times on time series plots (see \code{\link{axis}} and \code{\link{axis.POSIXct}})}
\item{LayoutMat}{(optional) [numeric] a matrix object specifying the location of the next N figures on the output device. Each value in the matrix must be 0 or a positive integer. If N is the largest positive integer in the matrix, then the integers \eqn{1, \dots, N-1} must also appear at least once in the matrix (see \code{\link{layout}})}
\item{LayoutWidths}{(optional) [numeric] a vector of values for the widths of columns on the device (see \code{\link{layout}})}
\item{LayoutHeights}{(optional) [numeric] a vector of values for the heights of rows on the device (see \code{\link{layout}})}
\item{verbose}{(optional) [boolean] indicating if the function is run in verbose mode or not, default = \code{TRUE}}
\item{...}{other parameters to be passed through to plotting functions}
\item{...}{(optional) other parameters to be passed through to plotting functions}
}
......@@ -48,28 +69,92 @@ Screen plot window.
}
\description{
Function which creates a screen plot giving an overview of the model outputs.
\details{
Different types of independent graphs are available (depending on the model, but always drawn in this order):
\itemize{
\item \code{"Precip"}: time series of total precipitation
\item \code{"PotEvap"}: time series of potential evapotranspiration
\item \code{"ActuEvap"}: time series of simulated actual evapotranspiration (overlaid to \code{"PotEvap"} if already drawn)
\item \code{"Temp"}: time series of temperature (plotted only if CemaNeige is used)
\item \code{"SnowPack"}: time series of snow water equivalent (plotted only if CemaNeige is used)
\item \code{"Flows"}: time series of simulated flows (and observed flows if provided)
\item \code{"Error"}: time series of simulated flows minus observed flows (and observed flows if provided)
\item \code{"Regime"}: centred 30-day rolling mean applied on interannual average of daily simulated flows (and observed flows if provided)
\item \code{"CorQQ"}: correlation plot between simulated and observed flows (only if observed flows provided)
\item \code{"CumFreq"}: cumulative frequency plot for simulated flows (and observed flows if provided)
}
Different dashboards of results including various graphs are available:
\itemize{
\item \code{"perf"}: corresponds to \code{"Error"}, \code{"Regime"}, \code{"CumFreq"} and \code{"CorQQ"}
\item \code{"ts"}: corresponds to \code{"Precip"}, \code{"PotEvap"}, \code{"Temp"}, \code{"SnowPack"} and \code{"Flows"}
\item \code{"synth"}: corresponds to \code{"Precip"}, \code{"Temp"}, \code{"SnowPack"}, \code{"Flows"}, \code{"Regime"}, \code{"CumFreq"} and \code{"CorQQ"}
\item \code{"all"}: corresponds to \code{"Precip"}, \code{"PotEvap"}, \code{"ActuEvap"}, \code{"Temp"}, \code{"SnowPack"}, \code{"Flows"}, \code{"Error"}, \code{"Regime"}, \code{"CumFreq"} and \code{"CorQQ"}
}
\details{
Dashboard of results including various graphs (depending on the model):\cr
(1) time series of total precipitation\cr
(2) time series of temperature (plotted only if CemaNeige is used)\cr
(3) time series of snow pack (plotted only if CemaNeige is used)\cr
(4) time series of simulated flows (and observed flows if provided)\cr
(5) interannual median monthly simulated flow (and observed flows if provided)\cr
(6) correlation plot between simulated and observed flows (if observed flows provided)\cr
(7) cumulative frequency plot for simulated flows (and observed flows if provided)
If several dashboards are selected, or if an independent graph is called with a dashboard, the graphical device will include all the requested graphs without redundancy.
}
\author{
Laurent Coron, Olivier Delaigue, Guillaume Thirel
}
\examples{
### See examples of RunModel_GR4J or RunModel_CemaNeigeGR4J functions
### see examples of RunModel_GR4J or RunModel_CemaNeigeGR4J functions
### to understand how the example datasets have been prepared
## loading examples dataset for GR4J and GR4J + CemaNeige
data(exampleSimPlot)
### Qobs and outputs from GR4J and GR4J + CemaNeige models
str(simGR4J, max.level = 1)
str(simCNGR4J, max.level = 1)
### default dashboard (which = "synth")
## GR models whithout CemaNeige
plot(simGR4J$OutputsModel, Qobs = simGR4J$Qobs)
## GR models whith CemaNeige ("Temp" and "SnowPack" added)
plot(simCNGR4J$OutputsModel, Qobs = simCNGR4J$Qobs)
### "Error" and "CorQQ" plots cannot be display whithout Qobs
plot(simGR4J$OutputsModel, which = "all", Qobs = simGR4J$Qobs)
plot(simGR4J$OutputsModel, which = "all", Qobs = NULL)
### complex plot arrangements
plot(simGR4J$OutputsModel, Qobs = simGR4J$Qobs,
which = c("Flows", "Regime", "CumFreq", "CorQQ"),
LayoutMat = matrix(c(1, 2, 3, 1, 4, 4), ncol = 2),
LayoutWidths = c(1.5, 1),
LayoutHeights = c(0.5, 1, 1))
### customize x-axis on time series plot
FunAxisTS <- function(x) {
axis.POSIXct(side = 1, x = x$DatesR,
at = pretty(x$DatesR, n = 10),
format = "\%Y-\%m-\%d")
}
plot(simGR4J$OutputsModel, Qobs = simGR4J$Qobs, AxisTS = FunAxisTS)
### add a main title
## the whole list of settable par's
opar <- par(no.readonly = TRUE)
## define outer margins and a title inside it
par(oma = c(0, 0, 3, 0))
plot(simGR4J$OutputsModel, Qobs = simGR4J$Qobs)
title(main = "GR4J outputs", outer = TRUE, line = 1.2, cex.main = 1.4)
## reset original par
par(opar)
}
src/airGR.c
View file @ 4bccfee2
......@@ -2,7 +2,7 @@
#include <stdlib.h> // for NULL
#include <R_ext/Rdynload.h>
/* FIXME:
/* FIXME:
Check these declarations against the C/Fortran source code.
*/
......@@ -11,18 +11,22 @@ extern void F77_NAME(frun_cemaneige)(int *, double *, double *, double *, double
extern void F77_NAME(frun_gr1a)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr2m)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr4h)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr5h)(int *, double *, double *, int *, double *, int *, double *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr4j)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr5j)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_gr6j)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
extern void F77_NAME(frun_pe_oudin)(int *, double *, double *, double *, double *);
static const R_FortranMethodDef FortranEntries[] = {
{"frun_cemaneige", (DL_FUNC) &F77_NAME(frun_cemaneige), 14},
{"frun_gr1a", (DL_FUNC) &F77_NAME(frun_gr1a), 11},
{"frun_gr2m", (DL_FUNC) &F77_NAME(frun_gr2m), 11},
{"frun_gr4h", (DL_FUNC) &F77_NAME(frun_gr4h), 11},
{"frun_gr5h", (DL_FUNC) &F77_NAME(frun_gr5h), 12},
{"frun_gr4j", (DL_FUNC) &F77_NAME(frun_gr4j), 11},
{"frun_gr5j", (DL_FUNC) &F77_NAME(frun_gr5j), 11},
{"frun_gr6j", (DL_FUNC) &F77_NAME(frun_gr6j), 11},
{"frun_pe_oudin", (DL_FUNC) &F77_NAME(frun_pe_oudin), 5},
{NULL, NULL, 0}
};
......
src/frun_CEMANEIGE.f → src/frun_CEMANEIGE.f90
View file @ 4bccfee2
!------------------------------------------------------------------------------
! 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: Coron, L.
! Further cleaning: Thirel, G.
!------------------------------------------------------------------------------
! Creation date: 2011
! Last modified: 22/11/2019
!------------------------------------------------------------------------------
! REFERENCES
! Riboust, P., Thirel, G., Le Moine, N. and Ribstein, P. (2019). Revisiting a
! simple degree-day model for integrating satellite data: Implementation of
! SWE-SCA hystereses. Journal of Hydrology and Hydromechanics, 67(1), 70–81,
! doi: 10.2478/johh-2018-0004.
!
! Valéry, A., Andréassian, V. and Perrin, C. (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, 517(0), 1166-1175, doi: 10.1016/j.jhydrol.2014.04.059.
!
! Valéry, A., Andréassian, V. and Perrin, C. (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, 517(0), 1176-1187,
! 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]
& 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 initialising = 2
& StateStart , ! [double] state variables used when the model run starts
& IsHyst , ! [logical] 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
logical, intent(in) :: IsHyst ! 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
!! 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.
!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
!--------------------------------------------------------------
!Initialisations
! Initializations
!--------------------------------------------------------------
!initialisation of constants
! initialization of constants
Tmelt=0.
MinSpeed=0.1
!initialisation of model states using StateStart
! initialization of model states using StateStart
G=StateStart(1)
eTG=StateStart(2)
Gthreshold=StateStart(3)
Glocalmax=StateStart(4)
Gratio=0.
PliqAndMelt=0.
!setting parameter values
! setting parameter values
CTG=Param(1)
Kf=Param(2)
IF (IsHyst) THEN
IF (IsHystBool) THEN
Gthreshold=StateStart(3)
Glocalmax=StateStart(4)
Gacc=Param(3)
prct=Param(4)
......@@ -75,33 +120,36 @@
ELSE
Gthreshold=0.9*MeanAnSolidPrecip
Glocalmax=-999.999
Gacc=-999.999
prct=-999.999
ENDIF
!initialisation of model outputs
c StateEnd = -999.999 !initialisation made in R
c Outputs = -999.999 !initialisation made in R
! initialization of model outputs
! StateEnd = -999.999 !initialization made in R
! Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
!Time loop
! Time loop
!--------------------------------------------------------------
DO k=1,LInputs
!SolidPrecip and LiquidPrecip
! SolidPrecip and LiquidPrecip
Pliq=(1.-InputsFracSolidPrecip(k))*InputsPrecip(k)
Psol=InputsFracSolidPrecip(k)*InputsPrecip(k)
!Snow pack volume before melt
! Snow pack volume before melt
Ginit=G
G=G+Psol
!Snow pack thermal state before melt
! Snow pack thermal state before melt
eTG=CTG*eTG + (1.-CTG)*InputsTemp(k)
IF(eTG.GT.0.) eTG=0.
!Potential melt
! Potential melt
IF(eTG.EQ.0..AND.InputsTemp(k).GT.Tmelt) THEN
PotMelt=Kf*(InputsTemp(k)-Tmelt)
IF(PotMelt.GT.G) PotMelt=G
......@@ -109,10 +157,10 @@ c Outputs = -999.999 !initialisation made in R
PotMelt=0.
ENDIF
IF (IsHyst) THEN
IF (IsHystBool) THEN
IF (Potmelt.GT.0.) THEN
IF (G.LT.Glocalmax.AND.Gratio.EQ.1.) Glocalmax=G !Update in case of potential melt and G lower than Gseuil
IF (G.LT.Glocalmax.AND.Gratio.EQ.1.) Glocalmax=G ! Update in case of potential melt and G lower than Gseuil
Gratio=MIN(G/Glocalmax,1.d0)
ENDIF
ELSE
......@@ -123,20 +171,19 @@ c Outputs = -999.999 !initialisation made in R
ENDIF
ENDIF
!Actual melt
! Actual melt
Melt=((1.-MinSpeed)*Gratio+MinSpeed)*PotMelt
!Update of snow pack volume
! Update of snow pack volume
G=G-Melt
IF (IsHyst) THEN
dG=G-Ginit !Melt in case of negative dG, accumulation otherwise
IF (IsHystBool) THEN
dG=G-Ginit ! Melt in case of negative dG, accumulation otherwise
IF (dG.GT.0.) THEN
Gratio = MIN(Gratio+(Psol-Melt)/Gacc,1.d0) !Psol - Melt = dG
IF (Gratio.EQ.1.) Glocalmax = Gthreshold
ENDIF
IF (dG.LT.0.) THEN
ELSE
Gratio=MIN(G/Glocalmax,1.d0)
ENDIF
ELSE
......@@ -148,22 +195,22 @@ c Outputs = -999.999 !initialisation made in R
ENDIF
!Water volume to pass to the hydrological model
! Water volume to pass to the hydrological model
PliqAndMelt=Pliq+Melt
!Storage of outputs
! Storage of outputs
DO I=1,NOutputs
IF(IndOutputs(I).EQ.1) Outputs(k,I)=Pliq ! Pliq ! observed liquid precipitation [mm/day]
IF(IndOutputs(I).EQ.2) Outputs(k,I)=Psol ! Psol ! observed solid precipitation [mm/day]
IF(IndOutputs(I).EQ.3) Outputs(k,I)=G ! SnowPack ! snow pack [mm]
IF(IndOutputs(I).EQ.4) Outputs(k,I)=eTG ! ThermalState ! thermal state [°C]
IF(IndOutputs(I).EQ.5) Outputs(k,I)=Gratio ! Gratio ! Gratio [-]
IF(IndOutputs(I).EQ.6) Outputs(k,I)=PotMelt ! PotMelt ! potential snow melt [mm/day]
IF(IndOutputs(I).EQ.7) Outputs(k,I)=Melt ! Melt ! melt [mm/day]
IF(IndOutputs(I).EQ.8) Outputs(k,I)=PliqAndMelt ! PliqAndMelt ! liquid precipitation + melt [mm/day]
IF(IndOutputs(I).EQ.9) Outputs(k,I)=InputsTemp(k) ! Temp ! air temperature [°C]
IF(IndOutputs(I).EQ.10) Outputs(k,I)=Gthreshold ! Gthreshold ! melt threshold [mm]
IF(IndOutputs(I).EQ.11) Outputs(k,I)=Glocalmax ! Glocalmax ! local melt threshold for hysteresis [mm]
IF(IndOutputs(I).EQ.1) Outputs(k,I)=Pliq ! Pliq ! [numeric] observed liquid precipitation [mm/time step]
IF(IndOutputs(I).EQ.2) Outputs(k,I)=Psol ! Psol ! [numeric] observed solid precipitation [mm/time step]
IF(IndOutputs(I).EQ.3) Outputs(k,I)=G ! SnowPack ! [numeric] snow pack [mm]
IF(IndOutputs(I).EQ.4) Outputs(k,I)=eTG ! ThermalState ! [numeric] thermal state [°C]
IF(IndOutputs(I).EQ.5) Outputs(k,I)=Gratio ! Gratio ! [numeric] Gratio [-]
IF(IndOutputs(I).EQ.6) Outputs(k,I)=PotMelt ! PotMelt ! [numeric] potential snow melt [mm/time step]
IF(IndOutputs(I).EQ.7) Outputs(k,I)=Melt ! Melt ! [numeric] melt [mm/time step]
IF(IndOutputs(I).EQ.8) Outputs(k,I)=PliqAndMelt ! PliqAndMelt ! [numeric] liquid precipitation + melt [mm/time step]
IF(IndOutputs(I).EQ.9) Outputs(k,I)=InputsTemp(k) ! Temp ! [numeric] air temperature [°C]
IF(IndOutputs(I).EQ.10) Outputs(k,I)=Gthreshold ! Gthreshold ! [numeric] melt threshold [mm]
IF(IndOutputs(I).EQ.11) Outputs(k,I)=Glocalmax ! Glocalmax ! [numeric] local melt threshold for hysteresis [mm]
ENDDO
ENDDO
......@@ -175,5 +222,5 @@ c Outputs = -999.999 !initialisation made in R
RETURN
ENDSUBROUTINE
END SUBROUTINE
src/frun_GR1A.f deleted 100644 → 0
View file @ 0967defb
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
c 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)
c Q = -999.999
c 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
c################################################################################################################################
C**********************************************************************
SUBROUTINE MOD_GR1A(Param,P0,P1,E1,Q,MISC)
C Calculation of streamflow on a single time step with the GR1A model
C Inputs:
C Param Vector of model parameters (Param(1) [-])
C P0 Value of rainfall during the previous time step [mm/year]
C P1 Value of rainfall during the current time step [mm/year]
C E1 Value of potential evapotranspiration during the current time step [mm/year]
C Outputs:
C Q Value of simulated flow at the catchment outlet for the current time step [mm/year]
C MISC Vector of model outputs for the time step [mm/year]
C**********************************************************************
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
C 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)
! fin speed-up
C 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
src/frun_GR1A.f90 0 → 100644
View file @ 4bccfee2
!------------------------------------------------------------------------------
! Subroutines relative to the annual GR1A model
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_GR1A.f
!------------------------------------------------------------------------------
! AUTHORS
! Original code: Mouelhi, S.
! Cleaning and formatting for airGR: Coron, L.
! Further cleaning: Thirel, G.
!------------------------------------------------------------------------------
! 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)
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 = StateStart ! 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
src/frun_GR2M.f deleted 100644 → 0
View file @ 0967defb
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.
!initilisation 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
c StateEnd = -999.999 !initialization made in R
c Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
!Time loop
!--------------------------------------------------------------
DO k=1,LInputs
P =InputsPrecip(k)
E =InputsPE(k)
c Q = -999.999
c 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
c################################################################################################################################
C**********************************************************************
SUBROUTINE MOD_GR2M(St,Param,P,E,Q,MISC)
C Calculation of streamflow on a single time step (month) with the GR2M model
C Inputs:
C St Vector of model states at the beginning of the time step [mm]
C Param Vector of model parameters (Param(1) [mm]; Param(2) [-])
C P Value of rainfall during the time step [mm/month]
C E Value of potential evapotranspiration during the time step [mm/month]
C Outputs:
C St Vector of model states at the end of the time step [mm]
C Q Value of simulated flow at the catchment outlet for the time step [mm/month]
C MISC Vector of model outputs for the time step [mm/month]
C**********************************************************************
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
C 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))
! fin 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))
! fin speed-up
AE = S1 - S2
C 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.)
! fin speed-up
P2=S2-St(1)
P3=P1+P2
C QR calculation (routing store)
R1=St(2)+P3
C Water exchange
R2=Param(2)*R1
EXCH = R2 - R1
C Total runoff
Q=R2*R2/(R2+60.)
C Updating store level
St(2)=R2-Q
C 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
src/frun_GR2M.f90 0 → 100644
View file @ 4bccfee2
!------------------------------------------------------------------------------
! Subroutines relative to the annual GR2M model
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_GR2M.f
!------------------------------------------------------------------------------
! AUTHORS
! Original code: Mouelhi, S.
! Cleaning and formatting for airGR: Coron, L.
! Further cleaning: Thirel, G.
!------------------------------------------------------------------------------
! Creation date: 2003
! Last modified: 16/04/2020
!------------------------------------------------------------------------------
! 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., Michel, C., Perrin, C. and Andréassian, V. (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])
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,S1,S2
doubleprecision :: P1,P2,P3,R1,R2,AE,AEXCH,PS
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
PS = P - P1
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
AEXCH = 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)=PS ! Ps ! [numeric] part of P filling the production store [mm/month]
MISC( 6)=AE ! AE ! [numeric] actual evapotranspiration [mm/month]
MISC( 7)=P2 ! Perc ! [numeric] percolation (P2) [mm/month]
MISC( 8)=P3 ! PR ! [numeric] P3=P1+P2 [mm/month]
MISC( 9)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
MISC(10)=AEXCH ! AEXCH ! [numeric] actual groundwater exchange (AEXCH) [mm/month]
MISC(11)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/month]
END SUBROUTINE
src/frun_GR4H.f deleted 100644 → 0
View file @ 0967defb
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
!--------------------------------------------------------------
!initilization 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
c StateEnd = -999.999 !initialization made in R
c Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
!Time loop
!--------------------------------------------------------------
DO k=1,LInputs
P1=InputsPrecip(k)
E =InputsPE(k)
c Q = -999.999
c 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
c################################################################################################################################
C**********************************************************************
SUBROUTINE MOD_GR4H(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
&MISC)
C Run on a single time step with the GR4H model
C Inputs:
C St Vector of model states in stores at the beginning of the time step [mm]
C StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
C StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
C OrdUH1 Vector of ordinates in UH1 [-]
C OrdUH2 Vector of ordinates in UH2 [-]
C Param Vector of model parameters [various units]
C P1 Value of rainfall during the time step [mm]
C E Value of potential evapotranspiration during the time step [mm]
C Outputs:
C St Vector of model states in stores at the beginning of the time step [mm]
C StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
C StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
C Q Value of simulated flow at the catchment outlet for the time step [mm/hour]
C MISC Vector of model outputs for the time step [mm/hour]
C**********************************************************************
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)
C 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))
! fin 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))
! fin speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
C 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))
! fin speed-up
St(1)=St(1)-PERC
PR=PR+PERC
C Split of effective rainfall into the two routing components
PRHU1=PR*B
PRHU2=PR*(1.-B)
C 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
C 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
C 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
! fin speed-up
C 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.))
! fin speed-up
St(2)=St(2)-QR
C Runoff from direct branch QD
AEXCH2=EXCH
IF((StUH2(1)+EXCH).LT.0.) AEXCH2=-StUH2(1)
QD=MAX(0.d0,StUH2(1)+EXCH)
C Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
C 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
src/frun_GR4H.f90 0 → 100644
View file @ 4bccfee2
!------------------------------------------------------------------------------
! Subroutines relative to the annual GR4H model
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_GR4H.f
!------------------------------------------------------------------------------
! AUTHORS
! Original code: Perrin, C.
! Cleaning and formatting for airGR: Coron, L.
! Further cleaning: Thirel, G.
!------------------------------------------------------------------------------
! Creation date: 2003
! Last modified: 25/11/2019
!------------------------------------------------------------------------------
! REFERENCES
! Mathevet, T. (2005). Quels modèles pluie-débit globaux pour le pas de temps
! horaire ? Développement empirique et comparaison de modèles sur un large
! échantillon de bassins versants. PhD thesis (in French), ENGREF - Cemagref
! Antony, Paris, France.
!
! 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 (in French), UPMC - Cemagref Antony, Paris, France.
!------------------------------------------------------------------------------
! 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])
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
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.
PS=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)=PN ! Pn ! [numeric] net rainfall [mm/hour]
MISC( 5)=PS ! Ps ! [numeric] part of Ps filling the production store [mm/hour]
MISC( 6)=AE ! AE ! [numeric] actual evapotranspiration [mm/hour]
MISC( 7)=PERC ! Perc ! [numeric] percolation (PERC) [mm/hour]
MISC( 8)=PR ! PR ! [numeric] PR=PN-PS+PERC [mm/hour]
MISC( 9)=StUH1(1) ! Q9 ! [numeric] outflow from UH1 (Q9) [mm/hour]
MISC(10)=StUH2(1) ! Q1 ! [numeric] outflow from UH2 (Q1) [mm/hour]
MISC(11)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
MISC(12)=EXCH ! Exch ! [numeric] potential semi-exchange between catchments (EXCH) [mm/hour]
MISC(13)=AEXCH1 ! AExch1 ! [numeric] actual exchange between catchments from routing store (AEXCH1) [mm/hour]
MISC(14)=AEXCH2 ! AExch2 ! [numeric] actual exchange between catchments from direct branch (after UH2) (AEXCH2) [mm/hour]
MISC(15)=AEXCH1+AEXCH2 ! AExch ! [numeric] actual total exchange between catchments (AEXCH1+AEXCH2) [mm/hour]
MISC(16)=QR ! QR ! [numeric] outflow from routing store (QR) [mm/hour]
MISC(17)=QD ! QD ! [numeric] outflow from UH2 branch after exchange (QD) [mm/hour]
MISC(18)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/hour]
END SUBROUTINE
src/frun_GR4J.f deleted 100644 → 0
View file @ 0967defb
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
c StateEnd = -999.999 !initialization made in R
c Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
!Time loop
!--------------------------------------------------------------
DO k=1,LInputs
P1=InputsPrecip(k)
E =InputsPE(k)
c Q = -999.999
c 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
c################################################################################################################################
C**********************************************************************
SUBROUTINE MOD_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
&MISC)
C Run on a single time step with the GR4J model
C Inputs:
C St Vector of model states in stores at the beginning of the time step [mm]
C StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
C StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
C OrdUH1 Vector of ordinates in UH1 [-]
C OrdUH2 Vector of ordinates in UH2 [-]
C Param Vector of model parameters [various units]
C P1 Value of rainfall during the time step [mm/day]
C E Value of potential evapotranspiration during the time step [mm/day]
C Outputs:
C St Vector of model states in stores at the end of the time step [mm]
C StUH1 Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
C StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
C Q Value of simulated flow at the catchment outlet for the time step [mm/day]
C MISC Vector of model outputs for the time step [mm/day]
C**********************************************************************
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)
C 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))
! fin 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))
! fin speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
C 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))
! fin speed-up
St(1)=St(1)-PERC
PR=PR+PERC
C Split of effective rainfall into the two routing components
PRHU1=PR*B
PRHU2=PR*(1.-B)
C 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
C 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
C 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
! fin speed-up
C 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.))
! fin speed-up
St(2)=St(2)-QR
C Runoff from direct branch QD
AEXCH2=EXCH
IF((StUH2(1)+EXCH).LT.0.) AEXCH2=-StUH2(1)
QD=MAX(0.d0,StUH2(1)+EXCH)
C Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
C 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
src/frun_GR4J.f90 0 → 100644
View file @ 4bccfee2
!------------------------------------------------------------------------------
! Subroutines relative to the annual GR4J model
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_GR4J.f
!------------------------------------------------------------------------------
! AUTHORS
! Original code: Perrin, C.
! Cleaning and formatting for airGR: Coron, L.
! Further cleaning: Thirel, G.
!------------------------------------------------------------------------------
! Creation date: 2000
! Last modified: 25/11/2019
!------------------------------------------------------------------------------
! REFERENCES
! Perrin, C., Michel, C. and Andréassian, V. (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])
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
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 ! [numeric] observed potential evapotranspiration [mm/day]
MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/day]
MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
MISC( 4)=PN ! Pn ! [numeric] net rainfall [mm/day]
MISC( 5)=PS ! Ps ! [numeric] part of Ps filling the production store [mm/day]
MISC( 6)=AE ! AE ! [numeric] actual evapotranspiration [mm/day]
MISC( 7)=PERC ! Perc ! [numeric] percolation (PERC) [mm/day]
MISC( 8)=PR ! PR ! [numeric] PR=PN-PS+PERC [mm/day]
MISC( 9)=StUH1(1) ! Q9 ! [numeric] outflow from UH1 (Q9) [mm/day]
MISC(10)=StUH2(1) ! Q1 ! [numeric] outflow from UH2 (Q1) [mm/day]
MISC(11)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
MISC(12)=EXCH ! Exch ! [numeric] potential semi-exchange between catchments (EXCH) [mm/day]
MISC(13)=AEXCH1 ! AExch1 ! [numeric] actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
MISC(14)=AEXCH2 ! AExch2 ! [numeric] actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
MISC(15)=AEXCH1+AEXCH2 ! AExch ! [numeric] actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
MISC(16)=QR ! QR ! [numeric] outflow from routing store (QR) [mm/day]
MISC(17)=QD ! QD ! [numeric] outflow from UH2 branch after exchange (QD) [mm/day]
MISC(18)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/day]
END SUBROUTINE
src/frun_GR5H.f90 0 → 100644
View file @ 4bccfee2
!------------------------------------------------------------------------------
! Subroutines relative to the annual GR5H model
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_GR5H.f
!------------------------------------------------------------------------------
! AUTHORS
! Original code: Le Moine, N., Ficchì, A.
! Cleaning and formatting for airGR: Coron, L.
! Further cleaning: Thirel, G.
!------------------------------------------------------------------------------
! Creation date: 2006
! Last modified: 26/11/2019
!------------------------------------------------------------------------------
! REFERENCES
! Ficchi, A. (2017). An adaptive hydrological model for multiple time-steps:
! Diagnostics and improvements based on fluxes consistency. PhD thesis,
! UPMC - Irstea Antony, Paris, France.
!
! Ficchi, A., Perrin, C. and Andréassian, V. (2019). Hydrological modelling at
! multiple sub-daily time steps: model improvement via flux-matching. Journal
! of Hydrology, 575, 1308-1327, doi: 10.1016/j.jhydrol.2019.05.084.
!------------------------------------------------------------------------------
! Quick description of public procedures:
! 1. frun_gr5h
! 2. MOD_GR5H
!------------------------------------------------------------------------------
SUBROUTINE frun_gr5h(LInputs,InputsPrecip,InputsPE,NParam,Param, &
NStates,StateStart,Imax,NOutputs,IndOutputs, &
Outputs,StateEnd)
! Subroutine that initializes GR5H, get its parameters, performs the call
! to the MOD_GR5H subroutine at each time step, and stores the final states
! Inputs
! LInputs ! Integer, length of input and output series
! InputsPrecip ! Vector of real, input series of total precipitation [mm/hour]
! InputsPE ! Vector of real, input series of potential evapotranspiration (PE) [mm/hour]
! NParam ! Integer, number of model parameters
! Param ! Vector of real, parameter set
! NStates ! Integer, number of state variables
! StateStart ! Vector of real, state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
! Imax ! Real, fixed capacity of the interception store [mm] (used only if IsIntStore >= 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 Unit Hydrograph (UH) storages [mm])
Implicit None
!! dummies
! in
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, dimension(LInputs), intent(in) :: InputsPrecip
doubleprecision, dimension(LInputs), intent(in) :: InputsPE
doubleprecision, dimension(NParam), intent(in) :: Param
doubleprecision, dimension(NStates), intent(in) :: StateStart
doubleprecision, intent(in) :: Imax ! interception capacity (fixed parameter) used only if IsIntStore >= 0
integer, dimension(NOutputs), intent(in) :: IndOutputs
! out
doubleprecision, dimension(NStates), intent(out) :: StateEnd
doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
!! locals
logical :: IsIntStore ! TRUE if interception store is used, FALSE otherwise
integer :: I,K
integer, parameter :: NH=480,NMISC=30
doubleprecision, dimension(3) :: St
doubleprecision, dimension(2*NH) :: StUH2, OrdUH2
doubleprecision, dimension(NMISC) :: MISC
doubleprecision :: D,P1,E,Q
IF (Imax .LT. 0.d0) THEN
IsIntStore = .FALSE.
ELSE
IsIntStore = .TRUE.
ENDIF
!--------------------------------------------------------------
! Initializations
!--------------------------------------------------------------
! initialization of model states to zero
St=0.
StUH2=0.
! initialization of model states using StateStart
St(1) = StateStart(1)
St(2) = StateStart(2)
IF (IsIntStore .EQV. .TRUE.) St(3) = StateStart(4)
DO I=1,2*NH
StUH2(I)=StateStart(7+NH+I)
ENDDO
! parameter values
! Param(1) : production store capacity (X1 - PROD) [mm]
! Param(2) : intercatchment exchange coefficient (X2 - CES1) [mm/h]
! Param(3) : routing store capacity (X3 - ROUT) [mm]
! Param(4) : time constant of unit hydrograph (X4 - TB) [hour]
! Param(5) : intercatchment exchange threshold (X5 - CES2) [-]
!computation of HU ordinates
OrdUH2 = 0.
D=1.25
CALL UH2_H(OrdUH2,Param(4),D)
! initialization of model outputs
Q = -999.999
MISC = -999.999
! StateEnd = -999.999 !initialization made in R
! Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
! Time loop
!--------------------------------------------------------------
DO k=1,LInputs
P1=InputsPrecip(k)
E =InputsPE(k)
! Q = -999.999
! MISC = -999.999
! model run on one time step
CALL MOD_GR5H(St,StUH2,OrdUH2,Param,IsIntStore,Imax,P1,E,Q,MISC)
! storage of outputs
DO I=1,NOutputs
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
! model states at the end of the run
StateEnd(1) = St(1)
StateEnd(2) = St(2)
StateEnd(4) = St(3)
DO K=1,2*NH
StateEnd(7+NH+K)=StUH2(K)
ENDDO
RETURN
ENDSUBROUTINE
!################################################################################################################################
!**********************************************************************
SUBROUTINE MOD_GR5H(St,StUH2,OrdUH2,Param,IsIntStore,Imax,P1,E,Q,MISC)
! Calculation of streamflow on a single time step (hour) with the GR5H model
! Inputs:
! St Vector of real, model states in stores at the beginning of the time step [mm]
! StUH2 Vector of real, model states in Unit Hydrograph 2 at the beginning of the time step [mm/h]
! OrdUH2 Vector of real, ordinates in UH2 [-]
! Param Vector of real, model parameters [various units]
! IsIntStore Logical, whether interception store is used
! Imax Real, value of interception store capacity, fixed parameter [mm]
! P1 Real, value of rainfall during the time step [mm/h]
! E Real, value of potential evapotranspiration during the time step [mm/h]
! Outputs:
! St Vector of real, model states in stores at the end of the time step [mm/h]
! StUH2 Vector of real, model states in Unit Hydrograph 2 at the end of the time step [mm/h]
! Q Real, value of simulated flow at the catchment outlet for the time step [mm/h]
! MISC Vector of real, model outputs for the time step [mm or mm/h]
!**********************************************************************
Implicit None
!! locals
integer, parameter :: NParam=5,NMISC=30,NH=480
doubleprecision :: A,EN,ES,PN,PR,PS,WS,EI
doubleprecision :: PERC,EXCH,QR,QD,Q1,Q9
doubleprecision :: AE,AEXCH1,AEXCH2
integer :: K
doubleprecision, parameter :: B=0.9
doubleprecision, parameter :: stored_val=759.69140625
doubleprecision :: expWS, TWS, Sr, Rr ! speed-up
!! dummies
! in
doubleprecision, dimension(NParam), intent(in) :: Param
logical, intent(in) :: IsIntStore
doubleprecision, intent(in) :: P1,E,Imax
doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
! inout
doubleprecision, dimension(3), intent(inout) :: St
doubleprecision, dimension(2*NH), intent(inout) :: StUH2
! out
doubleprecision, intent(out) :: Q
doubleprecision, dimension(NMISC), intent(out) :: MISC
A=Param(1)
! Interception and production store
IF (IsIntStore .EQV. .TRUE.) THEN
! MODIFIED - A. Ficchi
! Calculation of interception fluxes [EI] and throughfall [PTH]
! & update of the Interception store state, St(3)
! Interception store calculation, with evaporation prior to throughfall
EI=MIN(E, P1+St(3))
PN=MAX(0.d0, P1-(Imax-St(3))-EI)
St(3)=St(3)+P1-EI-PN
EN=MAX(0.d0, E-EI)
! Production (SMA) store, saving the total actual evaporation including evaporation from interception store (EI)
IF(EN.GT.0) THEN
WS=EN/A
IF(WS.GT.13)WS=13.
expWS=exp(2.*WS)
TWS = (expWS - 1.)/(expWS + 1.)
Sr = St(1)/A
ES=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
St(1)=St(1)-ES
AE=ES+EI
ELSE
AE=EI
ES = 0.
ENDIF
IF(PN.GT.0.) THEN
WS=PN/A
IF(WS.GT.13) WS=13.
expWS=exp(2.*WS)
TWS = (expWS - 1.)/(expWS + 1.)
Sr = St(1)/A
PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
PR=PN-PS
St(1)=St(1)+PS
ELSE
PS=0.
PR=0.
ENDIF
ENDIF
IF (IsIntStore .EQV. .FALSE.) THEN
IF(P1.LE.E) THEN
EN=E-P1
PN=0.
WS=EN/A
IF(WS.GT.13.) WS=13.
! speed-up
expWS=exp(2.*WS)
TWS = (expWS - 1.)/(expWS + 1.)
Sr = St(1)/A
ES=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
! ES=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
! end speed-up
AE=ES+P1
EI = P1
St(1)=St(1)-ES
PS=0.
PR=0.
ELSE
EN=0.
ES=0.
AE=E
EI = 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
ENDIF
! Percolation from production store
IF(St(1).LT.0.) St(1)=0.
! speed-up
! (21/4)**4 = 759.69140625 = stored_val
Sr = St(1)/Param(1)
Sr = Sr * Sr
Sr = Sr * Sr
PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/stored_val)))
! PERC=X(2)*(1.-(1.+(X(2)/(21./4.*Param(1)))**4.)**(-0.25))
! end speed-up
St(1)=St(1)-PERC
PR=PR+PERC
! Convolution of unit hydrograph UH2
DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1.)))
StUH2(K)=StUH2(K+1)+OrdUH2(K)*PR
ENDDO
StUH2(2*NH)=OrdUH2(2*NH)*PR
! Split of unit hydrograph first component into the two routing components
Q9=StUH2(1)*B
Q1=StUH2(1)*(1.-B)
! Potential intercatchment semi-exchange
EXCH=Param(2)*(St(2)/Param(3)-Param(5))
! Routing store
AEXCH1=EXCH
IF((St(2)+Q9+EXCH).LT.0.) AEXCH1=-St(2)-Q9
St(2)=St(2)+Q9+EXCH
IF(St(2).LT.0.) St(2)=0.
! speed-up
Rr = St(2)/Param(3)
Rr = Rr * Rr
Rr = Rr * Rr
QR = St(2)*(1.-1./SQRT(SQRT(1.+Rr)))
! QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
! end speed-up
St(2)=St(2)-QR
! Runoff from direct branch QD
AEXCH2=EXCH
IF((Q1+EXCH).LT.0.) AEXCH2=-Q1
QD=MAX(0.d0,Q1+EXCH)
! Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
! Variables storage
MISC( 1)=E ! PE ! [numeric] observed potential evapotranspiration [mm/h]
MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/h]
MISC( 3)=St(3) ! Interc ! [numeric] interception store level (St(3)) [mm]
MISC( 4)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
MISC( 5)=PN ! Pn ! [numeric] net rainfall [mm/h]
MISC( 6)=PS ! Ps ! [numeric] part of Ps filling the production store [mm/h]
MISC( 7)=AE ! AE ! [numeric] actual evapotranspiration [mm/h]
MISC( 8)=EI ! EI ! [numeric] evapotranspiration from rainfall neutralisation or interception store [mm/h]
MISC( 9)=ES ! ES ! [numeric] evapotranspiration from production store [mm/h]
MISC(10)=PERC ! Perc ! [numeric] percolation (PERC) [mm/h]
MISC(11)=PR ! PR ! [numeric] PR=PN-PS+PERC [mm/h]
MISC(12)=Q9 ! Q9 ! [numeric] outflow from UH1 (Q9) [mm/h]
MISC(13)=Q1 ! Q1 ! [numeric] outflow from UH2 (Q1) [mm/h]
MISC(14)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
MISC(15)=EXCH ! Exch ! [numeric] potential semi-exchange between catchments (EXCH) [mm/h]
MISC(16)=AEXCH1 ! AExch1 ! [numeric] actual exchange between catchments from branch 1 (AEXCH1) [mm/h]
MISC(17)=AEXCH2 ! AExch2 ! [numeric] actual exchange between catchments from branch 2 (AEXCH2) [mm/h]
MISC(18)=AEXCH1+AEXCH2 ! AExch ! [numeric] actual total exchange between catchments (AEXCH1+AEXCH2) [mm/h]
MISC(19)=QR ! QR ! [numeric] outflow from routing store (QR) [mm/h]
MISC(20)=QD ! QD ! [numeric] outflow from UH2 branch after exchange (QD) [mm/h]
MISC(21)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/h]
ENDSUBROUTINE
src/frun_GR5J.f deleted 100644 → 0
View file @ 0967defb
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
!--------------------------------------------------------------
!initilisation of model states to zero
St=0.
StUH2=0.
!initilisation 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
c StateEnd = -999.999 !initialization made in R
c Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
!Time loop
!--------------------------------------------------------------
DO k=1,LInputs
P1=InputsPrecip(k)
E =InputsPE(k)
c Q = -999.999
c 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
c################################################################################################################################
C**********************************************************************
SUBROUTINE MOD_GR5J(St,StUH2,OrdUH2,Param,P1,E,Q,
&MISC)
C Run on a single time step with the GR5J model
C Inputs:
C St Vector of model states in stores at the beginning of the time step [mm]
C StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
C OrdUH2 Vector of ordinates in UH2 [-]
C Param Vector of model parameters [various units]
C P1 Value of rainfall during the time step [mm]
C E Value of potential evapotranspiration during the time step [mm]
C Outputs:
C St Vector of model states in stores at the end of the time step [mm]
C StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
C Q Value of simulated flow at the catchment outlet for the time step [mm/day]
C MISC Vector of model outputs for the time step [mm or mm/day]
C**********************************************************************
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)
C 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))
! fin 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))
! fin speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
C 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))
! fin speed-up
St(1)=St(1)-PERC
PR=PR+PERC
C 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
C Split of unit hydrograph first component into the two routing components
Q9=StUH2(1)*B
Q1=StUH2(1)*(1.-B)
C Potential intercatchment semi-exchange
EXCH=Param(2)*(St(2)/Param(3)-Param(5))
C 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.))
! fin speed-up
St(2)=St(2)-QR
C Runoff from direct branch QD
AEXCH2=EXCH
IF((Q1+EXCH).LT.0.) AEXCH2=-Q1
QD=MAX(0.d0,Q1+EXCH)
C Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
C 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
src/frun_GR5J.f90 0 → 100644
View file @ 4bccfee2
!------------------------------------------------------------------------------
! Subroutines relative to the annual GR5J model
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_GR5J.f
!------------------------------------------------------------------------------
! AUTHORS
! Original code: Le Moine, N.
! Cleaning and formatting for airGR: Coron, L.
! Further cleaning: Thirel, G.
!------------------------------------------------------------------------------
! 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 (in French), UPMC - Cemagref Antony, Paris, France.
!
! Pushpalatha, R., Perrin, C., Le Moine, N., Mathevet, T. and Andréassian, V.
! (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])
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
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 ! [numeric] observed potential evapotranspiration [mm/day]
MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/day]
MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
MISC( 4)=PN ! Pn ! [numeric] net rainfall [mm/day]
MISC( 5)=PS ! Ps ! [numeric] part of Ps filling the production store [mm/day]
MISC( 6)=AE ! AE ! [numeric] actual evapotranspiration [mm/day]
MISC( 7)=PERC ! Perc ! [numeric] percolation (PERC) [mm/day]
MISC( 8)=PR ! PR ! [numeric] PR=PN-PS+PERC [mm/day]
MISC( 9)=Q9 ! Q9 ! [numeric] outflow from UH1 (Q9) [mm/day]
MISC(10)=Q1 ! Q1 ! [numeric] outflow from UH2 (Q1) [mm/day]
MISC(11)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
MISC(12)=EXCH ! Exch ! [numeric] potential semi-exchange between catchments (EXCH) [mm/day]
MISC(13)=AEXCH1 ! AExch1 ! [numeric] actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
MISC(14)=AEXCH2 ! AExch2 ! [numeric] actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
MISC(15)=AEXCH1+AEXCH2 ! AExch ! [numeric] actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
MISC(16)=QR ! QR ! [numeric] outflow from routing store (QR) [mm/day]
MISC(17)=QD ! QD ! [numeric] outflow from UH2 branch after exchange (QD) [mm/day]
MISC(18)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/day]
END SUBROUTINE
src/frun_GR6J.f deleted 100644 → 0
View file @ 0967defb
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])
!DEC$ ATTRIBUTES DLLEXPORT :: frun_GR6J
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(3), StUH1(NH), StUH2(2*NH)
doubleprecision OrdUH1(NH), OrdUH2(2*NH)
doubleprecision MISC(NMISC)
doubleprecision D
doubleprecision P1,E,Q
integer I,K
!--------------------------------------------------------------
!Initializations
!--------------------------------------------------------------
!initilisation 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)
St(3) = StateStart(3)
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 - CES1) [mm/day]
!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) [-]
!Param(6) : time constant of exponential store (X6 - EXP) [day]
!computation of HU 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
c StateEnd = -999.999 !initialization made in R
c Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
!Time loop
!--------------------------------------------------------------
DO k=1,LInputs
P1=InputsPrecip(k)
E =InputsPE(k)
c Q = -999.999
c MISC = -999.999
!model run on one time step
CALL MOD_GR6J(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)
StateEnd(3) = St(3)
DO K=1,NH
StateEnd(7+K)=StUH1(K)
ENDDO
DO K=1,2*NH
StateEnd(7+NH+K)=StUH2(K)
ENDDO
RETURN
ENDSUBROUTINE
c################################################################################################################################
C**********************************************************************
SUBROUTINE MOD_GR6J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
&MISC)
C Run on a single time step with the GR6J model
C Inputs:
C St Vector of model states in stores at the beginning of the time step [mm]
C StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
C StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
C OrdUH1 Vector of ordinates in UH1 [-]
C OrdUH2 Vector of ordinates in UH2 [-]
C Param Vector of model parameters [various units]
C P1 Value of rainfall during the time step [mm]
C E Value of potential evapotranspiration during the time step [mm]
C Outputs:
C St Vector of model states in stores at the end of the time step [mm]
C StUH1 Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
C StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
C Q Value of simulated flow at the catchment outlet for the time step [mm/day]
C MISC Vector of model outputs for the time step [mm or mm/day]
C**********************************************************************
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
A=Param(1)
C 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))
! fin 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))
! fin speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
C 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))
! fin speed-up
St(1)=St(1)-PERC
PR=PR+PERC
C Split of effective rainfall into the two routing components
PRUH1=PR*B
PRUH2=PR*(1.-B)
C 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
ENDDO
StUH1(NH)=OrdUH1(NH)*PRUH1
C 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
ENDDO
StUH2(2*NH)=OrdUH2(2*NH)*PRUH2
C Potential intercatchment semi-exchange
EXCH=Param(2)*(St(2)/Param(3)-Param(5))
C Routing store
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.
! 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.))
! fin speed-up
St(2)=St(2)-QR
C 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.7.)THEN
QRExp=St(3)+Param(6)/EXP(AR)
GOTO 3
ENDIF
IF(AR.LT.-7.)THEN
QRExp=Param(6)*EXP(AR)
GOTO 3
ENDIF
QRExp=Param(6)*LOG(EXP(AR)+1.)
3 CONTINUE
St(3)=St(3)-QRExp
C Runoff from direct branch QD
AEXCH2=EXCH
IF((StUH2(1)+EXCH).LT.0) AEXCH2=-StUH2(1)
QD=MAX(0.d0,StUH2(1)+EXCH)
C Total runoff
Q=QR+QD+QRExp
IF(Q.LT.0.) Q=0.
C 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 third-exchange between catchments (EXCH) [mm/day]
MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from routing store (AEXCH1) [mm/day]
MISC(14)=AEXCH2 ! AExch2 ! actual exchange between catchments from direct branch (after UH2) (AEXCH2) [mm/day]
MISC(15)=AEXCH1+AEXCH2+EXCH ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2+EXCH) [mm/day]
MISC(16)=QR ! QR ! outflow from routing store (QR) [mm/day]
MISC(17)=QRExp ! QRExp ! outflow from exponential store (QRExp) [mm/day]
MISC(18)=St(3) ! Exp ! exponential store level (St(3)) (negative) [mm]
MISC(19)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/day]
MISC(20)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
ENDSUBROUTINE
src/frun_GR6J.f90 0 → 100644
View file @ 4bccfee2
!------------------------------------------------------------------------------
! Subroutines relative to the annual GR6J model
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_GR6J.f
!------------------------------------------------------------------------------
! AUTHORS
! Original code: Pushpalatha, R.
! Cleaning and formatting for airGR: Coron, L.
! Further cleaning: Thirel, G.
!------------------------------------------------------------------------------
! Creation date: 2010
! Last modified: 25/11/2019
!------------------------------------------------------------------------------
! REFERENCES
! Pushpalatha, R., Perrin, C., Le Moine, N., Mathevet, T. and Andréassian, V.
! (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])
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(3) :: 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)
St(3) = StateStart(3)
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 - CES1) [mm/day]
! 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) [-]
! Param(6) : time constant of exponential store (X6 - EXP) [day]
! computation of HU 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_GR6J(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)
StateEnd(3) = St(3)
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_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 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 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
!! locals
integer, parameter :: NParam=6,NMISC=30,NH=20
doubleprecision :: A,EN,ER,PN,PR,PS,WS,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)
! 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
PRUH1=PR*B
PRUH2=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)*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
ENDDO
StUH2(2*NH)=OrdUH2(2*NH)*PRUH2
! Potential intercatchment semi-exchange
EXCH=Param(2)*(St(2)/Param(3)-Param(5))
! Routing store
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.
! 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
! 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.7.) THEN
QRExp=St(3)+Param(6)/EXP(AR)
ELSEIF(AR.LT.-7.) THEN
QRExp=Param(6)*EXP(AR)
ELSE
QRExp=Param(6)*LOG(EXP(AR)+1.)
ENDIF
St(3)=St(3)-QRExp
! 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+QRExp
IF(Q.LT.0.) Q=0.
! Variables storage
MISC( 1)=E ! PE ! [numeric] observed potential evapotranspiration [mm/day]
MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/day]
MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
MISC( 4)=PN ! Pn ! [numeric] net rainfall [mm/day]
MISC( 5)=PS ! Ps ! [numeric] part of Ps filling the production store [mm/day]
MISC( 6)=AE ! AE ! [numeric] actual evapotranspiration [mm/day]
MISC( 7)=PERC ! Perc ! [numeric] percolation (PERC) [mm/day]
MISC( 8)=PR ! PR ! [numeric] PR=PN-PS+PERC [mm/day]
MISC( 9)=StUH1(1) ! Q9 ! [numeric] outflow from UH1 (Q9) [mm/day]
MISC(10)=StUH2(1) ! Q1 ! [numeric] outflow from UH2 (Q1) [mm/day]
MISC(11)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
MISC(12)=EXCH ! Exch ! [numeric] potential third-exchange between catchments (EXCH) [mm/day]
MISC(13)=AEXCH1 ! AExch1 ! [numeric] actual exchange between catchments from routing store (AEXCH1) [mm/day]
MISC(14)=AEXCH2 ! AExch2 ! [numeric] actual exchange between catchments from direct branch (after UH2) (AEXCH2) [mm/day]
MISC(15)=AEXCH1+AEXCH2+EXCH ! AExch ! [numeric] actual total exchange between catchments (AEXCH1+AEXCH2+EXCH) [mm/day]
MISC(16)=QR ! QR ! [numeric] outflow from routing store (QR) [mm/day]
MISC(17)=QRExp ! QRExp ! [numeric] outflow from exponential store (QRExp) [mm/day]
MISC(18)=St(3) ! Exp ! [numeric] exponential store level (St(3)) (negative) [mm]
MISC(19)=QD ! QD ! [numeric] outflow from UH2 branch after exchange (QD) [mm/day]
MISC(20)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/day]
END SUBROUTINE
src/frun_PE.f90 0 → 100644
View file @ 4bccfee2
!------------------------------------------------------------------------------
! Subroutines relative to the Oudin potential evapotranspiration (PE) formula
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_PE.f90
!------------------------------------------------------------------------------
! AUTHORS
! Original code: Oudin, L.
! Cleaning and formatting for airGR: Bourgin, F.
! Further cleaning: Delaigue, O., Thirel, G.
!------------------------------------------------------------------------------
! Creation date: 2004
! Last modified: 20/10/2020
!------------------------------------------------------------------------------
! REFERENCES
! Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V.,
! Anctil, F. and Loumagne, C. (2005). Which potential evapotranspiration
! input for a lumped rainfall-runoff model? Part 2 - Towards a simple and
! efficient potential evapotranspiration model for rainfall-runoff modelling.
! Journal of Hydrology, 303(1-4), 290-306, doi: 10.1016/j.jhydrol.2004.08.026.
!------------------------------------------------------------------------------
! Quick description of public procedures:
! 1. frun_pe_oudin
! 2. PE_OUDIN
!------------------------------------------------------------------------------
!*******************************************************************************
SUBROUTINE frun_pe_oudin(LInputs,InputsLAT,InputsTemp,InputsJJ,OutputsPE)
!*******************************************************************************
! Subroutine that performs the call to the PE_OUDIN subroutine at each time step,
! and stores the final values
! Inputs
! LInputs ! Integer, length of input and output series
! InputsLAT ! Vector of real, input series of latitude [rad]
! InputsTemp ! Vector of real, input series of air mean temperature [degC]
! InputsJJ ! Vector of real, input series of Julian day [-]
! Outputs
! OutputsPE ! Vector of real, output series of potential evapotranspiration (PE) [mm/time step]
Implicit None
!! dummies
! in
integer, intent(in) :: LInputs
doubleprecision, dimension(LInputs), intent(in) :: InputsLAT
doubleprecision, dimension(LInputs), intent(in) :: InputsTemp
doubleprecision, dimension(LInputs), intent(in) :: InputsJJ
! out
doubleprecision, dimension(LInputs), intent(out) :: OutputsPE
!! locals
integer :: k
doubleprecision :: FI, tt, jj, PEoud
!--------------------------------------------------------------
! Time loop
!--------------------------------------------------------------
DO k = 1, LInputs
tt = InputsTemp(k)
jj = InputsJJ(k)
FI = InputsLAT(k)!
!model run on one time step
CALL PE_OUDIN(FI, tt, jj, PEoud)
!storage of outputs
OutputsPE(k) = PEoud
ENDDO
RETURN
ENDSUBROUTINE
!################################################################################################################################
!*******************************************************************************
SUBROUTINE PE_OUDIN(FI,DT,JD,DPE)
!*******************************************************************************
! Calculation of potential evapotranspiration (DPE) on a single time step
! using air temperature and daily extra-atmospheric global radiation
! (that depends only on Julian day)
!
! The PE formula is described in:
! Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V.,
! Anctil, F. and Loumagne, C., 2005. Which potential evapotranspiration
! input for a rainfall-runoff model? Part 2 - Towards a simple and
! efficient PE model for rainfall-runoff modelling. Journal of Hydrology
! 303(1-4), 290-306.
!
! For the calculation of extra-atmospheric global radiation, see Appendix C of
! the article by Morton, F.I., 1983. Operational estimates of areal
! evapotranspiration and their significance to the science and practice
! of hydrology. Journal of Hydrology 66 (1/4), 1-76.
!
!***************************************************************
! Inputs
! FI ! Latitude [rad]
! DT ! Air Temperature [degC]
! JD ! Julian day [-]
!
! Outputs
! DPE ! Potential evapotranspiration [mm/time step]
!***************************************************************
IMPLICIT NONE
!! dummies
! in
doubleprecision, intent(in) :: FI, DT, JD
! out
doubleprecision, intent(out) :: DPE
!! locals
doubleprecision :: COSFI, TETA, COSTETA, COSGZ, GZ, COSGZ2
doubleprecision :: SINGZ, COSOM, COSOM2, SINOM, COSPZ, OM, GE
doubleprecision :: ETA
! Calculation of extra-atmospheric global radiation (Appendix C in Morton
! (1983), Eq. C-6 to C-11, p.60-61)
COSFI=COS(FI)
! TETA: Declination of the sun in radians
TETA=0.4093*SIN(JD/58.1-1.405)
COSTETA=COS(TETA)
COSGZ=MAX(0.001d0,COS(FI-TETA))
! GZ: Noon angular zenith distance of the sun
GZ=ACOS(COSGZ)
COSGZ2=COSGZ*COSGZ
IF(COSGZ2.GE.1.) THEN
SINGZ=0.
ELSE
SINGZ=SQRT(1.-COSGZ2)
ENDIF
COSOM=1.-COSGZ/COSFI/COSTETA
IF(COSOM.LT.-1.) COSOM=-1.
IF(COSOM.GT.1.) COSOM=1.
COSOM2=COSOM*COSOM
IF(COSOM2.GE.1.) THEN
SINOM=0.
ELSE
SINOM=SQRT(1.-COSOM2)
ENDIF
OM=ACOS(COSOM)
! PZ: Average angular zenith distance of the sun
COSPZ=COSGZ+COSFI*COSTETA*(SINOM/OM-1.)
IF(COSPZ.LT.0.001) COSPZ=0.001
! ETA: Radius vector of the sun
ETA=1.+COS(JD/58.1)/30.
! GE: extra-atmospheric global radiation
GE=446.*OM*COSPZ*ETA
! Daily PE by Oudin et al. (2006) formula:
DPE=MAX(0.d0,GE*(DT+5.)/100./28.5)
RETURN
END SUBROUTINE PE_OUDIN
!*******************************************************************************