Commit 66712dbf authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Mise au propre des noms de variables

parent c63f3cfe
......@@ -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
......
......@@ -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
......
......@@ -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