-
Notifications
You must be signed in to change notification settings - Fork 0
/
etpot.f
327 lines (274 loc) · 12.3 KB
/
etpot.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
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
subroutine etpot
!! ~ ~ ~ PURPOSE ~ ~ ~
!! this subroutine calculates potential evapotranspiration using one
!! of three methods. If Penman-Monteith is being used, potential plant
!! transpiration is also calculated.
!! ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! laiday(:) |m**2/m**2 |leaf area index
!! albday |none |albedo for the day in HRU
!! cht(:) |m |canopy height
!! co2(:) |ppmv |CO2 concentration
!! gsi(:) |m/s |maximum stomatal conductance
!! hru_ra(:) |MJ/m^2 |solar radiation for the day in HRU
!! hru_rmx(:) |MJ/m^2 |maximum possible radiation for the day in HRU
!! hru_sub(:) |none |subbasin in which HRU is located
!! icr(:) |none |sequence number of crop grown within the
!! |current year
!! idplt(:) |none |land cover code from crop.dat
!! igro(:) |none |land cover status code
!! |0 no land cover currently growing
!! |1 land cover growing
!! ihru |none |HRU number
!! ipet |none |code for potential ET method
!! |0 Priestley-Taylor method
!! |1 Penman/Monteith method
!! |2 Hargreaves method
!! |3 read in PET values
!! nro(:) |none |sequence number of year in rotation
!! petmeas |mm H2O |potential ET value read in for day
!! rhd(:) |none |relative humidity for the day in HRU
!! sno_hru(:) |mm H2O |amount of water in snow in HRU on current day
!! sub_elev(:)|m |elevation of HRU
!! tmn(:) |deg C |minimum air temperature on current day in HRU
!! tmpav(:) |deg C |average air temperature on current day in HRU
!! tmx(:) |deg C |maximum air temperature on current day for HRU
!! u10(:) |m/s |wind speed (measured at 10 meters above
!! |surface)
!! vpd2(:) |(m/s)*(1/kPa) |rate of decline in stomatal conductance per
!! | |unit increase in vapor pressure deficit
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ep_max |mm H2O |maximum amount of transpiration (plant et)
!! |that can occur on current day in HRU
!! pet_day |mm H2O |potential evapotranspiration on current day in
!! |HRU
!! vpd |kPa |vapor pressure deficit
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! chz |cm |vegetation height
!! d |cm |displacement height for plant type
!! dlt |kPa/deg C |slope of the saturation vapor pressure-
!! |temperature curve
!! ea |kPa |saturated vapor pressure
!! ed |kPa |actual vapor pressure
!! fvpd |kPa |amount of vapro pressure deficit over
!! |threshold value
!! gma |kPa/deg C |psychrometric constant
!! j |none |HRU number
!! pb |kPa |mean atmospheric pressure
!! pet_alpha |none |alpha factor in Priestley-Taylor PET equation
!! ralb |MJ/m2 |net incoming radiation for PET
!! ralb1 |MJ/m2 |net incoming radiation
!! ramm |MJ/m2 |extraterrestrial radiation
!! rbo |none |net emissivity
!! rc |s/m |canopy resistance
!! rho |MJ/(m3*kPa) |K1*0.622*xl*rho/pb
!! rn |MJ/m2 |net radiation
!! rn_pet |MJ/m2 |net radiation for continuous crop cover
!! rout |MJ/m2 |outgoing radiation
!! rto |none |cloud cover factor
!! rv |s/m |aerodynamic resistance to sensible heat and
!! |vapor transfer
!! tk |deg K |average air temperature on current day for HRU
!! uzz |m/s |wind speed at height zz
!! xl |MJ/kg |latent heat of vaporization
!! xx |kPa |difference between vpd and vpthreshold
!! zom |cm |roughness length for momentum transfer
!! zov |cm |roughness length for vapor transfer
!! zz |cm |height at which wind speed is determined
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!! Intrinsic: Log, Sqrt, Max, Min
!! SWAT: Ee
!! ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~
use parm
integer :: j
real :: tk, pb, gma, xl, ea, ed, dlt, ramm, ralb1, ralb, xx
real :: rbo, rto, rn, uzz, zz, zom, zov, rv, rn_pet, fvpd
real :: rc, rho, rout, d, chz, gsi_adj, pet_alpha
!! initialize local variables
j = 0
j = ihru
tk = 0.
tk = tmpav(j) + 273.15
!! calculate mean barometric pressure
pb = 0.
pb = 101.3 - sub_elev(hru_sub(j)) *
& (0.01152 - 0.544e-6 * sub_elev(hru_sub(j)))
!! calculate latent heat of vaporization
xl = 0.
xl = 2.501 - 2.361e-3 * tmpav(j)
!! calculate psychrometric constant
gma = 0.
gma = 1.013e-3 * pb / (0.622 * xl)
!! calculate saturation vapor pressure, actual vapor pressure and
!! vapor pressure deficit
ea = 0.
ed = 0.
vpd = 0.
ea = Ee(tmpav(j))
ed = ea * rhd(j)
vpd = ea - ed
!!calculate the slope of the saturation vapor pressure curve
dlt = 0.
dlt = 4098. * ea / (tmpav(j) + 237.3)**2
!! DETERMINE POTENTIAL ET
select case (ipet)
case (0) !! PRIESTLEY-TAYLOR POTENTIAL EVAPOTRANSPIRATION METHOD
!! net radiation
!! calculate net short-wave radiation for PET
ralb = 0.
if (sno_hru(j) <= .5) then
ralb = hru_ra(j) * (1.0 - 0.23)
else
ralb = hru_ra(j) * (1.0 - 0.8)
end if
!! calculate net long-wave radiation
!! net emissivity equation 2.2.20 in SWAT manual
rbo = 0.
rbo = -(0.34 - 0.139 * Sqrt(ed))
!! cloud cover factor equation 2.2.19
rto = 0.
if (hru_rmx(j) < 1.e-4) then
rto = 0.
else
rto = 0.9 * (hru_ra(j) / hru_rmx(j)) + 0.1
end if
!! net long-wave radiation equation 2.2.21
rout = 0.
rout = rbo * rto * 4.9e-9 * (tk**4)
!! calculate net radiation
rn_pet = 0.
rn_pet = ralb + rout
!! net radiation
pet_alpha = 1.28
pet_day = pet_alpha * (dlt / (dlt + gma)) * rn_pet / xl
pet_day = Max(0., pet_day)
case (1) !! PENMAN-MONTEITH POTENTIAL EVAPOTRANSPIRATION METHOD
!! net radiation
!! calculate net short-wave radiation for PET
ralb = 0.
if (sno_hru(j) <= .5) then
ralb = hru_ra(j) * (1.0 - 0.23)
else
ralb = hru_ra(j) * (1.0 - 0.8)
end if
!! calculate net short-wave radiation for max plant ET
ralb1 = 0.
ralb1 = hru_ra(j) * (1.0 - albday)
!! calculate net long-wave radiation
!! net emissivity equation 2.2.20 in SWAT manual
rbo = 0.
rbo = -(0.34 - 0.139 * Sqrt(ed))
!! cloud cover factor equation 2.2.19
rto = 0.
if (hru_rmx(j) < 1.e-4) then
rto = 0.
else
rto = 0.9 * (hru_ra(j) / hru_rmx(j)) + 0.1
end if
!! net long-wave radiation equation 2.2.21
rout = 0.
rout = rbo * rto * 4.9e-9 * (tk**4)
!! calculate net radiation
rn = 0.
rn_pet = 0.
rn = ralb1 + rout
rn_pet = ralb + rout
!! net radiation
rho = 0.
rho = 1710. - 6.85 * tmpav(j)
if (u10(j) < 0.01) u10(j) = 0.01
!! potential ET: reference crop alfalfa at 40 cm height
rv = 0.
rc = 0.
rv = 114. / (u10(j) * (170./1000.)**0.2)
rc = 49. / (1.4 - 0.4 * co2(hru_sub(j)) / 330.)
pet_day = (dlt * rn_pet + gma * rho * vpd / rv) /
& (xl * (dlt + gma * (1. + rc / rv)))
pet_day = Max(0., pet_day)
!! maximum plant ET
if (igro(j) <= 0) then
ep_max = 0.0
else
!! determine wind speed and height of wind speed measurement
!! adjust to 100 cm (1 m) above canopy if necessary
uzz = 0.
zz = 0.
if (cht(j) <= 1.0) then
zz = 170.0
else
zz = cht(j) * 100. + 100.
end if
uzz = u10(j) * (zz/1000.)**0.2
!! calculate canopy height in cm
chz = 0.
if (cht(j) < 0.01) then
chz = 1.
else
chz = cht(j) * 100.
end if
!! calculate roughness length for momentum transfer
zom = 0.
if (chz <= 200.) then
zom = 0.123 * chz
else
zom = 0.058 * chz**1.19
end if
!! calculate roughness length for vapor transfer
zov = 0.
zov = 0.1 * zom
!! calculate zero-plane displacement of wind profile
d = 0.
d = 0.667 * chz
!! calculate aerodynamic resistance
rv = Log((zz - d) / zom) * Log((zz - d) / zov)
rv = rv / ((0.41)**2 * uzz)
!! adjust stomatal conductivity for low vapor pressure
!! this adjustment will lower maximum plant ET for plants
!! sensitive to very low vapor pressure
xx = 0.
fvpd = 0.
xx = vpd - 1.
if (xx > 0.0) then
fvpd = Max(0.1,1.0 - vpd2(idplt(j)) * xx)
else
fvpd = 1.0
end if
gsi_adj = 0.
gsi_adj = gsi(idplt(j)) * fvpd
!! calculate canopy resistance
rc = 0.
rc = 1. / gsi_adj !single leaf resistance
rc = rc / (0.5 * (laiday(j) + 0.01)
& * (1.4 - 0.4 * co2(hru_sub(j)) / 330.))
!! calculate maximum plant ET
ep_max = (dlt * rn + gma * rho * vpd / rv) /
& (xl * (dlt + gma * (1. + rc / rv)))
if (ep_max < 0.) ep_max = 0.
ep_max = Min(ep_max, pet_day)
end if
case (2) !! HARGREAVES POTENTIAL EVAPOTRANSPIRATION METHOD
!! extraterrestrial radiation
!! 37.59 is coefficient in equation 2.2.6 !!extraterrestrial
!! 30.00 is coefficient in equation 2.2.7 !!max at surface
ramm = 0.
ramm = hru_rmx(j) * 37.59 / 30.
if (tmx(j) > tmn(j)) then
pet_day = harg_petco(hru_sub(j))*(ramm / xl)*(tmpav(j) + 17.8)*
& (tmx(j) - tmn(j))**0.5
pet_day = Max(0., pet_day)
else
pet_day = 0.
endif
case (3) !! READ IN PET VALUES
pet_day = petmeas
end select
return
end