frun_GR2M.f 6.25 KB
Newer Older
Delaigue Olivier's avatar
Delaigue Olivier committed
1
2
3
4


      SUBROUTINE frun_GR2M(
                                 !inputs
5
6
7
8
9
10
11
12
13
     &                             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
Delaigue Olivier's avatar
Delaigue Olivier committed
14
                                 !outputs
15
16
     &                             Outputs      , ! [double]  output series
     &                             StateEnd     ) ! [double]  state variables at the end of the model run (store levels [mm])
Delaigue Olivier's avatar
Delaigue Olivier committed
17
18


19
      !DEC$ ATTRIBUTES DLLEXPORT :: frun_GR2M
Delaigue Olivier's avatar
Delaigue Olivier committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35


      Implicit None
      !### input and output variables
      integer, intent(in) :: LInputs,NParam,NStates,NOutputs
      doubleprecision, dimension(LInputs)  :: InputsPrecip
      doubleprecision, dimension(LInputs)  :: InputsPE
      doubleprecision, dimension(NParam)   :: Param
      doubleprecision, dimension(NStates)  :: StateStart
      doubleprecision, dimension(NStates)  :: StateEnd
      integer, dimension(NOutputs) :: IndOutputs
      doubleprecision, dimension(LInputs,NOutputs) :: Outputs

      !parameters, internal states and variables
      integer NMISC
      parameter (NMISC=30)
36
      doubleprecision St(2)
Delaigue Olivier's avatar
Delaigue Olivier committed
37
38
39
40
41
      doubleprecision MISC(NMISC)
      doubleprecision P,E,Q
      integer I,K

      !--------------------------------------------------------------
42
      !Initializations
Delaigue Olivier's avatar
Delaigue Olivier committed
43
44
      !--------------------------------------------------------------

45
46
      !initialization of model states to zero
      St=0.
Delaigue Olivier's avatar
Delaigue Olivier committed
47
48

      !initilisation of model states using StateStart
49
50
      St(1)=StateStart(1)
      St(2)=StateStart(2)
Delaigue Olivier's avatar
Delaigue Olivier committed
51
52
53

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

56
      !initialization of model outputs
Delaigue Olivier's avatar
Delaigue Olivier committed
57
58
      Q = -999.999
      MISC = -999.999
59
60
c      StateEnd = -999.999 !initialization made in R
c      Outputs = -999.999  !initialization made in R
Delaigue Olivier's avatar
Delaigue Olivier committed
61
62
63
64
65
66
67
68
69
70
71



      !--------------------------------------------------------------
      !Time loop
      !--------------------------------------------------------------
      DO k=1,LInputs
        P =InputsPrecip(k)
        E =InputsPE(k)
c        Q = -999.999
c        MISC = -999.999
72
73
        !model run on one time step
        CALL MOD_GR2M(St,Param,P,E,Q,MISC)
Delaigue Olivier's avatar
Delaigue Olivier committed
74
75
76
77
78
79
        !storage of outputs
        DO I=1,NOutputs
        Outputs(k,I)=MISC(IndOutputs(I))
        ENDDO
      ENDDO
      !model states at the end of the run
80
81
      StateEnd(1)=St(1)
      StateEnd(2)=St(2)
Delaigue Olivier's avatar
Delaigue Olivier committed
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

      RETURN

      ENDSUBROUTINE





c################################################################################################################################




C**********************************************************************
97
98
      SUBROUTINE MOD_GR2M(St,Param,P,E,Q,MISC)
C Calculation of streamflow on a single time step (month) with the GR2M model
Delaigue Olivier's avatar
Delaigue Olivier committed
99
C Inputs:
100
101
102
103
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]
Delaigue Olivier's avatar
Delaigue Olivier committed
104
C Outputs:
105
106
107
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]
Delaigue Olivier's avatar
Delaigue Olivier committed
108
109
110
111
112
C**********************************************************************
      Implicit None
      INTEGER NMISC,NParam
      PARAMETER (NMISC=30)
      PARAMETER (NParam=2)
