frun_GR2M.f90 8.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
!------------------------------------------------------------------------------
!    Subroutines relative to the annual GR2M model
!------------------------------------------------------------------------------
! TITLE   : airGR
! PROJECT : airGR
! FILE    : frun_GR2M.f
!------------------------------------------------------------------------------
! AUTHORS
! Original code: S. Mouelhi
! Cleaning and formatting for airGR: L. Coron
! Further cleaning: G. Thirel
!------------------------------------------------------------------------------
! Creation date: 2003
14
! Last modified: 16/04/2020
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
!------------------------------------------------------------------------------
! REFERENCES
! Mouelhi S. (2003). Vers une chaîne cohérente de modèles pluie-débit 
! conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et 
! journalier. PhD thesis (in French), ENGREF, Cemagref Antony, France.
!
! Mouelhi, S., C. Michel, C. Perrin and V. Andréassian (2006). Stepwise 
! development of a two-parameter monthly water balance model. Journal of 
! Hydrology, 318(1-4), 200-214. doi:10.1016/j.jhydrol.2005.06.014.
!------------------------------------------------------------------------------
! Quick description of public procedures:
!         1. frun_gr2m
!         2. MOD_GR2M
!------------------------------------------------------------------------------


31
32
33
      SUBROUTINE frun_gr2m(LInputs,InputsPrecip,InputsPE,NParam,Param, &
                           NStates,StateStart,NOutputs,IndOutputs, &
                           Outputs,StateEnd)
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
! Subroutine that initializes GR2M, get its parameters, performs the call 
! to the MOD_GR2M subroutine at each time step, and stores the final states
! Inputs
!       LInputs      ! Integer, length of input and output series
!       InputsPrecip ! Vector of real, input series of total precipitation [mm/month]
!       InputsPE     ! Vector of real, input series of potential evapotranspiration (PE) [mm/month]
!       NParam       ! Integer, number of model parameters
!       Param        ! Vector of real, parameter set
!       NStates      ! Integer, number of state variables
!       StateStart   ! Vector of real, state variables used when the model run starts (store levels [mm])
!       NOutputs     ! Integer, number of output series
!       IndOutputs   ! Vector of integer, indices of output series
! Outputs      
!       Outputs      ! Vector of real, output series
!       StateEnd     ! Vector of real, state variables at the end of the model run (store levels [mm])
Delaigue Olivier's avatar
Delaigue Olivier committed
49
50


51
      !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr2m
Delaigue Olivier's avatar
Delaigue Olivier committed
52
53
54


      Implicit None
55
56
57

      !! dummies
      ! in
Delaigue Olivier's avatar
Delaigue Olivier committed
58
      integer, intent(in) :: LInputs,NParam,NStates,NOutputs
59
60
61
62
63
64
65
      doubleprecision, dimension(LInputs), intent(in) :: InputsPrecip
      doubleprecision, dimension(LInputs), intent(in) :: InputsPE
      doubleprecision, dimension(NParam),  intent(in) :: Param
      doubleprecision, dimension(NStates), intent(in) :: StateStart
      integer, dimension(NOutputs),        intent(in) :: IndOutputs
      ! out
      doubleprecision, dimension(NStates), intent(out) :: StateEnd
66
      doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
67
68
69
70
71
72
73
      
      !! locals
      integer :: I,K
      integer, parameter :: NMISC=30
      doubleprecision, dimension(2) :: St
      doubleprecision, dimension(NMISC) :: MISC
      doubleprecision :: P,E,Q
Delaigue Olivier's avatar
Delaigue Olivier committed
74
75

      !--------------------------------------------------------------
76
      ! Initializations
Delaigue Olivier's avatar
Delaigue Olivier committed
77
78
      !--------------------------------------------------------------

79
      ! initialization of model states to zero
80
      St=0.
Delaigue Olivier's avatar
Delaigue Olivier committed
81

82
      ! initialization of model states using StateStart
83
84
      St(1)=StateStart(1)
      St(2)=StateStart(2)
Delaigue Olivier's avatar
Delaigue Olivier committed
85

86
87
88
      ! parameter values
      ! Param(1) : production store capacity [mm]
      ! Param(2) : groundwater exchange coefficient [-]
Delaigue Olivier's avatar
Delaigue Olivier committed
89

90
      ! initialization of model outputs
Delaigue Olivier's avatar
Delaigue Olivier committed
91
92
      Q = -999.999
      MISC = -999.999
93
94
!      StateEnd = -999.999 ! initialization made in R
!      Outputs = -999.999  ! initialization made in R
Delaigue Olivier's avatar
Delaigue Olivier committed
95
96
97
98



      !--------------------------------------------------------------
99
      ! Time loop
Delaigue Olivier's avatar
Delaigue Olivier committed
100
101
102
103
      !--------------------------------------------------------------
      DO k=1,LInputs
        P =InputsPrecip(k)
        E =InputsPE(k)
104
105
!        Q = -999.999
!        MISC = -999.999
106
107
        !model run on one time step
        CALL MOD_GR2M(St,Param,P,E,Q,MISC)
Delaigue Olivier's avatar
Delaigue Olivier committed
108
109
110
111
112
        !storage of outputs
        DO I=1,NOutputs
        Outputs(k,I)=MISC(IndOutputs(I))
        ENDDO
      ENDDO
113
      ! model states at the end of the run
114
115
      StateEnd(1)=St(1)
      StateEnd(2)=St(2)
