From 66712dbf83b5c267f86064c1480c5835685783a6 Mon Sep 17 00:00:00 2001
From: Olivier Delaigue <olivier.delaigue@irstea.fr>
Date: Thu, 14 Apr 2016 16:03:53 +0200
Subject: [PATCH] Mise au propre des noms de variables
---
src/frun_GR1A.f | 63 +++++++------
src/frun_GR2M.f | 121 ++++++++++++-------------
src/frun_GR4H.f | 232 ++++++++++++++++++++++++++++++------------------
src/frun_GR4J.f | 230 ++++++++++++++++++++++++-----------------------
src/frun_GR5J.f | 184 +++++++++++++++++++-------------------
src/frun_GR6J.f | 208 ++++++++++++++++++++++++-------------------
6 files changed, 562 insertions(+), 476 deletions(-)
diff --git a/src/frun_GR1A.f b/src/frun_GR1A.f
index f9d34dae..62668b20 100644
--- a/src/frun_GR1A.f
+++ b/src/frun_GR1A.f
@@ -2,21 +2,21 @@
SUBROUTINE frun_GR1A(
!inputs
- & LInputs , ! [numeric] length of input and output series
- & InputsPrecip , ! [numeric] input series of total precipitation [mm]
- & InputsPE , ! [numeric] input series PE [mm]
- & NParam , ! [numeric] number of model parameter
- & Param , ! [numeric] parameter set
- & NStates , ! [numeric] number of state variables used for model initialising
- & StateStart , ! [numeric] state variables used when the model run starts (none here)
- & NOutputs , ! [numeric] number of output series
- & IndOutputs , ! [numeric] indices of output series
+ & 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 , ! [numeric] output series
- & StateEnd ) ! [numeric] state variables at the end of the model run (none here)
+ & Outputs , ! [double] output series
+ & StateEnd ) ! [double] state variables at the end of the model run (none here)
- !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr1a
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_GR1A
Implicit None
@@ -38,16 +38,16 @@
integer I,K
!--------------------------------------------------------------
- !Initialisations
+ !Initializations
!--------------------------------------------------------------
!parameter values
- !Param(1) : fonction parameter [mm]
+ !Param(1) : PE adjustment factor [-]
- !initialisation of model outputs
+ !initialization of model outputs
Q = -999.999
MISC = -999.999
-c Outputs = -999.999 !initialisation made in R
+c Outputs = -999.999 !initialization made in R
@@ -60,7 +60,7 @@ c Outputs = -999.999 !initialisation made in R
E1=InputsPE(k)
c Q = -999.999
c MISC = -999.999
- !model run on one time-step
+ !model run on one time step
CALL MOD_GR1A(Param,P0,P1,E1,Q,MISC)
!storage of outputs
DO I=1,NOutputs
@@ -83,15 +83,15 @@ c###############################################################################
C**********************************************************************
SUBROUTINE MOD_GR1A(Param,P0,P1,E1,Q,MISC)
-C Run on a single time-step with the GR4H model
+C Calculation of streamflow on a single time step with the GR1A model
C Inputs:
-C Param Vector of model parameters [mm]
-C P0 Value of rainfall during the previous time-step [mm]
-C P1 Value of rainfall during the time-step [mm]
-C E1 Value of potential evapotranspiration during the time-step [mm]
+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 time-step [mm]
-C MISC Vector of model outputs for the time-step [mm]
+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]
C**********************************************************************
Implicit None
INTEGER NMISC,NParam
@@ -101,14 +101,19 @@ C**********************************************************************
DOUBLEPRECISION MISC(NMISC)
DOUBLEPRECISION P0,P1,E1,Q
+ DOUBLEPRECISION tt ! speed-up
+
C Runoff
- Q=P1*(1.-1./(1.+((0.7*P1+0.3*P0)/Param(1)/E1)**2.)**0.5)
-
+ ! 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] potential evapotranspiration [mm/year]
- MISC( 2)=P1 ! Precip ! [numeric] total precipitation [mm/year]
- MISC( 3)=Q ! Qsim ! [numeric] outflow at catchment outlet [mm/year]
+ MISC( 1)=E1 ! PE ! [numeric] observed potential evapotranspiration [mm/year]
+ MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/year]
+ MISC( 3)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/year]
ENDSUBROUTINE
diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f
index b70245c0..2a15353c 100644
--- a/src/frun_GR2M.f
+++ b/src/frun_GR2M.f
@@ -2,21 +2,21 @@
SUBROUTINE frun_GR2M(
!inputs
- & LInputs , ! [numeric] length of input and output series
- & InputsPrecip , ! [numeric] input series of total precipitation [mm]
- & InputsPE , ! [numeric] input series PE [mm]
- & NParam , ! [numeric] number of model parameter
- & Param , ! [numeric] parameter set
- & NStates , ! [numeric] number of state variables used for model initialising
- & StateStart , ! [numeric] state variables used when the model run starts (reservoir levels [mm])
- & NOutputs , ! [numeric] number of output series
- & IndOutputs , ! [numeric] indices of output series
+ & 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 , ! [numeric] output series
- & StateEnd ) ! [numeric] state variables at the end of the model run (reservoir levels [mm])
+ & Outputs , ! [double] output series
+ & StateEnd ) ! [double] state variables at the end of the model run (store levels [mm])
- !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr2m
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_GR2M
Implicit None
@@ -33,31 +33,31 @@
!parameters, internal states and variables
integer NMISC
parameter (NMISC=30)
- doubleprecision X(2)
+ doubleprecision St(2)
doubleprecision MISC(NMISC)
doubleprecision P,E,Q
integer I,K
!--------------------------------------------------------------
- !Initialisations
+ !Initializations
!--------------------------------------------------------------
- !initilisation of model states to zero
- X=0.
+ !initialization of model states to zero
+ St=0.
!initilisation of model states using StateStart
- X(1)=StateStart(1)
- X(2)=StateStart(2)
+ St(1)=StateStart(1)
+ St(2)=StateStart(2)
!parameter values
!Param(1) : production store capacity [mm]
- !Param(2) : groundwater exchange coefficient [mm/month]
+ !Param(2) : groundwater exchange coefficient [-]
- !initialisation of model outputs
+ !initialization of model outputs
Q = -999.999
MISC = -999.999
-c StateEnd = -999.999 !initialisation made in R
-c Outputs = -999.999 !initialisation made in R
+c StateEnd = -999.999 !initialization made in R
+c Outputs = -999.999 !initialization made in R
@@ -69,16 +69,16 @@ c Outputs = -999.999 !initialisation made in R
E =InputsPE(k)
c Q = -999.999
c MISC = -999.999
- !model run on one time-step
- CALL MOD_GR2M(X,Param,P,E,Q,MISC)
+ !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)=X(1)
- StateEnd(2)=X(2)
+ StateEnd(1)=St(1)
+ StateEnd(2)=St(2)
RETURN
@@ -94,28 +94,28 @@ c###############################################################################
C**********************************************************************
- SUBROUTINE MOD_GR2M(X,Param,P,E,Q,MISC)
-C Run on a single time-step with the GR4H model
+ 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 X Vector of model states at the beginning of the time-step [mm]
-C Param Vector of model parameters [mixed units]
-C P Value of rainfall during the time-step [mm]
-C E Value of potential evapotranspiration during the time-step [mm]
+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 X 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]
-C MISC Vector of model outputs for the time-step [mm]
+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]
C**********************************************************************
Implicit None
INTEGER NMISC,NParam
PARAMETER (NMISC=30)
PARAMETER (NParam=2)
- DOUBLEPRECISION X(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
+ DOUBLEPRECISION P1,P2,P3,R1,R2,AE,EXCH
DOUBLEPRECISION TWS, Sr, Rr ! speed-up
@@ -125,11 +125,11 @@ C Production store
! speed-up
TWS = tanHyp(WS)
- S1=(X(1)+Param(1)*TWS)/(1.+X(1)/Param(1)*TWS)
+ 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+X(1)-S1
+ P1=P+St(1)-S1
WS=E/Param(1)
IF(WS.GT.13)WS=13
@@ -138,54 +138,43 @@ C Production store
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.
- X(1)=S2/Sr**(1./3.)
+ St(1)=S2/Sr**(1./3.)
! X(1)=S2/(1+(S2/Param(1))**3.)**(1./3.)
! fin speed-up
- P2=S2-X(1)
+ P2=S2-St(1)
P3=P1+P2
C QR calculation (routing store)
- R1=X(2)+P3
+ R1=St(2)+P3
C Water exchange
- R2=Param(2)*R1
+ R2=Param(2)*R1
+ EXCH = R2 - R1
C Total runoff
Q=R2*R2/(R2+60.)
C Updating store level
- X(2)=R2-Q
+ St(2)=R2-Q
C Variables storage
- MISC( 1)=E ! PE ! [numeric] potential evapotranspiration [mm/month]
- MISC( 2)=P1 ! Precip ! [numeric] total precipitation [mm/month]
- MISC( 3)=X(1) ! Prod ! [numeric] production store level (X(2)) [mm]
- MISC( 4)=X(2) ! Rout ! [numeric] routing store level (X(1)) [mm]
- MISC( 5)=Q ! Qsim ! [numeric] outflow at catchment outlet [mm/month]
-
-
-C MISC( 1)=E ! PE ! potential evapotranspiration [mm/d]
-C MISC( 2)=P1 ! Precip ! total precipitation [mm/d]
-C MISC( 3)=X(2) ! Prod ! production store level (X(2)) [mm]
-C MISC( 4)=AE ! AE ! actual evapotranspiration [mm/d]
-C MISC( 5)=PERC ! Perc ! percolation (PERC) [mm]
-C MISC( 6)=PR ! PR ! PR=PN-PS+PERC [mm]
-C MISC( 7)=X(8) ! Q9 ! outflow from HU1 (Q9) [mm/d]
-C MISC( 8)=X(8+NH) ! Q1 ! outflow from HU2 (Q1) [mm/d]
-C MISC( 9)=X(1) ! Rout ! routing store level (X(1)) [mm]
-C MISC(10)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/d]
-C MISC(11)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/d]
-C MISC(12)=QR ! QR ! outflow from routing store (QR) [mm/d]
-C MISC(13)=QD ! QD ! outflow from HU2 branch after exchange (QD) [mm/d]
-C MISC(14)=Q ! Qsim ! outflow at catchment outlet [mm/d]
-
+ MISC( 1)=E ! PE ! [numeric] observed potential evapotranspiration [mm/month]
+ MISC( 2)=P1 ! Precip ! [numeric] observed total precipitation [mm/month]
+ MISC( 3)=AE ! AE ! [numeric] actual evapotranspiration [mm/month]
+ MISC( 4)=P2 ! P2 ! [numeric] percolation (P2) [mm/month]
+ MISC( 5)=P3 ! P3 ! [numeric] P3=P1+P2 [mm/month]
+ MISC( 6)=EXCH ! EXCH ! [numeric] groundwater exchange (EXCH) [mm/month]
+ MISC( 7)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm]
+ MISC( 8)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm]
+ MISC( 9)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/month]
ENDSUBROUTINE
diff --git a/src/frun_GR4H.f b/src/frun_GR4H.f
index 773bdd47..3c473319 100644
--- a/src/frun_GR4H.f
+++ b/src/frun_GR4H.f
@@ -2,21 +2,21 @@
SUBROUTINE frun_GR4H(
!inputs
- & LInputs , ! [numeric] length of input and output series
- & InputsPrecip , ! [numeric] input series of total precipitation [mm]
- & InputsPE , ! [numeric] input series PE [mm]
- & NParam , ! [numeric] number of model parameter
- & Param , ! [numeric] parameter set
- & NStates , ! [numeric] number of state variables used for model initialising
- & StateStart , ! [numeric] state variables used when the model run starts (reservoir levels [mm] and HU)
- & NOutputs , ! [numeric] number of output series
- & IndOutputs , ! [numeric] indices of output series
+ & 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 , ! [numeric] output series
- & StateEnd ) ! [numeric] state variables at the end of the model run (reservoir levels [mm] and HU)
+ & 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
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_GR4H
Implicit None
@@ -31,43 +31,53 @@
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
!parameters, internal states and variables
- integer NPX,NH,NMISC
- parameter (NPX=14,NH=480,NMISC=30)
- doubleprecision X(5*NH+7),XV(3*NPX+5*NH)
+ 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
!--------------------------------------------------------------
- !Initialisations
+ !Initializations
!--------------------------------------------------------------
- !initilisation of model states to zero
- X=0.
- XV=0.
+ !initilization of model states using StateStart
+ St=0.
+ StUH1=0.
+ StUH2=0.
- !initilisation of model states using StateStart
- DO I=1,3*NH
- X(I)=StateStart(I)
+ !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 constant (X2 - CES) [mm/h]
+ !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) [h]
+ !Param(4) : time constant of unit hydrograph (X4 - TB) [hour]
- !computation of HU ordinates
+ !computation of UH ordinates
+ OrdUH1 = 0.
+ OrdUH2 = 0.
+
D=1.25
- CALL HU1_H(XV,Param(4),D)
- CALL HU2_H(XV,Param(4),D)
+ CALL UH1_H(OrdUH1,Param(4),D)
+ CALL UH2_H(OrdUH2,Param(4),D)
- !initialisation of model outputs
+ !initialization of model outputs
Q = -999.999
MISC = -999.999
-c StateEnd = -999.999 !initialisation made in R
-c Outputs = -999.999 !initialisation made in R
+c StateEnd = -999.999 !initialization made in R
+c Outputs = -999.999 !initialization made in R
@@ -79,16 +89,21 @@ c Outputs = -999.999 !initialisation made in R
E =InputsPE(k)
c Q = -999.999
c MISC = -999.999
- !model run on one time-step
- CALL MOD_GR4H(X,XV,Param,P1,E,Q,MISC)
+ !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
- DO K=1,3*NH
- StateEnd(K)=X(K)
+ 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
@@ -105,25 +120,31 @@ c###############################################################################
C**********************************************************************
- SUBROUTINE MOD_GR4H(X,XV,Param,P1,E,Q,MISC)
-C Run on a single time-step with the GR4H model
+ 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 X Vector of model states at the beginning of the time-step [mm]
-C XV Vector of model states at the beginning of the time-step [mm]
-C Param Vector of model parameters [mixed units]
-C P1 Value of rainfall during the time-step [mm]
-C E Value of potential evapotranspiration during the time-step [mm]
+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 X Vector of model states at the end of the time-step [mm]
-C XV 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]
-C MISC Vector of model outputs for the time-step [mm]
+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]
+C MISC Vector of model outputs for the time step [mm]
C**********************************************************************
Implicit None
- INTEGER NPX,NH,NMISC,NParam
- PARAMETER (NPX=14,NH=480,NMISC=30)
+ INTEGER NH,NMISC,NParam
+ PARAMETER (NH=480,NMISC=30)
PARAMETER (NParam=4)
- DOUBLEPRECISION X(5*NH+7),XV(3*NPX+5*NH)
+ DOUBLEPRECISION St(2),StUH1(NH)
+ DOUBLEPRECISION StUH2(2*NH),OrdUH1(NH),OrdUH2(2*NH)
DOUBLEPRECISION Param(NParam)
DOUBLEPRECISION MISC(NMISC)
DOUBLEPRECISION P1,E,Q
@@ -134,19 +155,27 @@ C**********************************************************************
DATA B/0.9/
+ DOUBLEPRECISION TWS, Sr, Rr ! speed-up
+
A=Param(1)
-C Production store
+C Interception and production store
IF(P1.LE.E) THEN
EN=E-P1
PN=0.
WS=EN/A
IF(WS.GT.13)WS=13.
- ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+
+ ! 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
- IF(X(2).LT.ER) AE=X(2)+P1
- X(2)=X(2)-ER
+ St(1)=St(1)-ER
PR=0.
ELSE
EN=0.
@@ -154,68 +183,97 @@ C Production store
PN=P1-E
WS=PN/A
IF(WS.GT.13)WS=13.
- PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
+
+ ! 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
- X(2)=X(2)+PS
+ St(1)=St(1)+PS
ENDIF
C Percolation from production store
- IF(X(2).LT.0.)X(2)=0.
- PERC=X(2)*(1.-(1.+(X(2)/(21./4.*Param(1)))**4.)**(-0.25))
- X(2)=X(2)-PERC
+ 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 Unit hydrograph HU1
- DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1)))
- X(7+K)=X(8+K)+XV(3*NPX+K)*PRHU1
+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
- X(7+NH)=XV(3*NPX+NH)*PRHU1
+ StUH1(NH)=OrdUH1(NH)*PRHU1
-C Unit hydrograph HU2
- DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1)))
- X(7+NH+K)=X(8+NH+K)+XV(3*NPX+NH+K)*PRHU2
+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
- X(7+3*NH)=XV(3*NPX+3*NH)*PRHU2
+ StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
C Potential intercatchment semi-exchange
- EXCH=Param(2)*(X(1)/Param(3))**3.5
+ ! 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((X(1)+X(8)+EXCH).LT.0) AEXCH1=-X(1)-X(8)
- X(1)=X(1)+X(8)+EXCH
- IF(X(1).LT.0.)X(1)=0.
- QR=X(1)*(1.-(1.+(X(1)/Param(3))**4.)**(-1./4.))
- X(1)=X(1)-QR
+ 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(1)=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((X(8+NH)+EXCH).LT.0) AEXCH2=-X(8+NH)
- QD=MAX(0.,X(8+NH)+EXCH)
+ IF((StUH2(1)+EXCH).LT.0) AEXCH2=-StUH2(1)
+ QD=MAX(0.,StUH2(1)+EXCH)
C Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
C Variables storage
- MISC( 1)=E ! PE ! [numeric] potential evapotranspiration [mm/h]
- MISC( 2)=P1 ! Precip ! [numeric] total precipitation [mm/h]
- MISC( 3)=X(2) ! Prod ! [numeric] production store level (X(2)) [mm]
- MISC( 4)=AE ! AE ! [numeric] actual evapotranspiration [mm/h]
- MISC( 5)=PERC ! Perc ! [numeric] percolation (PERC) [mm]
- MISC( 6)=PR ! PR ! [numeric] PR=PN-PS+PERC [mm]
- MISC( 7)=X(8) ! Q9 ! [numeric] outflow from HU1 (Q9) [mm/h]
- MISC( 8)=X(8+NH) ! Q1 ! [numeric] outflow from HU2 (Q1) [mm/h]
- MISC( 9)=X(1) ! Rout ! [numeric] routing store level (X(1)) [mm]
- MISC(10)=EXCH ! Exch ! [numeric] potential semi-exchange between catchments (EXCH) [mm/h]
- MISC(11)=AEXCH1+AEXCH2 ! AExch ! [numeric] actual total exchange between catchments (AEXCH1+AEXCH2) [mm/h]
- MISC(12)=QR ! QR ! [numeric] outflow from routing store (QR) [mm/h]
- MISC(13)=QD ! QD ! [numeric] outflow from HU2 branch after exchange (QD) [mm/h]
- MISC(14)=Q ! Qsim ! [numeric] outflow at catchment outlet [mm/h]
+ 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]
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f
index 8db23f83..29302aeb 100644
--- a/src/frun_GR4J.f
+++ b/src/frun_GR4J.f
@@ -3,20 +3,20 @@
SUBROUTINE frun_GR4J(
!inputs
& LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm]
- & InputsPE , ! [double] input series PE [mm]
- & NParam , ! [integer] number of model parameter
+ & 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 used for model initialising
- & StateStart , ! [double] state variables used when the model run starts (reservoir levels [mm] and HU)
+ & 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 (reservoir levels [mm] and HU)
+ & 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
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_GR4J
Implicit None
@@ -31,43 +31,53 @@
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
!parameters, internal states and variables
- integer NPX,NH,NMISC
- parameter (NPX=14,NH=20,NMISC=30)
- doubleprecision X(5*NH+7),XV(3*NPX+5*NH)
+ 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
!--------------------------------------------------------------
- !Initialisations
+ !Initializations
!--------------------------------------------------------------
- !initilisation of model states to zero
- X=0.
- XV=0.
+ !initialization of model states to zero
+ St=0.
+ StUH1=0.
+ StUH2=0.
- !initilisation of model states using StateStart
- DO I=1,3*NH
- X(I)=StateStart(I)
+ !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 constant (X2 - CES) [mm/d]
+ !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) [d]
+ !Param(4) : time constant of unit hydrograph (X4 - TB) [day]
- !computation of HU ordinates
+ !computation of UH ordinates
+ OrdUH1 = 0.
+ OrdUH2 = 0.
+
D=2.5
- CALL HU1(XV,Param(4),D)
- CALL HU2(XV,Param(4),D)
+ CALL UH1(OrdUH1,Param(4),D)
+ CALL UH2(OrdUH2,Param(4),D)
- !initialisation of model outputs
+ !initialization of model outputs
Q = -999.999
MISC = -999.999
-c StateEnd = -999.999 !initialisation made in R
-c Outputs = -999.999 !initialisation made in R
+c StateEnd = -999.999 !initialization made in R
+c Outputs = -999.999 !initialization made in R
@@ -79,16 +89,21 @@ c Outputs = -999.999 !initialisation made in R
E =InputsPE(k)
c Q = -999.999
c MISC = -999.999
- !model run on one time-step
- CALL MOD_GR4J(X,XV,Param,P1,E,Q,MISC)
+ !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
- DO K=1,3*NH
- StateEnd(K)=X(K)
+ 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
@@ -105,25 +120,31 @@ c###############################################################################
C**********************************************************************
- SUBROUTINE MOD_GR4J(X,XV,Param,P1,E,Q,MISC)
-C Run on a single time-step with the GR4J model
+ 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 X Vector of model states at the beginning of the time-step [mm]
-C XV Vector of model states at the beginning of the time-step [mm]
-C Param Vector of model parameters [mixed units]
-C P1 Value of rainfall during the time-step [mm]
-C E Value of potential evapotranspiration during the time-step [mm]
+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 X Vector of model states at the end of the time-step [mm]
-C XV 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]
-C MISC Vector of model outputs for the time-step [mm]
+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]
+C MISC Vector of model outputs for the time step [mm]
C**********************************************************************
Implicit None
- INTEGER NPX,NH,NMISC,NParam
- PARAMETER (NPX=14,NH=20,NMISC=30)
+ INTEGER NH,NMISC,NParam
+ PARAMETER (NH=20,NMISC=30)
PARAMETER (NParam=4)
- DOUBLEPRECISION X(5*NH+7),XV(3*NPX+5*NH)
+ DOUBLEPRECISION St(2),StUH1(NH),StUH2(2*NH)
+ DOUBLEPRECISION OrdUH1(NH),OrdUH2(2*NH)
DOUBLEPRECISION Param(NParam)
DOUBLEPRECISION MISC(NMISC)
DOUBLEPRECISION P1,E,Q
@@ -133,126 +154,117 @@ C**********************************************************************
INTEGER K
DATA B/0.9/
-
DOUBLEPRECISION TWS, Sr, Rr ! speed-up
A=Param(1)
-C Production store
+C Interception and production store
IF(P1.LE.E) THEN
- EN=E-P1
- PN=0.
- WS=EN/A
- IF(WS.GT.13)WS=13.
-
+ EN=E-P1
+ PN=0.
+ WS=EN/A
+ IF(WS.GT.13)WS=13.
! speed-up
- TWS = tanHyp(WS)
- Sr = X(2)/A
- ER=X(2)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+ 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
- IF(X(2).LT.ER) AE=X(2)+P1
- X(2)=X(2)-ER
- PR=0.
+ 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.
-
+ EN=0.
+ AE=E
+ PN=P1-E
+ WS=PN/A
+ IF(WS.GT.13)WS=13.
! speed-up
- TWS = tanHyp(WS)
- Sr = X(2)/A
- PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
+ 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
- X(2)=X(2)+PS
+ PR=PN-PS
+ St(1)=St(1)+PS
ENDIF
C Percolation from production store
- IF(X(2).LT.0.)X(2)=0.
-
+ IF(St(1).LT.0.)St(1)=0.
! speed-up
! (9/4)**4 = 25.62891
- Sr = X(2)/Param(1)
+ Sr = St(1)/Param(1)
Sr = Sr * Sr
Sr = Sr * Sr
- PERC=X(2)*(1.-1./SQRT(SQRT(1.+Sr/25.62891)))
+ 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
-
- X(2)=X(2)-PERC
+ 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 Unit hydrograph HU1
- DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1)))
- X(7+K)=X(8+K)+XV(3*NPX+K)*PRHU1
+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
- X(7+NH)=XV(3*NPX+NH)*PRHU1
+ StUH1(NH)=OrdUH1(NH)*PRHU1
-C Unit hydrograph HU2
- DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1)))
- X(7+NH+K)=X(8+NH+K)+XV(3*NPX+NH+K)*PRHU2
+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
- X(7+3*NH)=XV(3*NPX+3*NH)*PRHU2
+ StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
C Potential intercatchment semi-exchange
! speed-up
- Rr = X(1)/Param(3)
+ 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((X(1)+X(8)+EXCH).LT.0) AEXCH1=-X(1)-X(8)
- X(1)=X(1)+X(8)+EXCH
- IF(X(1).LT.0.)X(1)=0.
-
+ 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 = X(1)/Param(3)
+ Rr = St(2)/Param(3)
Rr = Rr * Rr
Rr = Rr * Rr
- QR=X(1)*(1.-1./SQRT(SQRT(1.+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
-
- X(1)=X(1)-QR
+ St(2)=St(2)-QR
C Runoff from direct branch QD
AEXCH2=EXCH
- IF((X(8+NH)+EXCH).LT.0) AEXCH2=-X(8+NH)
- QD=MAX(0.,X(8+NH)+EXCH)
+ IF((StUH2(1)+EXCH).LT.0) AEXCH2=-StUH2(1)
+ QD=MAX(0.,StUH2(1)+EXCH)
C Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
C Variables storage
- MISC( 1)=E ! PE ! potential evapotranspiration [mm/d]
- MISC( 2)=P1 ! Precip ! total precipitation [mm/d]
- MISC( 3)=X(2) ! Prod ! production store level (X(2)) [mm]
- MISC( 4)=AE ! AE ! actual evapotranspiration [mm/d]
- MISC( 5)=PERC ! Perc ! percolation (PERC) [mm]
- MISC( 6)=PR ! PR ! PR=PN-PS+PERC [mm]
- MISC( 7)=X(8) ! Q9 ! outflow from HU1 (Q9) [mm/d]
- MISC( 8)=X(8+NH) ! Q1 ! outflow from HU2 (Q1) [mm/d]
- MISC( 9)=X(1) ! Rout ! routing store level (X(1)) [mm]
- MISC(10)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/d]
- MISC(11)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/d]
- MISC(12)=QR ! QR ! outflow from routing store (QR) [mm/d]
- MISC(13)=QD ! QD ! outflow from HU2 branch after exchange (QD) [mm/d]
- MISC(14)=Q ! Qsim ! outflow at catchment outlet [mm/d]
+ 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)=AE ! AE ! actual evapotranspiration [mm/day]
+ MISC( 5)=PERC ! Perc ! percolation (PERC) [mm/day]
+ MISC( 6)=PR ! PR ! PR=PN-PS+PERC [mm/day]
+ MISC( 7)=StUH1(1) ! Q9 ! outflow from UH1 (Q9) [mm/day]
+ MISC( 8)=StUH2(1) ! Q1 ! outflow from UH2 (Q1) [mm/day]
+ MISC( 9)=St(2) ! Rout ! routing store level (St(2)) [mm]
+ MISC(10)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
+ MISC(11)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
+ MISC(12)=QR ! QR ! outflow from routing store (QR) [mm/day]
+ MISC(13)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/day]
+ MISC(14)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f
index b8fb3d02..3dacbb39 100644
--- a/src/frun_GR5J.f
+++ b/src/frun_GR5J.f
@@ -3,20 +3,20 @@
SUBROUTINE frun_GR5J(
!inputs
& LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm]
- & InputsPE , ! [double] input series PE [mm]
- & NParam , ! [integer] number of model parameter
+ & 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 used for model initialising
- & StateStart , ! [double] state variables used when the model run starts (reservoir levels [mm] and HU)
+ & 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 (reservoir levels [mm] and HU)
+ & 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
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_GR5J
Implicit None
@@ -31,44 +31,49 @@
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
!parameters, internal states and variables
- integer NPX,NH,NMISC
- parameter (NPX=14,NH=20,NMISC=30)
- doubleprecision X(5*NH+7),XV(3*NPX+5*NH)
+ 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
!--------------------------------------------------------------
- !Initialisations
+ !Initializations
!--------------------------------------------------------------
!initilisation of model states to zero
- X=0.
- XV=0.
+ St=0.
+ StUH2=0.
!initilisation of model states using StateStart
- DO I=1,3*NH
- X(I)=StateStart(I)
+ St(1) = StateStart(1)
+ St(2) = StateStart(2)
+
+ DO I=1,2*NH
+ StUH2(I)=StateStart(7+I)
ENDDO
!parameter values
!Param(1) : production store capacity (X1 - PROD) [mm]
- !Param(2) : intercatchment exchange constant (X2 - CES1) [mm/d]
+ !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) [d]
- !Param(5) : intercatchment exchange constant (X5 - CES2) [-]
+ !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 HU1(XV,Param(4),D)
- CALL HU2(XV,Param(4),D)
+ CALL UH2(OrdUH2,Param(4),D)
- !initialisation of model outputs
+ !initialization of model outputs
Q = -999.999
MISC = -999.999
-c StateEnd = -999.999 !initialisation made in R
-c Outputs = -999.999 !initialisation made in R
+c StateEnd = -999.999 !initialization made in R
+c Outputs = -999.999 !initialization made in R
@@ -80,16 +85,18 @@ c Outputs = -999.999 !initialisation made in R
E =InputsPE(k)
c Q = -999.999
c MISC = -999.999
- !model run on one time-step
- CALL MOD_GR5J(X,XV,Param,P1,E,Q,MISC)
+ !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
- DO K=1,3*NH
- StateEnd(K)=X(K)
+ StateEnd(1) = St(1)
+ StateEnd(2) = St(2)
+ DO K=1,2*NH
+ StateEnd(7+K)=StUH2(K)
ENDDO
RETURN
@@ -106,57 +113,56 @@ c###############################################################################
C**********************************************************************
- SUBROUTINE MOD_GR5J(X,XV,Param,P1,E,Q,MISC)
-C Run on a single time-step with the GR5J model
+ 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 X Vector of model states at the beginning of the time-step [mm]
-C XV Vector of model states at the beginning of the time-step [mm]
-C Param Vector of model parameters [mixed units]
-C P1 Value of rainfall during the time-step [mm]
-C E Value of potential evapotranspiration during the time-step [mm]
+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 X Vector of model states at the end of the time-step [mm]
-C XV 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]
-C MISC Vector of model outputs for the time-step [mm]
+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 NPX,NH,NMISC,NParam
- PARAMETER (NPX=14,NH=20,NMISC=30)
+ INTEGER NH,NMISC,NParam
+ PARAMETER (NH=20,NMISC=30)
PARAMETER (NParam=5)
- DOUBLEPRECISION X(5*NH+7),XV(3*NPX+5*NH)
+ 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,PRHU1,PRHU2,EXCH,QR,QD
+ DOUBLEPRECISION PERC,PRHU1,PRHU2,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 Production store
+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 = X(2)/A
- ER=X(2)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+ 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
- IF(X(2).LT.ER) AE=X(2)+P1
- X(2)=X(2)-ER
+ St(1)=St(1)-ER
PR=0.
ELSE
EN=0.
@@ -164,92 +170,86 @@ C Production store
PN=P1-E
WS=PN/A
IF(WS.GT.13)WS=13.
-
! speed-up
TWS = tanHyp(WS)
- Sr = X(2)/A
+ 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
- X(2)=X(2)+PS
+ St(1)=St(1)+PS
ENDIF
C Percolation from production store
- IF(X(2).LT.0.)X(2)=0.
+ IF(St(1).LT.0.)St(1)=0.
! speed-up
! (9/4)**4 = 25.62890625
- Sr = X(2)/Param(1)
+ Sr = St(1)/Param(1)
Sr = Sr * Sr
Sr = Sr * Sr
- PERC=X(2)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
+ 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
- X(2)=X(2)-PERC
+ St(1)=St(1)-PERC
PR=PR+PERC
- PRHU1=PR*B
- PRHU2=PR*(1.-B)
-
-C Unit hydrograph HU1
- DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1)))
- X(7+K)=X(8+K)+XV(3*NPX+K)*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)*PR
ENDDO
- X(7+NH)=XV(3*NPX+NH)*PRHU1
+ StUH2(2*NH)=OrdUH2(2*NH)*PR
-C Unit hydrograph HU2
- DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1)))
- X(7+NH+K)=X(8+NH+K)+XV(3*NPX+NH+K)*PRHU2
- ENDDO
- X(7+3*NH)=XV(3*NPX+3*NH)*PRHU2
+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)*(X(1)/Param(3)-Param(5))
+ EXCH=Param(2)*(St(2)/Param(3)-Param(5))
C Routing store
AEXCH1=EXCH
- IF((X(1)+X(8)+EXCH).LT.0) AEXCH1=-X(1)-X(8)
- X(1)=X(1)+X(8)+EXCH
- IF(X(1).LT.0.)X(1)=0.
+ 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 = X(1)/Param(3)
+ Rr = St(2)/Param(3)
Rr = Rr * Rr
Rr = Rr * Rr
- QR=X(1)*(1.-1./SQRT(SQRT(1.+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
- X(1)=X(1)-QR
+ St(2)=St(2)-QR
C Runoff from direct branch QD
AEXCH2=EXCH
- IF((X(8+NH)+EXCH).LT.0) AEXCH2=-X(8+NH)
- QD=MAX(0.,X(8+NH)+EXCH)
+ IF((Q1+EXCH).LT.0) AEXCH2=-Q1
+ QD=MAX(0.,Q1+EXCH)
C Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
C Variables storage
- MISC( 1)=E ! PE ! potential evapotranspiration [mm/d]
- MISC( 2)=P1 ! Precip ! total precipitation [mm/d]
- MISC( 3)=X(2) ! Prod ! production store level (X(2)) [mm]
- MISC( 4)=AE ! AE ! actual evapotranspiration [mm/d]
- MISC( 5)=PERC ! Perc ! percolation (PERC) [mm]
- MISC( 6)=PR ! PR ! PR=PN-PS+PERC [mm]
- MISC( 7)=X(8) ! Q9 ! outflow from HU1 (Q9) [mm/d]
- MISC( 8)=X(8+NH) ! Q1 ! outflow from HU2 (Q1) [mm/d]
- MISC( 9)=X(1) ! Rout ! routing store level (X(1)) [mm]
- MISC(10)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/d]
- MISC(11)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/d]
- MISC(12)=QR ! QR ! outflow from routing store (QR) [mm/d]
- MISC(13)=QD ! QD ! outflow from HU2 branch after exchange (QD) [mm/d]
- MISC(14)=Q ! Qsim ! outflow at catchment outlet [mm/d]
+ 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)=AE ! AE ! actual evapotranspiration [mm/day]
+ MISC( 5)=PERC ! Perc ! percolation (PERC) [mm/day]
+ MISC( 6)=PR ! PR ! PR=PN-PS+PERC [mm/day]
+ MISC( 7)=Q9 ! Q9 ! outflow from first branch (Q9) [mm/day]
+ MISC( 8)=Q1 ! Q1 ! outflow from second branch (Q1) [mm/day]
+ MISC( 9)=St(2) ! Rout ! routing store level (St(2)) [mm]
+ MISC(10)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
+ MISC(11)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
+ MISC(12)=QR ! QR ! outflow from routing store (QR) [mm/day]
+ MISC(13)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/day]
+ MISC(14)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f
index 06ac8253..df27a4cd 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f
@@ -3,20 +3,20 @@
SUBROUTINE frun_GR6J(
!inputs
& LInputs , ! [integer] length of input and output series
- & InputsPrecip , ! [double] input series of total precipitation [mm]
- & InputsPE , ! [double] input series PE [mm]
- & NParam , ! [integer] number of model parameter
+ & 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 used for model initialising
- & StateStart , ! [double] state variables used when the model run starts (reservoir levels [mm] and HU)
+ & 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 (reservoir levels [mm] and HU)
+ & 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
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_GR6J
Implicit None
@@ -31,45 +31,56 @@
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
!parameters, internal states and variables
- integer NPX,NH,NMISC
- parameter (NPX=14,NH=20,NMISC=30)
- doubleprecision X(5*NH+7),XV(3*NPX+5*NH)
+ 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
!--------------------------------------------------------------
- !Initialisations
+ !Initializations
!--------------------------------------------------------------
!initilisation of model states to zero
- X=0.
- XV=0.
-
- !initilisation of model states using StateStart
- DO I=1,3*NH
- X(I)=StateStart(I)
+ 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 constant (X2 - CES1) [mm/d]
+ !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) [d]
- !Param(5) : intercatchment exchange constant (X5 - CES2) [-]
- !Param(6) : time constant of exponential store (X6 - EXP) [d]
+ !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 HU1(XV,Param(4),D)
- CALL HU2(XV,Param(4),D)
-
- !initialisation of model outputs
+ 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 !initialisation made in R
-c Outputs = -999.999 !initialisation made in R
+c StateEnd = -999.999 !initialization made in R
+c Outputs = -999.999 !initialization made in R
@@ -81,16 +92,22 @@ c Outputs = -999.999 !initialisation made in R
E =InputsPE(k)
c Q = -999.999
c MISC = -999.999
- !model run on one time-step
- CALL MOD_GR6J(X,XV,Param,P1,E,Q,MISC)
+ !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
- DO K=1,3*NH
- StateEnd(K)=X(K)
+ 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
@@ -107,25 +124,31 @@ c###############################################################################
C**********************************************************************
- SUBROUTINE MOD_GR6J(X,XV,Param,P1,E,Q,MISC)
-C Run on a single time-step with the GR6J model
+ 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 X Vector of model states at the beginning of the time-step [mm]
-C XV Vector of model states at the beginning of the time-step [mm]
-C Param Vector of model parameters [mixed units]
-C P1 Value of rainfall during the time-step [mm]
-C E Value of potential evapotranspiration during the time-step [mm]
+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 X Vector of model states at the end of the time-step [mm]
-C XV 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]
-C MISC Vector of model outputs for the time-step [mm]
+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]
+C MISC Vector of model outputs for the time step [mm]
C**********************************************************************
Implicit None
- INTEGER NPX,NH,NMISC,NParam
- PARAMETER (NPX=14,NH=20,NMISC=30)
+ INTEGER NH,NMISC,NParam
+ PARAMETER (NH=20,NMISC=30)
PARAMETER (NParam=6)
- DOUBLEPRECISION X(5*NH+7),XV(3*NPX+5*NH)
+ DOUBLEPRECISION St(3),StUH1(NH)
+ DOUBLEPRECISION StUH2(2*NH),OrdUH1(NH),OrdUH2(2*NH)
DOUBLEPRECISION Param(NParam)
DOUBLEPRECISION MISC(NMISC)
DOUBLEPRECISION P1,E,Q
@@ -136,7 +159,6 @@ C**********************************************************************
DATA B/0.9/
DATA C/0.4/
-
DOUBLEPRECISION TWS, Sr, Rr ! speed-up
A=Param(1)
@@ -151,14 +173,13 @@ C Production store
! speed-up
TWS = tanHyp(WS)
- Sr = X(2)/A
- ER=X(2)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+ 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
- IF(X(2).LT.ER) AE=X(2)+P1
- X(2)=X(2)-ER
+ St(1)=St(1)-ER
PR=0.
ELSE
EN=0.
@@ -169,72 +190,73 @@ C Production store
! speed-up
TWS = tanHyp(WS)
- Sr = X(2)/A
+ 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
- X(2)=X(2)+PS
+ St(1)=St(1)+PS
ENDIF
C Percolation from production store
- IF(X(2).LT.0.)X(2)=0.
+ IF(St(1).LT.0.)St(1)=0.
! speed-up
! (9/4)**4 = 25.62890625
- Sr = X(2)/Param(1)
+ Sr = St(1)/Param(1)
Sr = Sr * Sr
Sr = Sr * Sr
- PERC=X(2)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
+ 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
- X(2)=X(2)-PERC
+ 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 Unit hydrograph HU1
- DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1)))
- X(7+K)=X(8+K)+XV(3*NPX+K)*PRHU1
+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
- X(7+NH)=XV(3*NPX+NH)*PRHU1
+ StUH1(NH)=OrdUH1(NH)*PRHU1
-C Unit hydrograph HU2
- DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1)))
- X(7+NH+K)=X(8+NH+K)+XV(3*NPX+NH+K)*PRHU2
+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
- X(7+3*NH)=XV(3*NPX+3*NH)*PRHU2
+ StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
C Potential intercatchment semi-exchange
- EXCH=Param(2)*(X(1)/Param(3)-Param(5))
+ EXCH=Param(2)*(St(2)/Param(3)-Param(5))
C Routing store
AEXCH1=EXCH
- IF((X(1)+X(8)+EXCH).LT.0) AEXCH1=-X(1)-X(8)
- X(1)=X(1)+(1-C)*X(8)+EXCH
- IF(X(1).LT.0.)X(1)=0.
+ IF((St(2)+StUH1(1)+EXCH).LT.0) AEXCH1=-St(2)-StUH1(1)
+ St(2)=St(2)+(1-C)*StUH1(1)+EXCH
+ IF(St(2).LT.0.)St(2)=0.
! speed-up
- Rr = X(1)/Param(3)
+ Rr = St(2)/Param(3)
Rr = Rr * Rr
Rr = Rr * Rr
- QR=X(1)*(1.-1./SQRT(SQRT(1.+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
- X(1)=X(1)-QR
+ St(2)=St(2)-QR
C Update of exponential store
- X(6)=X(6)+C*X(8)+EXCH
- AR=X(6)/Param(6)
+ 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
- QR1=X(6)+Param(6)/EXP(AR)
+ QR1=St(3)+Param(6)/EXP(AR)
GOTO 3
ENDIF
@@ -246,34 +268,34 @@ C Update of exponential store
QR1=Param(6)*LOG(EXP(AR)+1.)
3 CONTINUE
- X(6)=X(6)-QR1
+ St(3)=St(3)-QR1
C Runoff from direct branch QD
AEXCH2=EXCH
- IF((X(8+NH)+EXCH).LT.0) AEXCH2=-X(8+NH)
- QD=MAX(0.,X(8+NH)+EXCH)
+ IF((StUH2(1)+EXCH).LT.0) AEXCH2=-StUH2(1)
+ QD=MAX(0.,StUH2(1)+EXCH)
C Total runoff
Q=QR+QD+QR1
IF(Q.LT.0.) Q=0.
C Variables storage
- MISC( 1)=E ! PE ! potential evapotranspiration [mm/d]
- MISC( 2)=P1 ! Precip ! total precipitation [mm/d]
- MISC( 3)=X(2) ! Prod ! production store level (X(2)) [mm]
- MISC( 4)=AE ! AE ! actual evapotranspiration [mm/d]
- MISC( 5)=PERC ! Perc ! percolation (PERC) [mm]
- MISC( 6)=PR ! PR ! PR=PN-PS+PERC [mm]
- MISC( 7)=X(8) ! Q9 ! outflow from HU1 (Q9) [mm/d]
- MISC( 8)=X(8+NH) ! Q1 ! outflow from HU2 (Q1) [mm/d]
- MISC( 9)=X(1) ! Rout ! routing store level (X(1)) [mm]
- MISC(10)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/d]
- MISC(11)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/d]
- MISC(12)=QR ! QR ! outflow from routing store (QR) [mm/d]
- MISC(13)=QR1 ! QR1 ! outflow from exponential store (QR1) [mm/d]
- MISC(14)=X(6) ! Exp ! exponential store level (X(6)) (negative) [mm]
- MISC(15)=QD ! QD ! outflow from HU2 branch after exchange (QD) [mm/d]
- MISC(16)=Q ! Qsim ! outflow at catchment outlet [mm/d]
+ 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)=AE ! AE ! actual evapotranspiration [mm/day]
+ MISC( 5)=PERC ! Perc ! percolation (PERC) [mm/day]
+ MISC( 6)=PR ! PR ! PR=PN-PS+PERC [mm/day]
+ MISC( 7)=StUH1(1) ! Q9 ! outflow from UH1 (Q9) [mm/day]
+ MISC( 8)=StUH2(1) ! Q1 ! outflow from UH2 (Q1) [mm/day]
+ MISC( 9)=St(2) ! Rout ! routing store level (St(2)) [mm]
+ MISC(10)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
+ MISC(11)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
+ MISC(12)=QR ! QR ! outflow from routing store (QR) [mm/day]
+ MISC(13)=QR1 ! QR1 ! outflow from exponential store (QR1) [mm/day]
+ MISC(14)=St(3) ! Exp ! exponential store level (St(3)) (negative) [mm]
+ MISC(15)=QD ! QD ! outflow from UH2 branch after exchange (QD) [mm/day]
+ MISC(16)=Q ! Qsim ! simulated outflow at catchment outlet [mm/day]
ENDSUBROUTINE
--
GitLab