diff --git a/DESCRIPTION b/DESCRIPTION
index aa8950b0408e46697b9e1b6b541a9266c652f9aa..a136852130f107196cf0d7051ab96a8bebeb4075 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.3.2.58
-Date: 2019-11-19
+Version: 1.3.2.59
+Date: 2019-11-20
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@irstea.fr"),
diff --git a/NEWS.md b/NEWS.md
index d5fe4de40df2e023d6736321e06bd7a83b7f4552..d6a5478552d99f2f2c8bcaae8bfc7504a5017f33 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
-### 1.3.2.58 Release Notes (2019-11-19)
+### 1.3.2.59 Release Notes (2019-11-19)
#### New features
@@ -16,10 +16,13 @@
- Fixed bug in <code>plot.OutputsModel()</code>. The function does not return any error message when <code>log_scale = TRUE</code>, <code>Qobs = NULL</code> and user want to draw flows time series.
- Fixed bug in <code>RunModel_*GR*()</code>. The functions do not return any error message anymore due to slightly negative values returned by GR4H, GR4J, GR5J or GR6J Fortran codes (the message was returned by <code>CreateIniStates()</code> when the final states were created). The <code>RunModel_*GR*()</code> functions now returns zero instead of these slightly negative values, except for the ExpStore where negatives values are allowed.
+#### Minor user-visible changes
+
+- Spurious flows set to <code>NA</code> into the <code>BasinObs</code> time series of the <code>L0123001</code> dataset.
#### CRAN-compatibility updates
-- Modification of the inputs in the <code>frun_cemaneige</code> Fortran subroutine to be consistent withe the <code>RunModel_CemaNeige()</code> function and the "airGR.c" file which registers native routines.
+- Cleaning of the Fortran codes (comment formatting).
____________________________________________________________________________________
diff --git a/src/frun_CEMANEIGE.f b/src/frun_CEMANEIGE.f
index 7b77eddb11a2573d75c4b0a1261d492ecc09bd9a..c4ca0aeaefb272a33635313f20d783d046c19ff5 100644
--- a/src/frun_CEMANEIGE.f
+++ b/src/frun_CEMANEIGE.f
@@ -3,7 +3,7 @@
SUBROUTINE frun_cemaneige(
- !inputs
+ ! inputs
& LInputs , ! [integer] length of input and output series
& InputsPrecip , ! [double] input series of total precipitation [mm/time step]
& InputsFracSolidPrecip, ! [double] input series of fraction of solid precipitation [0-1]
@@ -11,12 +11,12 @@
& 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 initialisation = 4
+ & NStates , ! [integer] number of state variables used for model initialization = 4
& StateStart , ! [double] state variables used when the model run starts
& IsHyst , ! [integer] whether we should use the linear hysteresis or not
& NOutputs , ! [integer] number of output series
& IndOutputs , ! [integer] indices of output series
- !outputs
+ ! outputs
& Outputs , ! [double] output series
& StateEnd ) ! [double] state variables at the end of the model run
@@ -25,7 +25,7 @@
Implicit None
- !### input and output variables
+ ! input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, intent(in) :: MeanAnSolidPrecip
doubleprecision, intent(in), dimension(LInputs) :: InputsPrecip
@@ -41,7 +41,7 @@
doubleprecision, intent(out), dimension(LInputs,NOutputs) ::
& Outputs
- !parameters, internal states and variables
+ ! parameters, internal states and variables
doubleprecision CTG,Kf
doubleprecision G,eTG,PliqAndMelt,dG,prct
doubleprecision Tmelt,Gthreshold,MinSpeed,Gacc,Glocalmax
@@ -53,20 +53,20 @@
!--------------------------------------------------------------
- !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)
Gratio=0.
PliqAndMelt=0.
- !setting parameter values
+ ! setting parameter values
CTG=Param(1)
Kf=Param(2)
IF (IsHystBool) THEN
@@ -82,31 +82,31 @@
Gthreshold=0.9*MeanAnSolidPrecip
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
@@ -128,10 +128,10 @@ 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 (IsHystBool) THEN
dG=G-Ginit !Melt in case of negative dG, accumulation otherwise
@@ -153,10 +153,10 @@ 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/time step]
IF(IndOutputs(I).EQ.2) Outputs(k,I)=Psol ! Psol ! observed solid precipitation [mm/time step]
diff --git a/src/frun_GR1A.f b/src/frun_GR1A.f
index 9d35501c2f1bd47c77e375200409a4a738147d00..00df0efe7400e8132b73eeab261a7c814a418475 100644
--- a/src/frun_GR1A.f
+++ b/src/frun_GR1A.f
@@ -1,7 +1,7 @@
SUBROUTINE frun_gr1a(
- !inputs
+ ! 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]
@@ -11,7 +11,7 @@
& 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
& Outputs , ! [double] output series
& StateEnd ) ! [double] state variables at the end of the model run (none here)
@@ -20,7 +20,7 @@
Implicit None
- !### input and output variables
+ ! input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, dimension(LInputs) :: InputsPrecip
doubleprecision, dimension(LInputs) :: InputsPE
@@ -30,7 +30,7 @@
integer, dimension(NOutputs) :: IndOutputs
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
- !parameters, internal states and variables
+ ! parameters, internal states and variables
integer NMISC
parameter (NMISC=3)
doubleprecision MISC(NMISC)
@@ -38,33 +38,33 @@
integer I,K
!--------------------------------------------------------------
- !Initializations
+ ! Initializations
!--------------------------------------------------------------
- !parameter values
- !Param(1) : PE adjustment factor [-]
+ ! parameter values
+ ! Param(1) : PE adjustment factor [-]
- !initialization of model outputs
+ ! initialization of model outputs
Q = -999.999
MISC = -999.999
-c Outputs = -999.999 !initialization made in R
+! Outputs = -999.999 ! initialization made in R
StateStart = -999.999 ! CRAN-compatibility updates
StateEnd = -999.999 ! CRAN-compatibility updates
!--------------------------------------------------------------
- !Time loop
+ ! 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
+! Q = -999.999
+! MISC = -999.999
+ ! model run on one time step
CALL MOD_GR1A(Param,P0,P1,E1,Q,MISC)
- !storage of outputs
+ ! storage of outputs
DO I=1,NOutputs
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
@@ -78,23 +78,23 @@ c MISC = -999.999
-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**********************************************************************
+! Calculation of streamflow on a single time step with the GR1A model
+! Inputs:
+! Param Vector of model parameters (Param(1) [-])
+! P0 Value of rainfall during the previous time step [mm/year]
+! P1 Value of rainfall during the current time step [mm/year]
+! E1 Value of potential evapotranspiration during the current time step [mm/year]
+! Outputs:
+! Q Value of simulated flow at the catchment outlet for the current time step [mm/year]
+! MISC Vector of model outputs for the time step [mm/year]
+!**********************************************************************
Implicit None
INTEGER NMISC,NParam
PARAMETER (NMISC=3)
@@ -105,14 +105,14 @@ C**********************************************************************
DOUBLEPRECISION tt ! speed-up
-C Runoff
+! 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
+ ! end speed-up
-C Variables storage
+! 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]
diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f
index 46f76304c06022e5795b86297f7c38ff65d43d6a..904d9befe5003a71a3fa6535d2b4ee83470390c3 100644
--- a/src/frun_GR2M.f
+++ b/src/frun_GR2M.f
@@ -1,7 +1,7 @@
SUBROUTINE frun_gr2m(
- !inputs
+ ! 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]
@@ -11,7 +11,7 @@
& 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
& Outputs , ! [double] output series
& StateEnd ) ! [double] state variables at the end of the model run (store levels [mm])
@@ -20,7 +20,7 @@
Implicit None
- !### input and output variables
+ ! input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, dimension(LInputs) :: InputsPrecip
doubleprecision, dimension(LInputs) :: InputsPE
@@ -30,7 +30,7 @@
integer, dimension(NOutputs) :: IndOutputs
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
- !parameters, internal states and variables
+ ! parameters, internal states and variables
integer NMISC
parameter (NMISC=30)
doubleprecision St(2)
@@ -39,36 +39,36 @@
integer I,K
!--------------------------------------------------------------
- !Initializations
+ ! Initializations
!--------------------------------------------------------------
- !initialization of model states to zero
+ ! initialization of model states to zero
St=0.
- !initilisation of model states using StateStart
+ ! 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 [-]
+ ! parameter values
+ ! Param(1) : production store capacity [mm]
+ ! Param(2) : groundwater exchange coefficient [-]
- !initialization of model outputs
+ ! 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
+! StateEnd = -999.999 ! initialization made in R
+! Outputs = -999.999 ! initialization made in R
!--------------------------------------------------------------
- !Time loop
+ ! Time loop
!--------------------------------------------------------------
DO k=1,LInputs
P =InputsPrecip(k)
E =InputsPE(k)
-c Q = -999.999
-c MISC = -999.999
+! Q = -999.999
+! MISC = -999.999
!model run on one time step
CALL MOD_GR2M(St,Param,P,E,Q,MISC)
!storage of outputs
@@ -76,7 +76,7 @@ c MISC = -999.999
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
- !model states at the end of the run
+ ! model states at the end of the run
StateEnd(1)=St(1)
StateEnd(2)=St(2)
@@ -88,24 +88,24 @@ c MISC = -999.999
-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**********************************************************************
+! Calculation of streamflow on a single time step (month) with the GR2M model
+! Inputs:
+! St Vector of model states at the beginning of the time step [mm]
+! Param Vector of model parameters (Param(1) [mm]; Param(2) [-])
+! P Value of rainfall during the time step [mm/month]
+! E Value of potential evapotranspiration during the time step [mm/month]
+! Outputs:
+! St Vector of model states at the end of the time step [mm]
+! Q Value of simulated flow at the catchment outlet for the time step [mm/month]
+! MISC Vector of model outputs for the time step [mm/month]
+!**********************************************************************
Implicit None
INTEGER NMISC,NParam
PARAMETER (NMISC=30)
@@ -119,7 +119,7 @@ C**********************************************************************
DOUBLEPRECISION TWS, Sr ! speed-up
-C Production store
+! Production store
WS=P/Param(1)
IF(WS.GT.13.)WS=13.
@@ -127,7 +127,7 @@ C Production store
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
+ ! end speed-up
P1=P+St(1)-S1
WS=E/Param(1)
@@ -137,35 +137,35 @@ C Production store
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
+ ! end speed-up
AE = S1 - S2
-C Percolation
+! 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
+ ! end speed-up
P2=S2-St(1)
P3=P1+P2
-C QR calculation (routing store)
+! QR calculation (routing store)
R1=St(2)+P3
-C Water exchange
+! Water exchange
R2=Param(2)*R1
EXCH = R2 - R1
-C Total runoff
+! Total runoff
Q=R2*R2/(R2+60.)
-C Updating store level
+! Updating store level
St(2)=R2-Q
-C Variables storage
+! 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]
diff --git a/src/frun_GR4H.f b/src/frun_GR4H.f
index 1b1193cfe4988bb28e03dbebfa741315903cc4c7..5ffe19e4cf3e35a267c70b8b73ec38bb2d77fd59 100644
--- a/src/frun_GR4H.f
+++ b/src/frun_GR4H.f
@@ -1,7 +1,7 @@
SUBROUTINE frun_gr4h(
- !inputs
+ ! 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]
@@ -11,7 +11,7 @@
& 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
& Outputs , ! [double] output series
& StateEnd ) ! [double] state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
@@ -20,7 +20,7 @@
Implicit None
- !### input and output variables
+ ! input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, dimension(LInputs) :: InputsPrecip
doubleprecision, dimension(LInputs) :: InputsPE
@@ -30,7 +30,7 @@
integer, dimension(NOutputs) :: IndOutputs
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
- !parameters, internal states and variables
+ ! parameters, internal states and variables
integer NH,NMISC
parameter (NH=480,NMISC=30)
doubleprecision St(2), StUH1(NH), StUH2(2*NH)
@@ -41,15 +41,15 @@
integer I,K
!--------------------------------------------------------------
- !Initializations
+ ! Initializations
!--------------------------------------------------------------
- !initilization of model states using StateStart
+ ! initialization of model states using StateStart
St=0.
StUH1=0.
StUH2=0.
- !initialization of model states using StateStart
+ ! initialization of model states using StateStart
St(1) = StateStart(1)
St(2) = StateStart(2)
DO I=1,NH
@@ -59,11 +59,11 @@
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]
+ ! 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.
@@ -73,30 +73,30 @@
CALL UH1_H(OrdUH1,Param(4),D)
CALL UH2_H(OrdUH2,Param(4),D)
- !initialization of model outputs
+ ! 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
+! StateEnd = -999.999 !initialization made in R
+! Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
- !Time loop
+ ! 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
+! 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
+ ! storage of outputs
DO I=1,NOutputs
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
- !model states at the end of the run
+ ! model states at the end of the run
StateEnd(1) = St(1)
StateEnd(2) = St(2)
DO K=1,NH
@@ -114,31 +114,31 @@ c MISC = -999.999
-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**********************************************************************
+! Run on a single time step with the GR4H model
+! Inputs:
+! St Vector of model states in stores at the beginning of the time step [mm]
+! StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
+! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+! OrdUH1 Vector of ordinates in UH1 [-]
+! OrdUH2 Vector of ordinates in UH2 [-]
+! Param Vector of model parameters [various units]
+! P1 Value of rainfall during the time step [mm]
+! E Value of potential evapotranspiration during the time step [mm]
+! Outputs:
+! St Vector of model states in stores at the beginning of the time step [mm]
+! StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
+! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+! Q Value of simulated flow at the catchment outlet for the time step [mm/hour]
+! MISC Vector of model outputs for the time step [mm/hour]
+!**********************************************************************
Implicit None
INTEGER NH,NMISC,NParam
PARAMETER (NH=480,NMISC=30)
@@ -160,7 +160,7 @@ C**********************************************************************
A=Param(1)
-C Interception and production store
+! Interception and production store
IF(P1.LE.E) THEN
EN=E-P1
PN=0.
@@ -172,7 +172,7 @@ C Interception and production store
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
+ ! end speed-up
AE=ER+P1
St(1)=St(1)-ER
@@ -189,13 +189,13 @@ C Interception and production store
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
+ ! end speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
-C Percolation from production store
+! Percolation from production store
IF(St(1).LT.0.)St(1)=0.
! speed-up
@@ -205,36 +205,36 @@ C Percolation from production store
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
+ ! end speed-up
St(1)=St(1)-PERC
PR=PR+PERC
-C Split of effective rainfall into the two routing components
+! Split of effective rainfall into the two routing components
PRHU1=PR*B
PRHU2=PR*(1.-B)
-C Convolution of unit hydrograph UH1
+! 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
+! 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
+! 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
+ ! end speed-up
-C Routing store
+! Routing store
AEXCH1=EXCH
IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
St(2)=St(2)+StUH1(1)+EXCH
@@ -246,20 +246,20 @@ C Routing store
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
+ ! end speed-up
St(2)=St(2)-QR
-C Runoff from direct branch QD
+! 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
+! Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
-C Variables storage
+! 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]
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f
index 892d8c07cb04315dcd1c1438597861fe21121300..14fb0c1d41405e61880b8f6f38329dc015e32ff3 100644
--- a/src/frun_GR4J.f
+++ b/src/frun_GR4J.f
@@ -1,7 +1,7 @@
SUBROUTINE frun_gr4j(
- !inputs
+ ! 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]
@@ -11,7 +11,7 @@
& 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
& Outputs , ! [double] output series
& StateEnd ) ! [double] state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
@@ -20,7 +20,7 @@
Implicit None
- !### input and output variables
+ ! input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, dimension(LInputs) :: InputsPrecip
doubleprecision, dimension(LInputs) :: InputsPE
@@ -30,7 +30,7 @@
integer, dimension(NOutputs) :: IndOutputs
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
- !parameters, internal states and variables
+ ! parameters, internal states and variables
integer NH,NMISC
parameter (NH=20,NMISC=30)
doubleprecision St(2), StUH1(NH), StUH2(2*NH)
@@ -41,15 +41,15 @@
integer I,K
!--------------------------------------------------------------
- !Initializations
+ ! Initializations
!--------------------------------------------------------------
- !initialization of model states to zero
+ ! initialization of model states to zero
St=0.
StUH1=0.
StUH2=0.
- !initialization of model states using StateStart
+ ! initialization of model states using StateStart
St(1) = StateStart(1)
St(2) = StateStart(2)
DO I=1,NH
@@ -59,13 +59,13 @@
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]
+ ! 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
+ ! computation of UH ordinates
OrdUH1 = 0.
OrdUH2 = 0.
@@ -73,30 +73,30 @@
CALL UH1(OrdUH1,Param(4),D)
CALL UH2(OrdUH2,Param(4),D)
- !initialization of model outputs
+ ! 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
+! StateEnd = -999.999 !initialization made in R
+! Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
- !Time loop
+ ! 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
+! 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
+ ! storage of outputs
DO I=1,NOutputs
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
- !model states at the end of the run
+ ! model states at the end of the run
StateEnd(1) = St(1)
StateEnd(2) = St(2)
DO K=1,NH
@@ -114,31 +114,31 @@ c MISC = -999.999
-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**********************************************************************
+! Run on a single time step with the GR4J model
+! Inputs:
+! St Vector of model states in stores at the beginning of the time step [mm]
+! StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
+! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+! OrdUH1 Vector of ordinates in UH1 [-]
+! OrdUH2 Vector of ordinates in UH2 [-]
+! Param Vector of model parameters [various units]
+! P1 Value of rainfall during the time step [mm/day]
+! E Value of potential evapotranspiration during the time step [mm/day]
+! Outputs:
+! St Vector of model states in stores at the end of the time step [mm]
+! StUH1 Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
+! StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
+! Q Value of simulated flow at the catchment outlet for the time step [mm/day]
+! MISC Vector of model outputs for the time step [mm/day]
+!**********************************************************************
Implicit None
INTEGER NH,NMISC,NParam
PARAMETER (NH=20,NMISC=30)
@@ -159,7 +159,7 @@ C**********************************************************************
A=Param(1)
-C Interception and production store
+! Interception and production store
IF(P1.LE.E) THEN
EN=E-P1
PN=0.
@@ -170,7 +170,7 @@ C Interception and production store
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
+ ! end speed-up
AE=ER+P1
St(1)=St(1)-ER
PS=0.
@@ -186,12 +186,12 @@ C Interception and production store
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
+ ! end speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
-C Percolation from production store
+! Percolation from production store
IF(St(1).LT.0.)St(1)=0.
! speed-up
! (9/4)**4 = 25.62891
@@ -200,35 +200,35 @@ C Percolation from production store
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
+ ! end speed-up
St(1)=St(1)-PERC
PR=PR+PERC
-C Split of effective rainfall into the two routing components
+! Split of effective rainfall into the two routing components
PRHU1=PR*B
PRHU2=PR*(1.-B)
-C Convolution of unit hydrograph UH1
+! 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
+! 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
+! 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
+ ! end speed-up
-C Routing store
+! Routing store
AEXCH1=EXCH
IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
St(2)=St(2)+StUH1(1)+EXCH
@@ -239,19 +239,19 @@ C Routing store
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
+ ! end speed-up
St(2)=St(2)-QR
-C Runoff from direct branch QD
+! 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
+! Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
-C Variables storage
+! 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]
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f
index f2750b629fef670b5dbee595e0b5f63c09a7f428..c9be386bafb934b81342e8c5b427a3708c37adc8 100644
--- a/src/frun_GR5J.f
+++ b/src/frun_GR5J.f
@@ -1,7 +1,7 @@
SUBROUTINE frun_gr5j(
- !inputs
+ ! 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]
@@ -11,7 +11,7 @@
& 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
& Outputs , ! [double] output series
& StateEnd ) ! [double] state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
@@ -20,7 +20,7 @@
Implicit None
- !### input and output variables
+ ! input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, dimension(LInputs) :: InputsPrecip
doubleprecision, dimension(LInputs) :: InputsPE
@@ -30,7 +30,7 @@
integer, dimension(NOutputs) :: IndOutputs
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
- !parameters, internal states and variables
+ ! parameters, internal states and variables
integer NH,NMISC
parameter (NH=20,NMISC=30)
doubleprecision St(2), StUH2(2*NH)
@@ -41,14 +41,14 @@
integer I,K
!--------------------------------------------------------------
- !Initializations
+ ! Initializations
!--------------------------------------------------------------
- !initilisation of model states to zero
+ ! initialization of model states to zero
St=0.
StUH2=0.
- !initilisation of model states using StateStart
+ ! initialization of model states using StateStart
St(1) = StateStart(1)
St(2) = StateStart(2)
@@ -56,12 +56,12 @@
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) [-]
+ ! 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.
@@ -69,30 +69,30 @@
D=2.5
CALL UH2(OrdUH2,Param(4),D)
- !initialization of model outputs
+ ! 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
+! StateEnd = -999.999 !initialization made in R
+! Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
- !Time loop
+ ! 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
+! 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
+ ! storage of outputs
DO I=1,NOutputs
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
- !model states at the end of the run
+ ! model states at the end of the run
StateEnd(1) = St(1)
StateEnd(2) = St(2)
DO K=1,2*NH
@@ -107,28 +107,28 @@ c MISC = -999.999
-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**********************************************************************
+! Run on a single time step with the GR5J model
+! Inputs:
+! St Vector of model states in stores at the beginning of the time step [mm]
+! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+! OrdUH2 Vector of ordinates in UH2 [-]
+! Param Vector of model parameters [various units]
+! P1 Value of rainfall during the time step [mm]
+! E Value of potential evapotranspiration during the time step [mm]
+! Outputs:
+! St Vector of model states in stores at the end of the time step [mm]
+! StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
+! Q Value of simulated flow at the catchment outlet for the time step [mm/day]
+! MISC Vector of model outputs for the time step [mm or mm/day]
+!**********************************************************************
Implicit None
INTEGER NH,NMISC,NParam
PARAMETER (NH=20,NMISC=30)
@@ -149,7 +149,7 @@ C**********************************************************************
A=Param(1)
-C Interception and production store
+! Interception and production store
IF(P1.LE.E) THEN
EN=E-P1
PN=0.
@@ -160,7 +160,7 @@ C Interception and production store
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
+ ! end speed-up
AE=ER+P1
St(1)=St(1)-ER
PS=0.
@@ -176,13 +176,13 @@ C Interception and production store
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
+ ! end speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
-C Percolation from production store
+! Percolation from production store
IF(St(1).LT.0.)St(1)=0.
! speed-up
@@ -192,26 +192,26 @@ C Percolation from production store
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
+ ! end speed-up
St(1)=St(1)-PERC
PR=PR+PERC
-C Convolution of unit hydrograph UH2
+! 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
+! 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
+! Potential intercatchment semi-exchange
EXCH=Param(2)*(St(2)/Param(3)-Param(5))
-C Routing store
+! Routing store
AEXCH1=EXCH
IF((St(2)+Q9+EXCH).LT.0.) AEXCH1=-St(2)-Q9
St(2)=St(2)+Q9+EXCH
@@ -223,20 +223,20 @@ C Routing store
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
+ ! end speed-up
St(2)=St(2)-QR
-C Runoff from direct branch QD
+! Runoff from direct branch QD
AEXCH2=EXCH
IF((Q1+EXCH).LT.0.) AEXCH2=-Q1
QD=MAX(0.d0,Q1+EXCH)
-C Total runoff
+! Total runoff
Q=QR+QD
IF(Q.LT.0.) Q=0.
-C Variables storage
+! 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]
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f
index c967d915eda7aa0a00449701eb28144c1eccb9c4..17927fd4dd436a89c847d1dc55311b7a0f40d935 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f
@@ -1,7 +1,7 @@
SUBROUTINE frun_gr6j(
- !inputs
+ ! 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]
@@ -11,7 +11,7 @@
& 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
& Outputs , ! [double] output series
& StateEnd ) ! [double] state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
@@ -20,7 +20,7 @@
Implicit None
- !### input and output variables
+ ! input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, dimension(LInputs) :: InputsPrecip
doubleprecision, dimension(LInputs) :: InputsPE
@@ -30,7 +30,7 @@
integer, dimension(NOutputs) :: IndOutputs
doubleprecision, dimension(LInputs,NOutputs) :: Outputs
- !parameters, internal states and variables
+ ! parameters, internal states and variables
integer NH,NMISC
parameter (NH=20,NMISC=30)
doubleprecision St(3), StUH1(NH), StUH2(2*NH)
@@ -41,15 +41,15 @@
integer I,K
!--------------------------------------------------------------
- !Initializations
+ ! Initializations
!--------------------------------------------------------------
- !initilisation of model states to zero
+ ! initialization of model states to zero
St=0.
StUH1=0.
StUH2=0.
- !initialization of model states using StateStart
+ ! initialization of model states using StateStart
St(1) = StateStart(1)
St(2) = StateStart(2)
St(3) = StateStart(3)
@@ -60,15 +60,15 @@
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]
+ ! 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
+ ! computation of HU ordinates
OrdUH1 = 0.
OrdUH2 = 0.
@@ -76,30 +76,30 @@
CALL UH1(OrdUH1,Param(4),D)
CALL UH2(OrdUH2,Param(4),D)
- !initialization of model outputs
+ ! 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
+! StateEnd = -999.999 !initialization made in R
+! Outputs = -999.999 !initialization made in R
!--------------------------------------------------------------
- !Time loop
+ ! 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
+! 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
+ ! storage of outputs
DO I=1,NOutputs
Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
- !model states at the end of the run
+ ! model states at the end of the run
StateEnd(1) = St(1)
StateEnd(2) = St(2)
StateEnd(3) = St(3)
@@ -118,31 +118,31 @@ c MISC = -999.999
-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**********************************************************************
+! Run on a single time step with the GR6J model
+! Inputs:
+! St Vector of model states in stores at the beginning of the time step [mm]
+! StUH1 Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
+! StUH2 Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+! OrdUH1 Vector of ordinates in UH1 [-]
+! OrdUH2 Vector of ordinates in UH2 [-]
+! Param Vector of model parameters [various units]
+! P1 Value of rainfall during the time step [mm]
+! E Value of potential evapotranspiration during the time step [mm]
+! Outputs:
+! St Vector of model states in stores at the end of the time step [mm]
+! StUH1 Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
+! StUH2 Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
+! Q Value of simulated flow at the catchment outlet for the time step [mm/day]
+! MISC Vector of model outputs for the time step [mm or mm/day]
+!**********************************************************************
Implicit None
INTEGER NH,NMISC,NParam
PARAMETER (NH=20,NMISC=30)
@@ -164,7 +164,7 @@ C**********************************************************************
A=Param(1)
-C Production store
+! Production store
IF(P1.LE.E) THEN
EN=E-P1
PN=0.
@@ -176,7 +176,7 @@ C Production store
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
+ ! end speed-up
AE=ER+P1
St(1)=St(1)-ER
@@ -194,13 +194,13 @@ C Production store
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
+ ! end speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
-C Percolation from production store
+! Percolation from production store
IF(St(1).LT.0.)St(1)=0.
! speed-up
! (9/4)**4 = 25.62890625
@@ -209,32 +209,32 @@ C Percolation from production store
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
+ ! end speed-up
St(1)=St(1)-PERC
PR=PR+PERC
-C Split of effective rainfall into the two routing components
+! Split of effective rainfall into the two routing components
PRUH1=PR*B
PRUH2=PR*(1.-B)
-C Convolution of unit hydrograph UH1
+! 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
+! 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
+! Potential intercatchment semi-exchange
EXCH=Param(2)*(St(2)/Param(3)-Param(5))
-C Routing store
+! 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
@@ -246,11 +246,11 @@ C Routing store
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
+ ! end speed-up
St(2)=St(2)-QR
-C Update of exponential store
+! Update of exponential store
St(3)=St(3)+C*StUH1(1)+EXCH
AR=St(3)/Param(6)
IF(AR.GT.33.)AR=33.
@@ -271,16 +271,16 @@ C Update of exponential store
St(3)=St(3)-QRExp
-C Runoff from direct branch QD
+! 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
+! Total runoff
Q=QR+QD+QRExp
IF(Q.LT.0.) Q=0.
-C Variables storage
+! 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]
diff --git a/src/utils_D.f b/src/utils_D.f
index b87484d8568d9d689d26e94241b3e57016d52e35..eaf0d6ff9cd49091dfc13bd34b4820721fa29d67 100644
--- a/src/utils_D.f
+++ b/src/utils_D.f
@@ -1,12 +1,12 @@
-C**********************************************************************
+!**********************************************************************
SUBROUTINE UH1(OrdUH1,C,D)
-C Computation of ordinates of GR unit hydrograph UH1 using successive differences on the S curve SS1
-C Inputs:
-C C: time constant
-C D: exponent
-C Outputs:
-C OrdUH1: NH ordinates of discrete hydrograph
-C**********************************************************************
+! Computation of ordinates of GR unit hydrograph UH1 using successive differences on the S curve SS1
+! Inputs:
+! C: time constant
+! D: exponent
+! Outputs:
+! OrdUH1: NH ordinates of discrete hydrograph
+!**********************************************************************
Implicit None
INTEGER NH
PARAMETER (NH=20)
@@ -20,15 +20,15 @@ C**********************************************************************
ENDSUBROUTINE
-C**********************************************************************
+!**********************************************************************
SUBROUTINE UH2(OrdUH2,C,D)
-C Computation of ordinates of GR unit hydrograph HU2 using successive differences on the S curve SS2
-C Inputs:
-C C: time constant
-C D: exponent
-C Outputs:
-C OrdUH2: 2*NH ordinates of discrete hydrograph
-C**********************************************************************
+! Computation of ordinates of GR unit hydrograph HU2 using successive differences on the S curve SS2
+! Inputs:
+! C: time constant
+! D: exponent
+! Outputs:
+! OrdUH2: 2*NH ordinates of discrete hydrograph
+!**********************************************************************
Implicit None
INTEGER NH
PARAMETER (NH=20)
@@ -42,16 +42,16 @@ C**********************************************************************
ENDSUBROUTINE
-C**********************************************************************
+!**********************************************************************
FUNCTION SS1(I,C,D)
-C Values of the S curve (cumulative HU curve) of GR unit hydrograph UH1
-C Inputs:
-C C: time constant
-C D: exponent
-C I: time-step
-C Outputs:
-C SS1: Values of the S curve for I
-C**********************************************************************
+! Values of the S curve (cumulative HU curve) of GR unit hydrograph UH1
+! Inputs:
+! C: time constant
+! D: exponent
+! I: time-step
+! Outputs:
+! SS1: Values of the S curve for I
+!**********************************************************************
Implicit None
DOUBLEPRECISION C,D,SS1
INTEGER I,FI
@@ -69,16 +69,16 @@ C**********************************************************************
ENDFUNCTION
-C**********************************************************************
+!**********************************************************************
FUNCTION SS2(I,C,D)
-C Values of the S curve (cumulative HU curve) of GR unit hydrograph UH2
-C Inputs:
-C C: time constant
-C D: exponent
-C I: time-step
-C Outputs:
-C SS2: Values of the S curve for I
-C**********************************************************************
+! Values of the S curve (cumulative HU curve) of GR unit hydrograph UH2
+! Inputs:
+! C: time constant
+! D: exponent
+! I: time-step
+! Outputs:
+! SS2: Values of the S curve for I
+!**********************************************************************
Implicit None
DOUBLEPRECISION C,D,SS2
INTEGER I,FI
@@ -100,10 +100,10 @@ C**********************************************************************
ENDFUNCTION
-C**********************************************************************
+!**********************************************************************
FUNCTION tanHyp(Val)
-C Computation of hyperbolic tangent
-C**********************************************************************
+! Computation of hyperbolic tangent
+!**********************************************************************
Implicit None
DOUBLEPRECISION Val,ValExp,tanHyp
diff --git a/src/utils_H.f b/src/utils_H.f
index 78352f3cc3ec5277e1a44f0a1b25e0e252610204..e3c2ada84d03cb62a29e943aa2b5b94bbf00e256 100644
--- a/src/utils_H.f
+++ b/src/utils_H.f
@@ -1,14 +1,14 @@
-C**********************************************************************
+!**********************************************************************
SUBROUTINE UH1_H(OrdUH1,C,D)
-C Computation of ordinates of GR unit hydrograph UH2 using successive differences on the S curve SS1
-C Inputs:
-C C: time constant
-C D: exponent
-C Outputs:
-C OrdUH1: NH ordinates of discrete hydrograph
-C**********************************************************************
+! Computation of ordinates of GR unit hydrograph UH2 using successive differences on the S curve SS1
+! Inputs:
+! C: time constant
+! D: exponent
+! Outputs:
+! OrdUH1: NH ordinates of discrete hydrograph
+!**********************************************************************
Implicit None
INTEGER NH
PARAMETER (NH=480)
@@ -22,15 +22,15 @@ C**********************************************************************
ENDSUBROUTINE
-C**********************************************************************
+!**********************************************************************
SUBROUTINE UH2_H(OrdUH2,C,D)
-C Computation of ordinates of GR unit hydrograph UH2 using successive differences on the S curve SS2
-C Inputs:
-C C: time constant
-C D: exponent
-C Outputs:
-C OrdUH2: 2*NH ordinates of discrete hydrograph
-C**********************************************************************
+! Computation of ordinates of GR unit hydrograph UH2 using successive differences on the S curve SS2
+! Inputs:
+! C: time constant
+! D: exponent
+! Outputs:
+! OrdUH2: 2*NH ordinates of discrete hydrograph
+!**********************************************************************
Implicit None
INTEGER NH
PARAMETER (NH=480)
@@ -44,16 +44,16 @@ C**********************************************************************
ENDSUBROUTINE
-C**********************************************************************
+!**********************************************************************
FUNCTION SS1_H(I,C,D)
-C Values of the S curve (cumulative HU curve) of GR unit hydrograph HU1
-C Inputs:
-C C: time constant
-C D: exponent
-C I: time-step
-C Outputs:
-C SS1: Values of the S curve for I
-C**********************************************************************
+! Values of the S curve (cumulative HU curve) of GR unit hydrograph HU1
+! Inputs:
+! C: time constant
+! D: exponent
+! I: time-step
+! Outputs:
+! SS1: Values of the S curve for I
+!**********************************************************************
Implicit None
DOUBLEPRECISION C,D,SS1_H
INTEGER I,FI
@@ -71,16 +71,16 @@ C**********************************************************************
ENDFUNCTION
-C**********************************************************************
+!**********************************************************************
FUNCTION SS2_H(I,C,D)
-C Values of the S curve (cumulative HU curve) of GR unit hydrograph HU2
-C Inputs:
-C C: time constant
-C D: exponent
-C I: time-step
-C Outputs:
-C SS2: Values of the S curve for I
-C**********************************************************************
+! Values of the S curve (cumulative HU curve) of GR unit hydrograph HU2
+! Inputs:
+! C: time constant
+! D: exponent
+! I: time-step
+! Outputs:
+! SS2: Values of the S curve for I
+!**********************************************************************
Implicit None
DOUBLEPRECISION C,D,SS2_H
INTEGER I,FI