Line data Source code
1 : !> \file mo_pet.f90
2 : !> \brief \copybrief mo_pet
3 : !> \details \copydetails mo_pet
4 :
5 : !> \brief Module for calculating reference/potential evapotranspiration [mm d-1]
6 : !> \details This module calculates PET [mm/d] based on one of the methods
7 : !! - Hargreaves-Samani (1982)
8 : !! - Priestly-Taylor (1972)
9 : !! - Penman-Monteith FAO (1998)
10 : !> \authors Matthias Zink, Christoph Schneider, Matthias Cuntz
11 : !> \date Apr 2014
12 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
13 : !! mHM is released under the LGPLv3+ license \license_note
14 : !> \ingroup f_mhm
15 : MODULE mo_pet
16 :
17 : ! This module is for the UFZ CHS mesoscale hydrologic model mHM.
18 :
19 : USE mo_kind, ONLY : i4, dp
20 :
21 : IMPLICIT NONE
22 :
23 : PUBLIC :: pet_hargreaves ! Hargreaves-Samani
24 : PUBLIC :: pet_priestly ! Priestley-Taylor
25 : PUBLIC :: pet_penman ! Penman-Monteith
26 : PUBLIC :: slope_satpressure
27 : PUBLIC :: extraterr_rad_approx
28 : PUBLIC :: sat_vap_pressure
29 :
30 :
31 : ! ------------------------------------------------------------------
32 :
33 : CONTAINS
34 :
35 : ! ------------------------------------------------------------------
36 :
37 : ! NAME
38 : ! pet_hargreaves
39 :
40 : ! PURPOSE
41 : !> \brief Reference Evapotranspiration after Hargreaves
42 :
43 : !> \details Calculates the Reference Evapotranspiration \f$ [mm\;d^{-1}] \f$ based on the Hargreaves-Samani
44 : !> (1982)
45 : !> model for a given cell by applying the equation
46 : !> \f[ PET = HarSamCoeff * R_a * (T_{avg} + HarSamConst) * \sqrt{ T_{max} - T_{min}} \f]
47 : !> where \f$ R_a\;[W\;m^{-2}]\f$ is the incoming solar radiation and
48 : !> \f$ T_{avg}, T_{max} \f$ and \f$ T_{min}\f$ \f$ [ ^{\circ}C]\f$ are the mean, maximum,
49 : !> and minimum daily temperatures at the given day, respectively.
50 : !> \note Hargreaves, G.H., and Samani, Z.A. (1982). "Estimating potential evapotranspiration."
51 : !> Tech. Note, J. Irrig. and drain. Engrg., ASCE, 108(3):225-230.
52 :
53 : ! INTENT(IN)
54 : !> \param[in] "real(dp) :: HarSamCoeff" coefficient of Hargreaves-Samani equation [-]
55 : !> \param[in] "real(dp) :: HarSamConst" constant of Hargreaves-Samani equation [-]
56 : !> \param[in] "real(dp) :: tavg" daily men temperature \f$[^{\circ}C]\f$
57 : !> \param[in] "real(dp) :: tmax" daily maximum of temp \f$[^{\circ}C]\f$
58 : !> \param[in] "real(dp) :: tmin" daily minimum of temp \f$[^{\circ}C]\f$
59 : !> \param[in] "real(dp) :: latitude" latitude of the cell for Ra estimation \f$[radians]\f$
60 : !> \param[in] "integer(i4) :: doy" day of year for Ra estimation
61 :
62 : ! RETURN
63 : !> \return real(dp) :: pet_hargreaves — Hargreaves-Samani pot. evapotranspiration [mm d-1]
64 :
65 : ! HISTORY
66 : !> \authors Matthias Zink
67 :
68 : !> \date Dec 2012
69 :
70 : ! Modifications:
71 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
72 :
73 893520 : elemental pure FUNCTION pet_hargreaves(HarSamCoeff, HarSamConst, tavg, tmax, tmin, latitude, doy)
74 :
75 : use mo_constants, only : deg2rad_dp
76 : use mo_utils, only : LE
77 :
78 : implicit none
79 :
80 : ! coefficient of Hargreaves-Samani equation [-]
81 : real(dp), intent(in) :: HarSamCoeff
82 :
83 : ! constant of Hargreaves-Samani equation [-]
84 : real(dp), intent(in) :: HarSamConst
85 :
86 : ! daily men temperature \f$[^{\circ}C]\f$
87 : real(dp), intent(in) :: tavg
88 :
89 : ! daily maximum of temp \f$[^{\circ}C]\f$
90 : real(dp), intent(in) :: tmax
91 :
92 : ! daily minimum of temp \f$[^{\circ}C]\f$
93 : real(dp), intent(in) :: tmin
94 :
95 : ! latitude of the cell for Ra estimation \f$[radians]\f$
96 : real(dp), intent(in) :: latitude
97 :
98 : ! day of year for Ra estimation
99 : integer(i4), intent(in) :: doy
100 :
101 : ! reference evapotranspiration in [mm d-1]
102 : real(dp) :: pet_hargreaves
103 :
104 : ! tmax-Tmin
105 : ! correction for shity input data (tmax<tmin) and to avoid numerical errors ! MZMZMZMZ
106 893520 : real(dp) :: delta_temp
107 :
108 :
109 893520 : delta_temp = tmax - tmin
110 893520 : if(LE(delta_temp, 0.0_dp) .or. LE(tavg, -HarSamConst)) then
111 : pet_hargreaves = 0.0_dp
112 : else
113 893520 : pet_hargreaves = HarSamCoeff * extraterr_rad_approx(doy, deg2rad_dp * latitude) * (tavg + HarSamConst) * sqrt(delta_temp)
114 : end if
115 :
116 893520 : END FUNCTION pet_hargreaves
117 :
118 :
119 : ! ------------------------------------------------------------------
120 :
121 : ! NAME
122 : ! pet_priestly
123 :
124 : ! PURPOSE
125 : !> \brief Reference Evapotranspiration after Priestly-Taylor
126 :
127 : !> \details Calculates the Reference Evapotranspiration \f$ [mm\;d^{-1}] \f$ based on the
128 : !> Priestly-Taylor (1972) model for every given cell by applying the equation
129 : !> \f[ PET = \alpha * \frac{\Delta}{(\gamma + \Delta)} * R_n \f]
130 : !> where \f$R_n\;[W\;m^{-2}]\f$ is the net solar radiation \f$\Delta = f(T_{avg})\f$ is the slope
131 : !> of the saturation-vapour pressure curve and \f$\alpha\f$ is a emperical coefficient.
132 :
133 : ! INTENT(IN)
134 : !> \param[in] "real(dp) :: PrieTayParam" Priestley-Taylor coefficient \f$ \alpha [-] \f$
135 : !> \param[in] "real(dp) :: Rn" net solar radiation \f$ [W\;m^{-2}] \f$
136 : !> \param[in] "real(dp) :: Tavg"
137 :
138 : ! RETURN
139 : !> \return real(dp) :: pet_priestly — Priestley-Taylor pot. evapotranspiration [mm d-1]
140 :
141 : ! HISTORY
142 : !> \authors Matthias Zink
143 :
144 : !> \date Apr 2014
145 :
146 : ! Modifications:
147 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
148 :
149 445536 : elemental pure FUNCTION pet_priestly(PrieTayParam, Rn, tavg)
150 :
151 893520 : use mo_constants, only : DaySecs, Psychro_dp, SpecHeatET_dp
152 :
153 : implicit none
154 :
155 : ! Priestley-Taylor coefficient \f$ \alpha [-] \f$
156 : real(dp), intent(in) :: PrieTayParam
157 :
158 : ! net solar radiation \f$ [W\;m^{-2}] \f$
159 : real(dp), intent(in) :: Rn
160 :
161 : real(dp), intent(in) :: Tavg
162 :
163 : ! reference evapotranspiration in [mm d-1]
164 : real(dp) :: pet_priestly
165 :
166 : ! save slope of saturation vapor pressure curve
167 445536 : real(dp) :: delta
168 :
169 :
170 445536 : delta = slope_satpressure(Tavg) ! slope of saturation vapor pressure curve
171 : ! in [mm d-1]
172 445536 : pet_priestly = PrieTayParam * delta / (Psychro_dp + delta) * (Rn * DaySecs / SpecHeatET_dp)
173 :
174 445536 : END FUNCTION pet_priestly
175 :
176 :
177 : ! ------------------------------------------------------------------
178 :
179 : ! NAME
180 : ! pet_penman
181 :
182 : ! PURPOSE
183 : !> \brief Reference Evapotranspiration after Penman-Monteith
184 :
185 : !> \details Calculates the reference evapotranspiration \f$ [mm\;d^{-1}] \f$ based on the
186 : !> Penman-Monteith model for every given cell by applying the equation
187 : !> \f[ PET = \frac{1}{\lambda} \cdot
188 : !> \frac{\Delta \cdot R_n + \rho \cdot c_p \cdot (e_s-e) \cdot \frac{a_sh}{r_a}}
189 : !> {\Delta + \gamma \cdot \frac{a_sh}{a_s} \cdot \left( 1 + \frac{r_s}{r_a} \right) } \f]
190 : !> where \f$R_n\;[W\;m^{-2}]\f$ is the net solar radiation,
191 : !> \f$\Delta\;[kPa\;K^{-1}]\f$ is the slope of the saturation-vapour pressure curve,
192 : !> \f$ \lambda\;[MJ\;kg^{-1}] \f$ is the latent heat of vaporization,
193 : !> \f$ (e_s-e)\;[kPa] \f$ is the vapour pressure deficit of the air,
194 : !> \f$ \rho\;[kg\;m^{-3}] \f$ is the mean atmospheric density,
195 : !> \f$ c_p=1005.0\;J\;kg^{-1}\;K^{-1} \f$ is the specific heat of the air,
196 : !> \f$ \gamma [kPa\;K^{-1}] \f$ is the psychrometric constant,
197 : !> \f$ r_s [s m^{-1}] \f$ is the bulk canopy resistance,
198 : !> \f$ r_a [s m^{-1}] \f$ is the aerodynamic resistance,
199 : !> \f$ a_s [1] \f$ is the fraction of one-sided leaf area covered by stomata
200 : !> (1 if stomata are on one side only, 2 if they are on both sides) and
201 : !> \f$ a_{sh} [-] \f$ is the fraction of projected area exchanging sensible heat with the air (2)
202 : !> Implementation refers to the so-called Penman-Montheith equation for transpiration.
203 : !> Adjusting the arguments \f$ a_{sh} \f$ and \f$ a_s \f$ we obtain the corrected MU equation (for details
204 : !> see Schymanski and Or, 2017). If \f$ a_{sh} = 1 = a_s \f$ Penman-Montheith equation for transpiration
205 : !> is preserved. For reproducing characteristics of symmetrical amphistomatous leaves use
206 : !> \f$ a_{sh} = 2 = a_s \f$, in which case the classic PM equation is only missing a factor
207 : !> of 2 in the nominator, as pointed out by Jarvis and McNaughton (1986, Eq. A9).
208 : !> These analytical solutions eliminated the non-linearity problem of the saturation vapour pressure curve,
209 : !> but they do not consider the dependency of the long-wave component of the soil surface or leaf energy balance
210 : !> (\f$ R_l \f$) on soil or leaf temperature (\f$ T_l \f$). We assume that net radiation
211 : !> equals the absorbed short-wave radiation, i.e. \f$ R_N = R_s \f$ (p.79 in Monteith and Unsworth, 2013).
212 :
213 : ! INTENT(IN)
214 : !> \param[in] "real(dp) :: net_rad" net radiation \f$[W m^{-2}]\f$
215 : !> \param[in] "real(dp) :: tavg" average daily temperature \f$[^{\circ}C]\f$
216 : !> \param[in] "real(dp) :: act_vap_pressure" actual vapur pressure \f$[kPa]\f$
217 : !> \param[in] "real(dp) :: aerodyn_resistance" aerodynmaical resistance \f$s\;m^{-1}\f$
218 : !> \param[in] "real(dp) :: bulksurface_resistance" bulk surface resistance \f$s\;m^{-1}\f$
219 : !> \param[in] "real(dp) :: a_s" fraction of one-sided leaf area covered by stomata \f$1\f$
220 : !> \param[in] "real(dp) :: a_sh" fraction of projected area exchanging sensible heat with the
221 : !> air \f$1\f$
222 :
223 : ! RETURN
224 : !> \return real(dp) :: pet_penman — Reference Evapotranspiration [mm d-1]
225 :
226 : ! HISTORY
227 : !> \authors Matthias Zink
228 :
229 : !> \date Apr 2014
230 :
231 : ! Modifications:
232 : ! Johannes Brenner Nov 2017 - include arguments a_s and a_sh to enable corrected MU approach
233 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
234 :
235 2673216 : elemental pure FUNCTION pet_penman(net_rad, tavg, act_vap_pressure, aerodyn_resistance, bulksurface_resistance, a_s, &
236 : a_sh)
237 :
238 445536 : use mo_constants, only : DaySecs, Psychro_dp, SpecHeatET_dp, cp0_dp, rho0_dp
239 :
240 : implicit none
241 :
242 : ! net radiation \f$[W m^{-2}]\f$
243 : real(dp), intent(in) :: net_rad
244 :
245 : ! average daily temperature \f$[^{\circ}C]\f$
246 : real(dp), intent(in) :: tavg
247 :
248 : ! actual vapur pressure \f$[kPa]\f$
249 : real(dp), intent(in) :: act_vap_pressure
250 :
251 : ! aerodynmaical resistance \f$s\;m^{-1}\f$
252 : real(dp), intent(in) :: aerodyn_resistance
253 :
254 : ! bulk surface resistance \f$s\;m^{-1}\f$
255 : real(dp), intent(in) :: bulksurface_resistance
256 :
257 : ! fraction of one-sided leaf area covered by stomata \f$1\f$
258 : real(dp), intent(in) :: a_s
259 :
260 : ! fraction of projected area exchanging sensible heat with the air \f$1\f$
261 : real(dp), intent(in) :: a_sh
262 :
263 : ! reference evapotranspiration in [mm d-1]
264 : real(dp) :: pet_penman
265 :
266 :
267 : pet_penman = DaySecs / SpecHeatET_dp * & ! conversion factor [W m-2] to [mm d-1]
268 : (slope_satpressure(tavg) * net_rad + &
269 : rho0_dp * cp0_dp * (sat_vap_pressure(tavg) - act_vap_pressure) * a_sh / aerodyn_resistance) / &
270 2673216 : (slope_satpressure(tavg) + Psychro_dp * a_sh / a_s * (1.0_dp + bulksurface_resistance / aerodyn_resistance))
271 :
272 2673216 : END FUNCTION pet_penman
273 :
274 : ! ------------------------------------------------------------------
275 :
276 : ! NAME
277 : ! extraterr_rad_approx
278 :
279 : ! PURPOSE
280 : !> \brief Approximation of extraterrestrial radiation
281 :
282 : !> \details Approximation of extraterrestrial radiation at the top of the atmosphere \f$ R_a \f$
283 : !> after Duffie and Beckman (1980).
284 : !> \f$ R_a \f$ is converted from \f$ [J\;m^{-2}\;d^{-1}] \f$ in \f$ [mm\;d^{-1}]\f$ .
285 : !> \f[ R_a = \frac{86400}{ \pi \cdot \lambda} \cdot E_0 \cdot
286 : !> d_r \cdot (\omega \cdot \sin(latitude) \cdot \sin(\delta) + \cos(latitude) \cdot \cos(\delta) \cdot
287 : !> \sin(\omega) \f]
288 : !> where \f$ E_0=1367\;J\;m^{-2}\;s^{-1} \f$ is the solar constant and
289 : !> It is dependent on the following sub equations:
290 : !> The relative distance Earth-Sun:
291 : !> \f[ d_r = 1 + 0.033 \cdot \cos \left( \frac{2 \cdot \pi \cdot doy}{365} \right) \f]
292 : !> in which doy is the day of the year.
293 : !> The solar declination [radians] defined by
294 : !> \f[ \delta = 0.4093 \cdot \sin\left( \frac{2 \cdot \pi \cdot doy}{365} - 1.405 \right) \f]
295 : !> The sunset hour angle [radians]:
296 : !> \f[ \omega = \arccos( - \tan(latitude) * \tan(\delta) ) \f]
297 : !> < \f$ \lambda = 2.45 \cdot 10^6\;J\;m^{-2}\;mm^{-1} \f$ is the latent heat of vaporization.
298 :
299 : ! INTENT(IN)
300 : !> \param[in] "integer(i4) :: doy" day of year [-]
301 : !> \param[in] "real(dp) :: latitude" latitude [rad]
302 :
303 : ! RETURN
304 : !> \return real(dp) :: extraterr_rad_approx — extraterrestrial radiation approximation \f$[W\;m^{-2}]\f$
305 :
306 : ! HISTORY
307 : !> \authors Matthias Zink
308 :
309 : !> \date Apr 2014
310 :
311 : ! Modifications:
312 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
313 :
314 893520 : elemental pure FUNCTION extraterr_rad_approx(doy, latitude)
315 :
316 2673216 : use mo_constants, only : DaySecs, YearDays, PI_D, SolarConst_dp, SpecHeatET_dp, TWOPI_D
317 : use mo_mhm_constants, only : DuffieDelta1, DuffieDelta2, DuffieDr
318 :
319 : implicit none
320 :
321 : ! day of year [-]
322 : integer(i4), intent(in) :: doy
323 :
324 : ! latitude [rad]
325 : real(dp), intent(in) :: latitude
326 :
327 : ! extraterrestrial radiation
328 : real(dp) :: extraterr_rad_approx
329 :
330 893520 : real(dp) :: dr, delta
331 :
332 893520 : real(dp) :: omega
333 :
334 893520 : real(dp) :: arg
335 :
336 :
337 : ! inverse relative distance Earth-Sun - correction for eccentricity of Earths orbit around the sun
338 893520 : dr = 1.0_dp + DuffieDr * cos(TWOPI_D * doy / YearDays)
339 : ! declination of the sun above the celestial equator in radians
340 893520 : delta = DuffieDelta1 * sin(TWOPI_D * doy / YearDays - DuffieDelta2)
341 :
342 : ! arccos(x) is only defined between PI and 0 (for x between -1 and 1)
343 : ! check limits
344 893520 : arg = - tan(latitude) * tan(delta)
345 0 : if(arg .lt. -1.0_dp) arg = -1.0_dp
346 893520 : if(arg .gt. 1.0_dp) arg = 1.0_dp
347 :
348 : ! sunrise hour angle in radians
349 893520 : omega = acos(arg)
350 :
351 : ! Ra - converted from [J m-2 d-1] in [mm d-1]
352 : extraterr_rad_approx = DaySecs / PI_D / SpecHeatET_dp * SolarConst_dp * &
353 893520 : dr * (omega * sin(latitude) * sin(delta) + cos(latitude) * cos(delta) * sin(omega))
354 :
355 893520 : end FUNCTION extraterr_rad_approx
356 :
357 :
358 : ! ------------------------------------------------------------------
359 :
360 : ! NAME
361 : ! slope_satpressure
362 :
363 : ! PURPOSE
364 : !> \brief slope of saturation vapour pressure curve
365 :
366 : !> \details slope of saturation vapour pressure curve after Tetens
367 : !> \f[ \Delta = \frac{0.6108 * e_s(T_a)}{e^(2 \cdot \log(T_a + 237.3))} \f]
368 :
369 : ! INTENT(IN)
370 : !> \param[in] "real(dp) :: tavg" average daily temperature \f$[^{\circ}C]\f$
371 :
372 : ! RETURN
373 : !> \return real(dp) :: slope_satpressure — slope of saturation vapour pressure curve
374 : !> \f$[kPa\;K{-1}]\f$
375 :
376 : ! HISTORY
377 : !> \authors Matthias Zink
378 :
379 : !> \date Apr 2014
380 :
381 : ! Modifications:
382 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
383 :
384 8668237 : elemental pure FUNCTION slope_satpressure(tavg)
385 :
386 893520 : use mo_mhm_constants, only : satpressureslope1, tetens_c3
387 :
388 : implicit none
389 :
390 : ! average daily temperature \f$[^{\circ}C]\f$
391 : real(dp), intent(in) :: tavg
392 :
393 : ! slope of saturation vapour pressure curve
394 : real(dp) :: slope_satpressure
395 :
396 :
397 8668237 : slope_satpressure = satpressureslope1 * sat_vap_pressure(tavg) / exp(2.0_dp * log(Tavg + tetens_c3))
398 :
399 8668237 : END FUNCTION slope_satpressure
400 :
401 : ! ------------------------------------------------------------------
402 :
403 : ! NAME
404 : ! sat_vap_pressure
405 :
406 : ! PURPOSE
407 : !> \brief calculation of the saturation vapour pressure
408 :
409 : !> \details Calculation of the saturation vapour pressure
410 : !> \f[ e_s(T_a) = 0.6108 \cdot \exp \left( \frac{17.27 \cdot T_a}{T_a + 237.3} \right) \f]
411 :
412 : ! INTENT(IN)
413 : !> \param[in] "real(dp) :: tavg" temperature [degC]
414 :
415 : ! RETURN
416 : !> \return real(dp) :: sat_vap_pressure — saturation vapour pressure [kPa]
417 :
418 : ! HISTORY
419 : !> \authors Matthias Zink
420 :
421 : !> \date Apr 2014
422 :
423 : ! Modifications:
424 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
425 :
426 11341453 : elemental pure FUNCTION sat_vap_pressure(tavg)
427 :
428 8668237 : use mo_mhm_constants, only : tetens_c1, tetens_c2, tetens_c3
429 :
430 : implicit none
431 :
432 : ! temperature [degC]
433 : real(dp), intent(in) :: tavg
434 :
435 : ! saturation vapour pressure [kPa]
436 : real(dp) :: sat_vap_pressure
437 :
438 :
439 11341453 : sat_vap_pressure = tetens_c1 * exp(tetens_c2 * tavg / (tavg + tetens_c3))
440 :
441 11341453 : END FUNCTION sat_vap_pressure
442 : !
443 : END MODULE mo_pet
|
