-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathetact.f
308 lines (279 loc) · 12 KB
/
etact.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
subroutine etact
!! ~ ~ ~ PURPOSE ~ ~ ~
!! this subroutine calculates potential plant transpiration for Priestley-
!! Taylor and Hargreaves ET methods, and potential and actual soil
!! evaporation. NO3 movement into surface soil layer due to evaporation
!! is also calculated.
!! ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! canstor(:) |mm H2O |amount of water held in canopy storage
!! elevb(:,:) |m |elevation at center of band in subbasin
!! elevb_fr(:,:)|none |fraction of subbasin area within elevation
!! |band
!! ep_max |mm H2O |maximum amount of transpiration (plant et)
!! |that can occur on current day in HRU
!! esco(:) |none |soil evaporation compensation factor
!! ihru |none |HRU number
!! ipet |none |code for potential ET method
!! |0 Priestley-Taylor method
!! |1 Penman/Monteith method
!! |2 Hargreaves method
!! laiday(:) |m**2/m**2 |leaf area index
!! pet_day |mm H2O |potential evapotranspiration on current day
!! |in HRU
!! pot_vol(:) |m**3 H2O |current volume of water stored in the
!! |depression/impounded area
!! sno_hru(:) |mm H2O |amount of water in snow in HRU on current day
!! snoeb(:,:) |mm H2O |snow water content in elevation band on
!! |current day
!! sol_cov(:) |kg/ha |amount of residue on soil surface
!! sol_fc(:,:) |mm H2O |amount of water available to plants in soil
!! |layer at field capacity (fc - wp water)
!! sol_nly(:) |none |number of soil layers in profile
!! sol_no3(:,:) |kg N/ha |amount of nitrogen stored in the nitrate
!! |pool
!! sol_st(:,:) |mm H2O |amount of water stored in the soil layer on
!! |current day
!! sol_z(:,:) |mm |depth to bottom of soil layer
!! tavband(:,:) |deg C |average temperature for the day in band in HRU
!! tmpav(:) |deg C |average air temperature on current day for
!! |HRU
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! canev |mm H2O |amount of water evaporated from canopy
!! |storage
!! ep_max |mm H2O |maximum amount of transpiration (plant et)
!! |that can occur on current day in HRU
!! es_day |mm H2O |actual amount of evaporation (soil et) that
!! |occurs on day in HRU
!! sno_hru(:) |mm H2O |amount of water in snow in HRU on current day
!! sno3up |kg N/ha |amount of nitrate moving upward in the soil
!! |profile in watershed
!! snoeb(:,:) |mm H2O |snow water content in elevation band on
!! |current day
!! snoev |mm H2O |amount of water in snow lost through
!! |sublimation on current day
!! sol_st(:,:) |mm H2O |amount of water stored in the soil layer on
!! |current day
!! sol_sw(:) |mm H2O |amount of water stored in the soil profile
!! |on current day
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! cej |
!! dep |mm |soil depth from which evaporation will occur
!! |in current soil layer
!! eaj |none |weighting factor to adjust PET for impact of
!! |plant cover
!! effnup |
!! eos1 |none |variable to hold intermediate calculation
!! |result
!! eosl |mm H2O |maximum amount of evaporation that can occur
!! |from soil profile
!! es_max |mm H2O |maximum amount of evaporation (soil et)
!! |that can occur on current day in HRU
!! esd |mm |maximum soil depth from which evaporation
!! |is allowed to occur
!! esleft |mm H2O |potenial soil evap that is still available
!! etco |
!! evz |
!! evzp |
!! ib |none |counter
!! j |none |HRU number
!! ly |none |counter
!! no3up |kg N/ha |amount of nitrate moving upward in profile
!! pet |mm H2O |amount of PET remaining after water stored
!! |in canopy is evaporated
!! sev |mm H2O |amount of evaporation from soil layer
!! sumsnoeb |mm H2O |amount of snow in elevation bands whose air
!! |temperature is greater than 0 degrees C
!! xx |none |variable to hold intermediate calculation
!! |result
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!! Intrinsic: Exp, Min, Max
!! SWAT: Expo
!! ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~
use parm
integer :: j, ib, ly
!! real, parameter :: esd = 500., etco = 0.80, effnup = 0.1
real :: esd, etco, effnup
real :: no3up, es_max, eos1, xx, cej, eaj, pet, esleft
real :: sumsnoeb, evzp, eosl, dep, evz, sev
j = 0
j = ihru
pet = 0.
pet = pet_day
!! added statements for test of real statement above
esd = 500.
etco = 0.80
effnup = 0.1
!! evaporate canopy storage first
!! canopy storage is calculated by the model only if the Green & Ampt
!! method is used to calculate surface runoff. The curve number methods
!! take canopy effects into account in the equations. For either of the
!! CN methods, canstor will always equal zero.
pet = pet - canstor(j)
if (pet < 0.) then
canstor(j) = -pet
canev = pet_day
pet = 0.
ep_max = 0.
es_max = 0.
else
canev = canstor(j)
canstor(j) = 0.
endif
if (pet > 1.0e-6) then
!! compute potential plant evap for methods other that Penman-Monteith
if (ipet /= 1) then
if (laiday(j) <= 3.0) then
ep_max = laiday(j) * pet / 3.
else
ep_max = pet
end if
if (ep_max < 0.) ep_max = 0.
end if
!! compute potential soil evaporation
cej = -5.e-5
eaj = 0.
es_max = 0.
eos1 = 0.
if (sno_hru(j) >= 0.5) then
eaj = 0.5
else
eaj = Exp(cej * (sol_cov(j)+ 0.1))
end if
es_max = pet * eaj
eos1 = pet / (es_max + ep_max + 1.e-10)
eos1 = es_max * eos1
es_max = Min(es_max, eos1)
es_max = Max(es_max, 0.)
! if (pot_vol(j) > 1.e-4) es_max = 0.
!! make sure maximum plant and soil ET doesn't exceed potential ET
!!if (pet_day < es_max + ep_max) then
!!es_max = pet_day - ep_max
if (pet < es_max + ep_max) then
es_max = pet * es_max / (es_max + ep_max)
ep_max = pet * ep_max / (es_max + ep_max)
end if
if (pet < es_max + ep_max) then
es_max = pet - ep_max - 1.0e-6
end if
!!end if
!! initialize soil evaporation variables
esleft = 0.
esleft = es_max
!! compute sublimation
if (elevb_fr(1,hru_sub(j)) <= 0.) then
!! compute sublimation without elevation bands
if (tmpav(j) > 0.) then
if (sno_hru(j) >= esleft) then
!! take all soil evap from snow cover
sno_hru(j) = sno_hru(j) - esleft
snoev = snoev + esleft
esleft = 0.
else
!! take all soil evap from snow cover then start taking from soil
esleft = esleft - sno_hru(j)
snoev = snoev + sno_hru(j)
sno_hru(j) = 0.
endif
endif
else
!! elevation bands
sumsnoeb = 0.
!! calculate air temp in elevation bands and sum snow
!! for elevation bands with temp > 0 deg C
do ib = 1, 10
if (elevb_fr(ib,hru_sub(j)) <= 0.) exit
if (tavband(ib,j) > 0.) then
sumsnoeb = sumsnoeb +
& snoeb(ib,j) * elevb_fr(ib,hru_sub(j))
end if
end do
!! compute sublimation from elevation bands
if (sumsnoeb >= esleft .and. sumsnoeb > 0.01) then
do ib = 1, 10
if (elevb_fr(ib,hru_sub(j)) <= 0.) exit
if (tavband(ib,j) > 0.) then
snoev = snoev + snoeb(ib,j) * (esleft / sumsnoeb) *
& elevb_fr(ib,hru_sub(j))
snoeb(ib,j) = snoeb(ib,j) - snoeb(ib,j) * (esleft /
& sumsnoeb)
end if
end do
else
do ib = 1, 10
if (elevb_fr(ib,hru_sub(j)) <= 0.) exit
if (tavband(ib,j) > 0.) then
snoev = snoev + snoeb(ib,j) * elevb_fr(ib,hru_sub(j))
snoeb(ib,j) = 0.
end if
end do
end if
esleft = esleft - snoev
sno_hru(j) = sno_hru(j) - snoev
endif
!! take soil evap from each soil layer
evzp = 0.
eosl = 0.
eosl = esleft
do ly = 1, sol_nly(j)
!! depth exceeds max depth for soil evap (esd)
dep = 0.
if (ly == 1) then
dep = sol_z(1,j)
else
dep = sol_z(ly-1,j)
endif
if (dep < esd) then
!! calculate evaporation from soil layer
evz = 0.
sev = 0.
xx = 0.
evz = eosl * sol_z(ly,j) / (sol_z(ly,j) + Exp(2.374 -
& .00713 * sol_z(ly,j)))
sev = evz - evzp * esco(j)
evzp = evz
if (sol_st(ly,j) < sol_fc(ly,j)) then
xx = 2.5 * (sol_st(ly,j) - sol_fc(ly,j)) / sol_fc(ly,j)
sev = sev * Expo(xx)
end if
sev = Min(sev, sol_st(ly,j) * etco)
if (sev < 0.) sev = 0.
if (sev > esleft) sev = esleft
!! adjust soil storage, potential evap
if (sol_st(ly,j) > sev) then
esleft = esleft - sev
sol_st(ly,j) = Max(1.e-6, sol_st(ly,j) - sev)
else
esleft = esleft - sol_st(ly,j)
sol_st(ly,j) = 0.
endif
endif
!! compute no3 flux from layer 2 to 1 by soil evaporation
if (ly == 2) then
no3up = 0.
no3up = effnup * sev * sol_no3(2,j) / (sol_st(2,j) + 1.e-6)
no3up = Min(no3up, sol_no3(2,j))
sno3up = sno3up + no3up * hru_dafr(j)
sol_no3(2,j) = sol_no3(2,j) - no3up
sol_no3(1,j) = sol_no3(1,j) + no3up
endif
end do
!! update total soil water content
sol_sw(j) = 0.
do ly = 1, sol_nly(j)
sol_sw(j) = sol_sw(j) + sol_st(ly,j)
end do
!! calculate actual amount of evaporation from soil
es_day = es_max - esleft
if (es_day < 0.) es_day = 0.
end if
return
end