113
      DOUBLEPRECISION St(2)
Delaigue Olivier's avatar
Delaigue Olivier committed
114
115
116
117
      DOUBLEPRECISION Param(NParam)
      DOUBLEPRECISION MISC(NMISC)
      DOUBLEPRECISION P,E,Q
      DOUBLEPRECISION WS,tanHyp,S1,S2
118
      DOUBLEPRECISION P1,P2,P3,R1,R2,AE,EXCH
Delaigue Olivier's avatar
Delaigue Olivier committed
119

120
121
      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
	  
Delaigue Olivier's avatar
Delaigue Olivier committed
122
123
C Production store
      WS=P/Param(1)  
124
      IF(WS.GT.13.)WS=13.
125
126
127
	  
 	  ! speed-up
	  TWS = tanHyp(WS)
128
	  S1=(St(1)+Param(1)*TWS)/(1.+St(1)/Param(1)*TWS)                 
129
130
131
	  ! S1=(X(1)+Param(1)*tanHyp(WS))/(1.+X(1)/Param(1)*tanHyp(WS))                 
  	  ! fin speed-up

132
	  P1=P+St(1)-S1                  
Delaigue Olivier's avatar
Delaigue Olivier committed
133
      WS=E/Param(1)         
134
      IF(WS.GT.13.)WS=13.
135
136
137
138
139
140
	  
 	  ! 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
141
	  AE = S1 - S2
Delaigue Olivier's avatar
Delaigue Olivier committed
142
143
                
C Percolation
144
145
146
 	  ! speed-up
	  Sr = S2/Param(1)
	  Sr = Sr * Sr * Sr + 1.
147
      St(1)=S2/Sr**(1./3.)
148
149
150
      ! X(1)=S2/(1+(S2/Param(1))**3.)**(1./3.)         
 	  ! fin speed-up

151
      P2=S2-St(1)  
Delaigue Olivier's avatar
Delaigue Olivier committed
152
153
154
      P3=P1+P2

C QR calculation (routing store)
155
      R1=St(2)+P3
Delaigue Olivier's avatar
Delaigue Olivier committed
156
157

C Water exchange
158
159
      R2=Param(2)*R1
	  EXCH = R2 - R1
Delaigue Olivier's avatar
Delaigue Olivier committed
160
161
162
163
164

C Total runoff
      Q=R2*R2/(R2+60.)

C Updating store level
165
      St(2)=R2-Q
Delaigue Olivier's avatar
Delaigue Olivier committed
166
167
168


C Variables storage
169
      MISC( 1)=E             ! PE     ! [numeric] observed potential evapotranspiration [mm/month]
170
      MISC( 2)=P             ! Precip ! [numeric] observed total precipitation  [mm/month]
171
      MISC( 3)=AE            ! AE     ! [numeric] actual evapotranspiration [mm/month]
172
173
174
      MISC( 4)=P1            ! Pn     ! [numeric] net rainfall (P1) [mm/month]
      MISC( 5)=P2            ! Perc   ! [numeric] percolation (P2) [mm/month]
      MISC( 6)=P3            ! PR     ! [numeric] P3=P1+P2 [mm/month]
175
176
177
178
      MISC( 7)=EXCH          ! EXCH   ! [numeric] groundwater exchange (EXCH) [mm/month]
      MISC( 8)=St(1)         ! Prod   ! [numeric] production store level (St(1)) [mm]
      MISC( 9)=St(2)         ! Rout   ! [numeric] routing store level (St(2)) [mm]
      MISC(10)=Q             ! Qsim   ! [numeric] simulated outflow at catchment outlet [mm/month]
Delaigue Olivier's avatar
Delaigue Olivier committed
179
180
181
182
183


      ENDSUBROUTINE