From df8031309e6dc1526058436651e964a6e4c0c6c9 Mon Sep 17 00:00:00 2001
From: Delaigue Olivier <olivier.delaigue@irstea.priv>
Date: Wed, 3 Apr 2019 15:38:36 +0200
Subject: [PATCH] v1.2.13.14 CLEAN: tabulation removed and variable statements
fixed in Fortran files
---
DESCRIPTION | 2 +-
NEWS.rmd | 8 ++-
src/frun_CEMANEIGE.f | 6 +--
src/frun_GR1A.f | 2 +-
src/frun_GR2M.f | 38 +++++++-------
src/frun_GR4H.f | 74 +++++++++++++-------------
src/frun_GR4J.f | 56 ++++++++++----------
src/frun_GR5J.f | 90 ++++++++++++++++----------------
src/frun_GR6J.f | 120 +++++++++++++++++++++----------------------
9 files changed, 200 insertions(+), 196 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 2e89f160..fabdb875 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.2.13.13
+Version: 1.2.13.14
Date: 2019-04-03
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
diff --git a/NEWS.rmd b/NEWS.rmd
index 80dcd671..74919746 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -13,8 +13,7 @@ output:
-### 1.2.13.13 Release Notes (2019-04-03)
-
+### 1.2.13.14 Release Notes (2019-04-03)
#### Deprecated and defunct
@@ -96,6 +95,11 @@ output:
- The order of authors has been updated in the DESCRIPTION and the CITATION files.
+
+#### CRAN-compatibility updates
+
+- Tabulation removed and variable statements fixed in Fortran files
+
____________________________________________________________________________________
diff --git a/src/frun_CEMANEIGE.f b/src/frun_CEMANEIGE.f
index 5d18df4d..23f7d26d 100644
--- a/src/frun_CEMANEIGE.f
+++ b/src/frun_CEMANEIGE.f
@@ -113,7 +113,7 @@ c Outputs = -999.999 !initialisation made in R
IF (IsHyst) THEN
IF (Potmelt.GT.0.) THEN
IF (G.LT.Glocalmax.AND.Gratio.EQ.1.) Glocalmax=G !Update in case of potential melt and G lower than Gseuil
- Gratio=MIN(G/Glocalmax,1.)
+ Gratio=MIN(G/Glocalmax,1.d0)
ENDIF
ELSE
IF(G.LT.Gthreshold) THEN
@@ -133,11 +133,11 @@ c Outputs = -999.999 !initialisation made in R
IF (dG.GT.0.) THEN
- Gratio = MIN(Gratio+(Psol-Melt)/Gacc,1.) !Psol - Melt = dG
+ Gratio = MIN(Gratio+(Psol-Melt)/Gacc,1.d0) !Psol - Melt = dG
IF (Gratio.EQ.1.) Glocalmax = Gthreshold
ENDIF
IF (dG.LT.0.) THEN
- Gratio=MIN(G/Glocalmax,1.)
+ Gratio=MIN(G/Glocalmax,1.d0)
ENDIF
ELSE
IF(G.LT.Gthreshold) THEN
diff --git a/src/frun_GR1A.f b/src/frun_GR1A.f
index b6655731..dc4d31a8 100644
--- a/src/frun_GR1A.f
+++ b/src/frun_GR1A.f
@@ -102,7 +102,7 @@ C**********************************************************************
DOUBLEPRECISION P0,P1,E1,Q
DOUBLEPRECISION tt ! speed-up
-
+
C Runoff
! speed-up
tt = (0.7*P1+0.3*P0)/Param(1)/E1
diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f
index 7a2102c1..a53b2c57 100644
--- a/src/frun_GR2M.f
+++ b/src/frun_GR2M.f
@@ -118,35 +118,35 @@ C**********************************************************************
DOUBLEPRECISION P1,P2,P3,R1,R2,AE,EXCH
DOUBLEPRECISION TWS, Sr, Rr ! speed-up
-
+
C Production store
WS=P/Param(1)
IF(WS.GT.13.)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- S1=(St(1)+Param(1)*TWS)/(1.+St(1)/Param(1)*TWS)
- ! S1=(X(1)+Param(1)*tanHyp(WS))/(1.+X(1)/Param(1)*tanHyp(WS))
- ! fin speed-up
-
- P1=P+St(1)-S1
+
+ ! speed-up
+ TWS = tanHyp(WS)
+ S1=(St(1)+Param(1)*TWS)/(1.+St(1)/Param(1)*TWS)
+ ! S1=(X(1)+Param(1)*tanHyp(WS))/(1.+X(1)/Param(1)*tanHyp(WS))
+ ! fin speed-up
+
+ P1=P+St(1)-S1
WS=E/Param(1)
IF(WS.GT.13.)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
+
+ ! speed-up
+ TWS = tanHyp(WS)
S2=S1*(1.-TWS)/(1.+(1.-S1/Param(1))*TWS)
! S2=S1*(1.-tanHyp(WS))/(1.+(1.-S1/Param(1))*tanHyp(WS))
- ! fin speed-up
- AE = S1 - S2
+ ! fin speed-up
+ AE = S1 - S2
C Percolation
- ! speed-up
- Sr = S2/Param(1)
- Sr = Sr * Sr * Sr + 1.
+ ! 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
+ ! fin speed-up
P2=S2-St(1)
P3=P1+P2
@@ -156,7 +156,7 @@ C QR calculation (routing store)
C Water exchange
R2=Param(2)*R1
- EXCH = R2 - R1
+ EXCH = R2 - R1
C Total runoff
Q=R2*R2/(R2+60.)
diff --git a/src/frun_GR4H.f b/src/frun_GR4H.f
index 4e711894..83732b42 100644
--- a/src/frun_GR4H.f
+++ b/src/frun_GR4H.f
@@ -50,8 +50,8 @@
StUH2=0.
!initialization of model states using StateStart
- St(1) = StateStart(1)
- St(2) = StateStart(2)
+ St(1) = StateStart(1)
+ St(2) = StateStart(2)
DO I=1,NH
StUH1(I)=StateStart(7+I)
ENDDO
@@ -66,9 +66,9 @@
!Param(4) : time constant of unit hydrograph (X4 - TB) [hour]
!computation of UH ordinates
- OrdUH1 = 0.
- OrdUH2 = 0.
-
+ OrdUH1 = 0.
+ OrdUH2 = 0.
+
D=1.25
CALL UH1_H(OrdUH1,Param(4),D)
CALL UH2_H(OrdUH2,Param(4),D)
@@ -97,8 +97,8 @@ c MISC = -999.999
ENDDO
ENDDO
!model states at the end of the run
- StateEnd(1) = St(1)
- StateEnd(2) = St(2)
+ StateEnd(1) = St(1)
+ StateEnd(2) = St(2)
DO K=1,NH
StateEnd(7+K)=StUH1(K)
ENDDO
@@ -156,7 +156,7 @@ C**********************************************************************
DATA B/0.9/
DOUBLEPRECISION TWS, Sr, Rr ! speed-up
-
+
A=Param(1)
@@ -166,14 +166,14 @@ C Interception and production store
PN=0.
WS=EN/A
IF(WS.GT.13.)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- Sr = St(1)/A
- ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+
+ ! 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
-
+ ! fin speed-up
+
AE=ER+P1
St(1)=St(1)-ER
PR=0.
@@ -183,14 +183,14 @@ C Interception and production store
PN=P1-E
WS=PN/A
IF(WS.GT.13.)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- Sr = St(1)/A
- PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
+
+ ! 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
-
+ ! fin speed-up
+
PR=PN-PS
St(1)=St(1)+PS
ENDIF
@@ -198,14 +198,14 @@ C Interception and production store
C Percolation from production store
IF(St(1).LT.0.)St(1)=0.
- ! speed-up
- ! (21/4)**4 = 759.69140625
- Sr = St(1)/Param(1)
- Sr = Sr * Sr
- Sr = Sr * Sr
+ ! 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
+ ! PERC=X(2)*(1.-(1.+(X(2)/(21./4.*Param(1)))**4.)**(-0.25))
+ ! fin speed-up
St(1)=St(1)-PERC
@@ -228,11 +228,11 @@ C Convolution of unit hydrograph UH2
StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
C Potential intercatchment semi-exchange
- ! speed-up
- Rr = St(2)/Param(3)
+ ! 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
+ ! fin speed-up
C Routing store
AEXCH1=EXCH
@@ -240,13 +240,13 @@ C Routing store
St(2)=St(2)+StUH1(1)+EXCH
IF(St(2).LT.0.)St(2)=0.
- ! speed-up
- Rr = St(2)/Param(3)
- Rr = Rr * Rr
- Rr = Rr * Rr
+ ! 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
+ ! fin speed-up
St(2)=St(2)-QR
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f
index cc44e959..99633faa 100644
--- a/src/frun_GR4J.f
+++ b/src/frun_GR4J.f
@@ -50,8 +50,8 @@
StUH2=0.
!initialization of model states using StateStart
- St(1) = StateStart(1)
- St(2) = StateStart(2)
+ St(1) = StateStart(1)
+ St(2) = StateStart(2)
DO I=1,NH
StUH1(I)=StateStart(7+I)
ENDDO
@@ -66,9 +66,9 @@
!Param(4) : time constant of unit hydrograph (X4 - TB) [day]
!computation of UH ordinates
- OrdUH1 = 0.
- OrdUH2 = 0.
-
+ OrdUH1 = 0.
+ OrdUH2 = 0.
+
D=2.5
CALL UH1(OrdUH1,Param(4),D)
CALL UH2(OrdUH2,Param(4),D)
@@ -97,8 +97,8 @@ c MISC = -999.999
ENDDO
ENDDO
!model states at the end of the run
- StateEnd(1) = St(1)
- StateEnd(2) = St(2)
+ StateEnd(1) = St(1)
+ StateEnd(2) = St(2)
DO K=1,NH
StateEnd(7+K)=StUH1(K)
ENDDO
@@ -165,15 +165,15 @@ C Interception and production store
PN=0.
WS=EN/A
IF(WS.GT.13.)WS=13.
- ! speed-up
+ ! 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
+ ! fin speed-up
AE=ER+P1
St(1)=St(1)-ER
- PS=0.
+ PS=0.
PR=0.
ELSE
EN=0.
@@ -181,26 +181,26 @@ C Interception and production store
PN=P1-E
WS=PN/A
IF(WS.GT.13.)WS=13.
- ! speed-up
+ ! 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
+ ! fin speed-up
PR=PN-PS
St(1)=St(1)+PS
ENDIF
C Percolation from production store
IF(St(1).LT.0.)St(1)=0.
- ! speed-up
- ! (9/4)**4 = 25.62891
- Sr = St(1)/Param(1)
- Sr = Sr * Sr
- Sr = Sr * Sr
+ ! speed-up
+ ! (9/4)**4 = 25.62891
+ Sr = St(1)/Param(1)
+ Sr = Sr * Sr
+ Sr = Sr * Sr
PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62891)))
- ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
- ! fin speed-up
+ ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
+ ! fin speed-up
St(1)=St(1)-PERC
PR=PR+PERC
@@ -222,24 +222,24 @@ C Convolution of unit hydrograph UH2
StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
C Potential intercatchment semi-exchange
- ! speed-up
- Rr = St(2)/Param(3)
+ ! 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
+ ! fin speed-up
C Routing store
AEXCH1=EXCH
IF((St(2)+StUH1(1)+EXCH).LT.0.) AEXCH1=-St(2)-StUH1(1)
St(2)=St(2)+StUH1(1)+EXCH
IF(St(2).LT.0.)St(2)=0.
- ! speed-up
- Rr = St(2)/Param(3)
- Rr = Rr * Rr
- Rr = Rr * Rr
+ ! 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
+ ! fin speed-up
St(2)=St(2)-QR
C Runoff from direct branch QD
@@ -264,7 +264,7 @@ C Variables storage
MISC(10)=StUH2(1) ! Q1 ! outflow from UH2 (Q1) [mm/day]
MISC(11)=St(2) ! Rout ! routing store level (St(2)) [mm]
MISC(12)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
- MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
+ MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
MISC(14)=AEXCH2 ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
MISC(15)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
MISC(16)=QR ! QR ! outflow from routing store (QR) [mm/day]
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f
index 19e1209c..74546f31 100644
--- a/src/frun_GR5J.f
+++ b/src/frun_GR5J.f
@@ -49,8 +49,8 @@
StUH2=0.
!initilisation of model states using StateStart
- St(1) = StateStart(1)
- St(2) = StateStart(2)
+ St(1) = StateStart(1)
+ St(2) = StateStart(2)
DO I=1,2*NH
StUH2(I)=StateStart(7+NH+I)
@@ -64,7 +64,7 @@
!Param(5) : intercatchment exchange threshold (X5 - CES2) [-]
!computation of HU ordinates
- OrdUH2 = 0.
+ OrdUH2 = 0.
D=2.5
CALL UH2(OrdUH2,Param(4),D)
@@ -93,8 +93,8 @@ c MISC = -999.999
ENDDO
ENDDO
!model states at the end of the run
- StateEnd(1) = St(1)
- StateEnd(2) = St(2)
+ StateEnd(1) = St(1)
+ StateEnd(2) = St(2)
DO K=1,2*NH
StateEnd(7+NH+K)=StUH2(K)
ENDDO
@@ -155,44 +155,44 @@ C Interception and production store
PN=0.
WS=EN/A
IF(WS.GT.13.)WS=13.
- ! speed-up
- TWS = tanHyp(WS)
- Sr = St(1)/A
- ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
- ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
- ! fin speed-up
- AE=ER+P1
- St(1)=St(1)-ER
- PS=0.
- PR=0.
+ ! speed-up
+ TWS = tanHyp(WS)
+ Sr = St(1)/A
+ ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+ ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+ ! fin speed-up
+ AE=ER+P1
+ St(1)=St(1)-ER
+ PS=0.
+ PR=0.
ELSE
- EN=0.
- AE=E
- PN=P1-E
- WS=PN/A
- IF(WS.GT.13.)WS=13.
- ! speed-up
- TWS = tanHyp(WS)
- Sr = St(1)/A
- PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
- ! PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
- ! fin speed-up
-
- PR=PN-PS
- St(1)=St(1)+PS
+ EN=0.
+ AE=E
+ PN=P1-E
+ WS=PN/A
+ IF(WS.GT.13.)WS=13.
+ ! speed-up
+ TWS = tanHyp(WS)
+ Sr = St(1)/A
+ PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
+ ! PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
+ ! fin speed-up
+
+ PR=PN-PS
+ St(1)=St(1)+PS
ENDIF
C Percolation from production store
IF(St(1).LT.0.)St(1)=0.
- ! speed-up
- ! (9/4)**4 = 25.62890625
- Sr = St(1)/Param(1)
- Sr = Sr * Sr
- Sr = Sr * Sr
+ ! speed-up
+ ! (9/4)**4 = 25.62890625
+ Sr = St(1)/Param(1)
+ Sr = Sr * Sr
+ Sr = Sr * Sr
PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
- ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
- ! fin speed-up
+ ! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
+ ! fin speed-up
St(1)=St(1)-PERC
@@ -200,7 +200,7 @@ C Percolation from production store
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
+ StUH2(K)=StUH2(K+1)+OrdUH2(K)*PR
ENDDO
StUH2(2*NH)=OrdUH2(2*NH)*PR
@@ -217,13 +217,13 @@ C Routing store
St(2)=St(2)+Q9+EXCH
IF(St(2).LT.0.)St(2)=0.
- ! speed-up
- Rr = St(2)/Param(3)
- Rr = Rr * Rr
- Rr = Rr * Rr
+ ! 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
+ ! fin speed-up
St(2)=St(2)-QR
@@ -245,11 +245,11 @@ C Variables storage
MISC( 6)=AE ! AE ! actual evapotranspiration [mm/day]
MISC( 7)=PERC ! Perc ! percolation (PERC) [mm/day]
MISC( 8)=PR ! PR ! PR=PN-PS+PERC [mm/day]
- MISC( 9)=Q9 ! Q9 ! outflow from UH1 (Q9) [mm/day]
- MISC(10)=Q1 ! Q1 ! outflow from UH2 (Q1) [mm/day]
+ MISC( 9)=Q9 ! Q9 ! outflow from UH1 (Q9) [mm/day]
+ MISC(10)=Q1 ! Q1 ! outflow from UH2 (Q1) [mm/day]
MISC(11)=St(2) ! Rout ! routing store level (St(2)) [mm]
MISC(12)=EXCH ! Exch ! potential semi-exchange between catchments (EXCH) [mm/day]
- MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
+ MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from branch 1 (AEXCH1) [mm/day]
MISC(14)=AEXCH2 ! AExch2 ! actual exchange between catchments from branch 2 (AEXCH2) [mm/day]
MISC(15)=AEXCH1+AEXCH2 ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
MISC(16)=QR ! QR ! outflow from routing store (QR) [mm/day]
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f
index 0c331a03..c174ed68 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f
@@ -50,14 +50,14 @@
StUH2=0.
!initialization of model states using StateStart
- St(1) = StateStart(1)
- St(2) = StateStart(2)
- St(3) = StateStart(3)
+ St(1) = StateStart(1)
+ St(2) = StateStart(2)
+ St(3) = StateStart(3)
DO I=1,NH
- StUH1(I)=StateStart(7+I)
+ StUH1(I)=StateStart(7+I)
ENDDO
DO I=1,2*NH
- StUH2(I)=StateStart(7+I+NH)
+ StUH2(I)=StateStart(7+I+NH)
ENDDO
!parameter values
@@ -69,13 +69,13 @@
!Param(6) : time constant of exponential store (X6 - EXP) [day]
!computation of HU ordinates
- OrdUH1 = 0.
- OrdUH2 = 0.
-
+ OrdUH1 = 0.
+ OrdUH2 = 0.
+
D=2.5
CALL UH1(OrdUH1,Param(4),D)
CALL UH2(OrdUH2,Param(4),D)
-
+
!initialization of model outputs
Q = -999.999
MISC = -999.999
@@ -96,18 +96,18 @@ c MISC = -999.999
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))
+ Outputs(k,I)=MISC(IndOutputs(I))
ENDDO
ENDDO
!model states at the end of the run
- StateEnd(1) = St(1)
- StateEnd(2) = St(2)
- StateEnd(3) = St(3)
+ StateEnd(1) = St(1)
+ StateEnd(2) = St(2)
+ StateEnd(3) = St(3)
DO K=1,NH
- StateEnd(7+K)=StUH1(K)
+ StateEnd(7+K)=StUH1(K)
ENDDO
DO K=1,2*NH
- StateEnd(7+NH+K)=StUH2(K)
+ StateEnd(7+NH+K)=StUH2(K)
ENDDO
RETURN
@@ -166,50 +166,50 @@ C**********************************************************************
C Production store
IF(P1.LE.E) THEN
- EN=E-P1
- PN=0.
- WS=EN/A
- IF(WS.GT.13)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- Sr = St(1)/A
- ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
- ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
- ! fin speed-up
-
- AE=ER+P1
- St(1)=St(1)-ER
- PS=0.
- PR=0.
+ EN=E-P1
+ PN=0.
+ WS=EN/A
+ IF(WS.GT.13)WS=13.
+
+ ! speed-up
+ TWS = tanHyp(WS)
+ Sr = St(1)/A
+ ER=St(1)*(2.-Sr)*TWS/(1.+(1.-Sr)*TWS)
+ ! ER=X(2)*(2.-X(2)/A)*tanHyp(WS)/(1.+(1.-X(2)/A)*tanHyp(WS))
+ ! fin speed-up
+
+ AE=ER+P1
+ St(1)=St(1)-ER
+ PS=0.
+ PR=0.
ELSE
- EN=0.
- AE=E
- PN=P1-E
- WS=PN/A
- IF(WS.GT.13)WS=13.
-
- ! speed-up
- TWS = tanHyp(WS)
- Sr = St(1)/A
- PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
- ! PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
- ! fin speed-up
-
- PR=PN-PS
- St(1)=St(1)+PS
+ EN=0.
+ AE=E
+ PN=P1-E
+ WS=PN/A
+ IF(WS.GT.13)WS=13.
+
+ ! speed-up
+ TWS = tanHyp(WS)
+ Sr = St(1)/A
+ PS=A*(1.-Sr*Sr)*TWS/(1.+Sr*TWS)
+ ! PS=A*(1.-(X(2)/A)**2.)*tanHyp(WS)/(1.+X(2)/A*tanHyp(WS))
+ ! fin speed-up
+
+ PR=PN-PS
+ St(1)=St(1)+PS
ENDIF
C Percolation from production store
IF(St(1).LT.0.)St(1)=0.
- ! speed-up
- ! (9/4)**4 = 25.62890625
- Sr = St(1)/Param(1)
- Sr = Sr * Sr
- Sr = Sr * Sr
+ ! speed-up
+ ! (9/4)**4 = 25.62890625
+ Sr = St(1)/Param(1)
+ Sr = Sr * Sr
+ Sr = Sr * Sr
PERC=St(1)*(1.-1./SQRT(SQRT(1.+Sr/25.62890625)))
! PERC=X(2)*(1.-(1.+(X(2)/(9./4.*Param(1)))**4.)**(-0.25))
- ! fin speed-up
+ ! fin speed-up
St(1)=St(1)-PERC
@@ -239,15 +239,15 @@ C Routing store
IF((St(2)+(1-C)*StUH1(1)+EXCH).LT.0) AEXCH1=-St(2)-(1-C)*StUH1(1)
St(2)=St(2)+(1-C)*StUH1(1)+EXCH
IF(St(2).LT.0.)St(2)=0.
-
- ! speed-up
- Rr = St(2)/Param(3)
- Rr = Rr * Rr
- Rr = Rr * Rr
+
+ ! 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
-
+ ! fin speed-up
+
St(2)=St(2)-QR
C Update of exponential store
@@ -293,7 +293,7 @@ C Variables storage
MISC(10)=StUH2(1) ! Q1 ! outflow from UH2 (Q1) [mm/day]
MISC(11)=St(2) ! Rout ! routing store level (St(2)) [mm]
MISC(12)=EXCH ! Exch ! potential third-exchange between catchments (EXCH) [mm/day]
- MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from routing store (AEXCH1) [mm/day]
+ MISC(13)=AEXCH1 ! AExch1 ! actual exchange between catchments from routing store (AEXCH1) [mm/day]
MISC(14)=AEXCH2 ! AExch2 ! actual exchange between catchments from direct branch (after UH2) (AEXCH2) [mm/day]
MISC(15)=AEXCH1+AEXCH2+EXCH ! AExch ! actual total exchange between catchments (AEXCH1+AEXCH2+EXCH) [mm/day]
MISC(16)=QR ! QR ! outflow from routing store (QR) [mm/day]
--
GitLab