Delaigue Olivier's avatar
Delaigue Olivier committed
116
117
118
119
120
121
122
123
124

      RETURN

      ENDSUBROUTINE





125
!################################################################################################################################
Delaigue Olivier's avatar
Delaigue Olivier committed
126
127
128
129




130
!**********************************************************************
131
      SUBROUTINE MOD_GR2M(St,Param,P,E,Q,MISC)
132
133
! Calculation of streamflow on a single time step (month) with the GR2M model
! Inputs:
134
135
136
137
!       St     Vector of real, model states at the beginning of the time step [mm]
!       Param  Vector of real, model parameters (Param(1) [mm]; Param(2) [-])
!       P      Real, value of rainfall during the time step [mm/month]
!       E      Real, value of potential evapotranspiration during the time step [mm/month]
138
! Outputs:
139
140
141
!       St     Vector of real, model states at the end of the time step [mm]
!       Q      Real, value of simulated flow at the catchment outlet for the time step [mm/month]
!       MISC   Vector of real, model outputs for the time step [mm/month]
142
!**********************************************************************
Delaigue Olivier's avatar
Delaigue Olivier committed
143
      Implicit None
144

145
146
      !! locals
      integer, parameter :: NParam=2,NMISC=30
147
      doubleprecision :: WS,S1,S2
148
      doubleprecision :: P1,P2,P3,R1,R2,AE,EXCH,PS
149
      doubleprecision :: expWS, TWS, Sr ! speed-up
150
151
152
153
154
155
156
157
158
159
160

      !! dummies
      ! in
      doubleprecision, dimension(NParam), intent(in) :: Param
      doubleprecision, intent(in) :: P,E
      ! inout
      doubleprecision, dimension(2), intent(inout) :: St
      ! out
      doubleprecision, intent(out) :: Q
      doubleprecision, dimension(NMISC), intent(out) :: MISC
      
161
! Production store
Delaigue Olivier's avatar
Delaigue Olivier committed
162
      WS=P/Param(1)  
163
      IF(WS.GT.13.) WS=13.
164
165

      ! speed-up
166
167
      expWS = exp(2.*WS)
      TWS = (expWS - 1.)/(expWS + 1.)
168
169
      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))                 
170
      ! end speed-up
171

172
173
      P1=P+St(1)-S1
      PS = P - P1      
Delaigue Olivier's avatar
Delaigue Olivier committed
174
      WS=E/Param(1)         
175
      IF(WS.GT.13.) WS=13.
176
177

      ! speed-up
178
179
      expWS = exp(2.*WS)
      TWS = (expWS - 1.)/(expWS + 1.)
180
181
      S2=S1*(1.-TWS)/(1.+(1.-S1/Param(1))*TWS)  
      ! S2=S1*(1.-tanHyp(WS))/(1.+(1.-S1/Param(1))*tanHyp(WS))  
182
      ! end speed-up
183
      AE = S1 - S2
Delaigue Olivier's avatar
Delaigue Olivier committed
184
                
185
! Percolation
186
187
188
      ! speed-up
      Sr = S2/Param(1)
      Sr = Sr * Sr * Sr + 1.
189
      St(1)=S2/Sr**(1./3.)
190
      ! X(1)=S2/(1+(S2/Param(1))**3.)**(1./3.)         
191
      ! end speed-up
192

193
      P2=S2-St(1)  
Delaigue Olivier's avatar
Delaigue Olivier committed
194
195
      P3=P1+P2

196
! QR calculation (routing store)
197
      R1=St(2)+P3
Delaigue Olivier's avatar
Delaigue Olivier committed
198

199
! Water exchange
200
      R2=Param(2)*R1
201
      EXCH = R2 - R1
Delaigue Olivier's avatar
Delaigue Olivier committed
202

203
! Total runoff
Delaigue Olivier's avatar
Delaigue Olivier committed
204
205
      Q=R2*R2/(R2+60.)

206
! Updating store level
207
      St(2)=R2-Q
Delaigue Olivier's avatar
Delaigue Olivier committed
208
209


210
! Variables storage
211
      MISC( 1)=E             ! PE     ! [numeric] observed potential evapotranspiration [mm/month]
212
      MISC( 2)=P             ! Precip ! [numeric] observed total precipitation  [mm/month]
213
      MISC( 3)=St(1)         ! Prod   ! [numeric] production store level (St(1)) [mm]
214
      MISC( 4)=P1            ! Pn     ! [numeric] net rainfall (P1) [mm/month]
215
216
217
218
219
220
221
      MISC( 5)=PS            ! Ps     ! [numeric] part of P filling the production store [mm/month]
      MISC( 6)=AE            ! AE     ! [numeric] actual evapotranspiration [mm/month]
      MISC( 7)=P2            ! Perc   ! [numeric] percolation (P2) [mm/month]
      MISC( 8)=P3            ! PR     ! [numeric] P3=P1+P2 [mm/month]
      MISC( 9)=St(2)         ! Rout   ! [numeric] routing store level (St(2)) [mm]
      MISC(10)=EXCH          ! EXCH   ! [numeric] groundwater exchange (EXCH) [mm/month]
      MISC(11)=Q             ! Qsim   ! [numeric] simulated outflow at catchment outlet [mm/month]
Delaigue Olivier's avatar
Delaigue Olivier committed
222
223


224
      END SUBROUTINE
Delaigue Olivier's avatar
Delaigue Olivier committed
225