Line data Source code
1 : !> \file mo_mpr_pet.f90
2 : !> \brief \copybrief mo_mpr_pet
3 : !> \details \copydetails mo_mpr_pet
4 :
5 : !> \brief MPR routine for PET.
6 : !> \details This module sets up pet correction factor at level-1 based on LAI
7 : !> \authors Mehmet Cuneyd Demirel, Simon Stisen
8 : !> \date May 2017
9 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
10 : !! mHM is released under the LGPLv3+ license \license_note
11 : !> \ingroup f_mpr
12 : module mo_mpr_pet
13 :
14 : use mo_kind, only : i4, dp
15 :
16 : implicit none
17 :
18 : PUBLIC :: pet_correctbyLAI ! estimate PET correction factor with distributed LAI
19 : PUBLIC :: pet_correctbyASP ! estimate PET correction factor with distributed Aspect
20 : PUBLIC :: priestley_taylor_alpha ! factor (alpha) for Presley-Taylor ET estimation
21 : !PUBLIC :: aerodynamical_resistance ! aerodynamical resistance (ra) for Penman-Monteith ET estimation
22 : PUBLIC :: bulksurface_resistance ! bulk surface (stomatal) resistance (rs) for Penman-Monteith ET estimation
23 :
24 : private
25 :
26 : contains
27 :
28 :
29 : ! ----------------------------------------------------------------------------
30 :
31 : ! NAME
32 : ! pet_correctbyLAI
33 :
34 : ! PURPOSE
35 : !> \brief estimate PET correction factor based on LAI at L1
36 :
37 : !> \details estimate PET correction factor based on LAI at L1 for a given
38 : !> Leaf Area Index field.
39 :
40 : !> Global parameters needed (see mhm_parameter.nml):
41 : !> Process Case 5:
42 : !> - param(1) = PET_a_forest
43 : !> - param(2) = PET_a_impervious
44 : !> - param(3) = PET_a_pervious
45 : !> - param(4) = PET_b
46 : !> - param(5) = PET_c
47 :
48 : !> Example DSF=PET_a+PET_b*(1-exp(PET_c*LAI))
49 :
50 : !> Similar to the crop coefficient concept Kc=a+b*(1-exp(c*LAI)) by Allen, R. G., L. S. Pereira,
51 : !> D. Raes, and M. Smith (1998), Crop evapotranspiration - Guidelines for computing crop water requirements.,
52 : !> FAO
53 : !> Irrigation and drainage paper 56. See Chapter 9, Equation 97 <http://www.fao.org/docrep/X0490E/x0490e0f.htm>
54 : !> Date: 17/5/2017
55 :
56 : ! INTENT(IN)
57 : !> \param[in] "real(dp), dimension(5) :: param" parameters
58 : !> \param[in] "real(dp) :: nodata" - nodata value
59 : !> \param[in] "integer(i4), dimension(:) :: LCOVER0" Land cover at level 0
60 : !> \param[in] "real(dp), dimension(:, :) :: LAI0" LAI at level-0
61 : !> \param[in] "logical, dimension(:, :) :: mask0" mask at L0
62 : !> \param[in] "integer(i4), dimension(:) :: cell_id0" Cell ids of hi res field
63 : !> \param[in] "integer(i4), dimension(:) :: upp_row_L1" Upper row of hi res block
64 : !> \param[in] "integer(i4), dimension(:) :: low_row_L1" Lower row of hi res block
65 : !> \param[in] "integer(i4), dimension(:) :: lef_col_L1" Left column of hi res block
66 : !> \param[in] "integer(i4), dimension(:) :: rig_col_L1" Right column of hi res block
67 : !> \param[in] "integer(i4), dimension(:) :: nL0_in_L1" Number of L0 cells within a L1 cel
68 :
69 : ! INTENT(INOUT)
70 : !> \param[inout] "real(dp), dimension(:, :) :: L1_petLAIcorFactor" pet cor factor at level-1
71 :
72 : ! HISTORY
73 : !> \authors M. Cuneyd Demirel and Simon Stisen from GEUS.dk
74 :
75 : !> \date May. 2017
76 :
77 : ! Modifications:
78 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
79 :
80 16 : subroutine pet_correctbyLAI(param, nodata, LCOVER0, LAI0, mask0, cell_id0, upp_row_L1, low_row_L1, lef_col_L1, &
81 8 : rig_col_L1, nL0_in_L1, L1_petLAIcorFactor)
82 :
83 : use mo_upscaling_operators, only : upscale_harmonic_mean
84 : !$ use omp_lib
85 :
86 : implicit none
87 :
88 : ! parameters
89 : real(dp), dimension(5), intent(in) :: param
90 :
91 : ! - nodata value
92 : real(dp), intent(in) :: nodata
93 :
94 : ! Land cover at level 0
95 : integer(i4), dimension(:), intent(in) :: LCOVER0
96 :
97 : ! LAI at level-0
98 : real(dp), dimension(:, :), intent(in) :: LAI0
99 :
100 : ! mask at L0
101 : logical, dimension(:, :), intent(in) :: mask0
102 :
103 : ! Cell ids of hi res field
104 : integer(i4), dimension(:), intent(in) :: cell_id0
105 :
106 : ! Upper row of hi res block
107 : integer(i4), dimension(:), intent(in) :: upp_row_L1
108 :
109 : ! Lower row of hi res block
110 : integer(i4), dimension(:), intent(in) :: low_row_L1
111 :
112 : ! Left column of hi res block
113 : integer(i4), dimension(:), intent(in) :: lef_col_L1
114 :
115 : ! Right column of hi res block
116 : integer(i4), dimension(:), intent(in) :: rig_col_L1
117 :
118 : ! Number of L0 cells within a L1 cel
119 : integer(i4), dimension(:), intent(in) :: nL0_in_L1
120 :
121 : ! pet cor factor at level-1
122 : real(dp), dimension(:, :), intent(inout) :: L1_petLAIcorFactor
123 :
124 : ! pet cor factor at level-0
125 2234212 : real(dp), dimension(size(LCOVER0, 1), size(LAI0, 2)) :: petLAIcorFactor_0
126 :
127 : ! loop index
128 : integer(i4) :: kk, tt
129 :
130 : ! loop index
131 : integer(i4) :: LL
132 :
133 :
134 : ! ------------------------------------------------------------------
135 : ! Estimate DSF=PET_a+PET_b*(1-exp(PET_c*LAI)) to correct PET as PET=DSF*PET
136 : ! ------------------------------------------------------------------
137 : !$OMP PARALLEL
138 : !$OMP DO PRIVATE( LL ) SCHEDULE( STATIC )
139 : ! need to be done for every landcover to get DSF
140 186184 : do kk = 1, size(LCOVER0, 1)
141 :
142 186180 : LL = LCOVER0(kk)
143 :
144 : ! TODO: memory order of arrays is not optimal, how to improve?
145 4 : select case(LL)
146 : case(1) ! forest
147 864240 : petLAIcorFactor_0(kk, :) = param(1) + (param(4) * (1.0_dp - exp(param(5) * LAI0(kk, :))))
148 : case(2) ! impervious
149 152412 : petLAIcorFactor_0(kk, :) = param(2) + (param(4) * (1.0_dp - exp(param(5) * LAI0(kk, :))))
150 : case(3) ! permeable
151 1481892 : petLAIcorFactor_0(kk, :) = param(3) + (param(4) * (1.0_dp - exp(param(5) * LAI0(kk, :))))
152 : end select
153 :
154 : end do
155 : !$OMP END DO
156 : !$OMP END PARALLEL
157 :
158 52 : do tt = 1, size(LAI0, 2)
159 48 : L1_petLAIcorFactor(:, tt) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
160 1732 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata, petLAIcorFactor_0(:, tt))
161 : end do
162 :
163 :
164 4 : end subroutine pet_correctbyLAI
165 :
166 : ! ----------------------------------------------------------------------------
167 :
168 : ! NAME
169 : ! pet_correctbyASP
170 :
171 : ! PURPOSE
172 : !> \brief correction of PET
173 :
174 : !> \details Correction of PET based on L0 aspect data.
175 :
176 :
177 : !> Global parameters needed (see mhm_parameter.nml):
178 : !> - param(1) = minCorrectionFactorPET
179 : !> - param(2) = maxCorrectionFactorPET
180 : !> - param(3) = aspectTresholdPET
181 :
182 : ! INTENT(IN)
183 : !> \param[in] "integer(i4), dimension(:) :: id0" Level 0 cell id
184 : !> \param[in] "real(dp), dimension(:) :: latitude_l0" latitude on l0
185 : !> \param[in] "real(dp), dimension(:) :: Asp0" [degree] Aspect at Level 0
186 : !> \param[in] "real(dp), dimension(3) :: param" process parameters
187 : !> \param[in] "real(dp) :: nodata" - no data value
188 :
189 : ! INTENT(OUT)
190 : !> \param[out] "real(dp), dimension(:) :: fAsp0" PET correction for Aspect
191 :
192 : ! HISTORY
193 : !> \authors Stephan Thober, Rohini Kumar
194 :
195 : !> \date Dec 2012
196 :
197 : ! Modifications:
198 : ! Juliane Mai Oct 2013 - OLD parametrization --> param(1) = minCorrectionFactorPET
199 : ! --> param(2) = maxCorrectionFactorPET
200 : ! --> param(3) = aspectTresholdPET
201 : ! -------------------------------
202 : ! maxCorrectionFactorPET = minCorrectionFactorPET + delta
203 : ! -------------------------------
204 : ! NEW parametrization
205 : ! --> param(1) = minCorrectionFactorPET
206 : ! --> param(2) = delta
207 : ! --> param(3) = aspectTresholdPET
208 : ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
209 : ! Stephan Thober Sep 2015 - Mapping L1 to Lo, latitude on L0
210 : ! Luis Samaniego Sep 2015 - PET correction on the southern hemisphere
211 : ! Matthias Zink Jun 2017 - renamed and moved from mo_multi_scale_param_reg.f90 to mo_mpr_pet.f90
212 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
213 :
214 396 : subroutine pet_correctbyASP(Id0, latitude_l0, Asp0, param, nodata, fAsp0)
215 : !$ use omp_lib
216 :
217 : implicit none
218 :
219 : ! Level 0 cell id
220 : integer(i4), dimension(:), intent(in) :: id0
221 :
222 : ! latitude on l0
223 : real(dp), dimension(:), intent(in) :: latitude_l0
224 :
225 : ! - no data value
226 : real(dp), intent(in) :: nodata
227 :
228 : ! process parameters
229 : real(dp), dimension(3), intent(in) :: param
230 :
231 : ! [degree] Aspect at Level 0
232 : real(dp), dimension(:), intent(in) :: Asp0
233 :
234 : ! PET correction for Aspect
235 : real(dp), dimension(:), intent(out) :: fAsp0
236 :
237 : ! PET correction for Aspect, south
238 9165068 : real(dp), dimension(size(id0, 1)) :: fAsp0S
239 :
240 198 : logical, dimension(size(id0, 1)) :: mask_north_hemisphere_l0
241 :
242 9165068 : real(dp) :: tmp_maxCorrectionFactorPET
243 :
244 :
245 9165266 : mask_north_hemisphere_l0 = merge(.TRUE., .FALSE., latitude_l0 .gt. 0.0_dp)
246 :
247 198 : tmp_maxCorrectionFactorPET = param(1) + param(2)
248 :
249 : ! for cells on the northern hemisphere
250 : !$OMP PARALLEL
251 0 : fAsp0 = merge(&
252 : param(1) + (tmp_maxCorrectionFactorPET - param(1)) / param(3) * asp0, &
253 : param(1) + (tmp_maxCorrectionFactorPET - param(1)) / (360._dp - param(3)) * (360._dp - Asp0), &
254 : ! ( asp0 < param(3) ) .and. mask_north_hemisphere_l0 )
255 9165068 : asp0 < param(3))
256 9165068 : fAsp0 = merge(fAsp0, nodata, Id0 /= int(nodata, i4))
257 : !$OMP END PARALLEL
258 :
259 : ! for cells on the southern hemisphere
260 : !$OMP PARALLEL
261 198 : fAsp0S = merge(&
262 : param(1) + (tmp_maxCorrectionFactorPET - param(1)) / (360._dp - param(3)) * (360._dp - Asp0), &
263 : param(1) + (tmp_maxCorrectionFactorPET - param(1)) / param(3) * asp0, &
264 9165068 : asp0 < param(3))
265 9165266 : fAsp0S = merge(fAsp0S, nodata, Id0 /= int(nodata, i4))
266 : !$OMP END PARALLEL
267 :
268 : !$OMP PARALLEL
269 9165266 : fAsp0 = merge(fAsp0, fAsp0S, mask_north_hemisphere_l0)
270 : !$OMP END PARALLEL
271 :
272 4 : end subroutine pet_correctbyASP
273 :
274 : ! ----------------------------------------------------------------------------
275 :
276 : ! NAME
277 : ! priestley_taylor_alpha
278 :
279 : ! PURPOSE
280 : !> \brief Regionalization of priestley taylor alpha
281 :
282 : !> \details estimation of priestley taylor alpha
283 : !> Global parameters needed (see mhm_parameter.nml):
284 : !> - param(1) = PriestleyTaylorCoeff
285 : !> - param(2) = PriestleyTaylorLAIcorr
286 :
287 : ! INTENT(IN)
288 : !> \param[in] "real(dp), dimension(:, :) :: LAI0" LAI at level-0
289 : !> \param[in] "real(dp), dimension(:) :: param" input parameter
290 : !> \param[in] "logical, dimension(:, :) :: mask0" mask at level 0
291 : !> \param[in] "real(dp) :: nodata" - nodata value
292 : !> \param[in] "integer(i4), dimension(:) :: cell_id0" Cell ids of hi res field
293 : !> \param[in] "integer(i4), dimension(:) :: nL0_in_L1" number of l0 cells within a l1 cell
294 : !> \param[in] "integer(i4), dimension(:) :: Upp_row_L1" upper row of a l1 cell in l0 grid
295 : !> \param[in] "integer(i4), dimension(:) :: Low_row_L1" lower row of a l1 cell in l0 grid
296 : !> \param[in] "integer(i4), dimension(:) :: Lef_col_L1" left col of a l1 cell in l0 grid
297 : !> \param[in] "integer(i4), dimension(:) :: Rig_col_L1" right col of a l1 cell in l0 grid
298 :
299 : ! INTENT(OUT)
300 : !> \param[out] "real(dp), dimension(:, :) :: priestley_taylor_alpha1" bulk surface resistance
301 :
302 : ! HISTORY
303 : !> \authors Matthias Zink
304 :
305 : !> \date Apr 2013
306 :
307 : ! Modifications:
308 : ! Matthias Zink Jun 2017 - moved from mo_multi_scale_param_reg.f90 to mo_mpr_pet.f90
309 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
310 :
311 3 : subroutine priestley_taylor_alpha(LAI0, param, mask0, nodata, cell_id0, nL0_in_L1, Upp_row_L1, Low_row_L1, Lef_col_L1, &
312 2 : Rig_col_L1, priestley_taylor_alpha1)
313 :
314 198 : use mo_upscaling_operators, only : upscale_arithmetic_mean
315 :
316 : implicit none
317 :
318 : ! LAI at level-0
319 : real(dp), dimension(:, :), intent(in) :: LAI0
320 :
321 : ! input parameter
322 : real(dp), dimension(:), intent(in) :: param
323 :
324 : ! mask at level 0
325 : logical, dimension(:, :), intent(in) :: mask0
326 :
327 : ! - nodata value
328 : real(dp), intent(in) :: nodata
329 :
330 : ! Cell ids of hi res field
331 : integer(i4), dimension(:), intent(in) :: cell_id0
332 :
333 : ! number of l0 cells within a l1 cell
334 : integer(i4), dimension(:), intent(in) :: nL0_in_L1
335 :
336 : ! upper row of a l1 cell in l0 grid
337 : integer(i4), dimension(:), intent(in) :: Upp_row_L1
338 :
339 : ! lower row of a l1 cell in l0 grid
340 : integer(i4), dimension(:), intent(in) :: Low_row_L1
341 :
342 : ! left col of a l1 cell in l0 grid
343 : integer(i4), dimension(:), intent(in) :: Lef_col_L1
344 :
345 : ! right col of a l1 cell in l0 grid
346 : integer(i4), dimension(:), intent(in) :: Rig_col_L1
347 :
348 : ! bulk surface resistance
349 : real(dp), dimension(:, :), intent(out) :: priestley_taylor_alpha1
350 :
351 : integer(i4) :: tt
352 :
353 : ! dim 1 = number of cells on level 0, time
354 1 : real(dp), dimension(:, :), allocatable :: priestley_taylor_alpha0
355 :
356 :
357 : ! initialize some things
358 558557 : allocate(priestley_taylor_alpha0 (size(LAI0, 1), size(LAI0, 2))) ; priestley_taylor_alpha0 = nodata
359 421 : priestley_taylor_alpha1 = nodata
360 : !
361 13 : do tt = 1, size(LAI0, 2)
362 558552 : priestley_taylor_alpha0(:, tt) = param(1) + param(2) * LAI0(:, tt)
363 :
364 : priestley_taylor_alpha1(:, tt) = upscale_arithmetic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
365 13 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata, priestley_taylor_alpha0(:, tt))
366 :
367 : end do
368 :
369 1 : end subroutine priestley_taylor_alpha
370 :
371 : ! ----------------------------------------------------------------------------
372 :
373 : ! NAME
374 : ! bulksurface_resistance
375 :
376 : ! PURPOSE
377 : !> \brief Regionalization of bulk surface resistance
378 :
379 : !> \details estimation of bulk surface resistance
380 : !> Global parameters needed (see mhm_parameter.nml):
381 : !> - param(1) = stomatal_resistance
382 :
383 : ! INTENT(IN)
384 : !> \param[in] "real(dp), dimension(:, :) :: LAI0" LAI at level-0
385 : !> \param[in] "real(dp) :: param" - global parameter
386 : !> \param[in] "logical, dimension(:, :) :: mask0" mask at level 0
387 : !> \param[in] "real(dp) :: nodata" - nodata value
388 : !> \param[in] "integer(i4), dimension(:) :: cell_id0" Cell ids of hi res field
389 : !> \param[in] "integer(i4), dimension(:) :: nL0_in_L1" number of l0 cells within a l1 cell
390 : !> \param[in] "integer(i4), dimension(:) :: Upp_row_L1" upper row of a l1 cell in l0 grid
391 : !> \param[in] "integer(i4), dimension(:) :: Low_row_L1" lower row of a l1 cell in l0 grid
392 : !> \param[in] "integer(i4), dimension(:) :: Lef_col_L1" left col of a l1 cell in l0 grid
393 : !> \param[in] "integer(i4), dimension(:) :: Rig_col_L1" right col of a l1 cell in l0 grid
394 :
395 : ! INTENT(OUT)
396 : !> \param[out] "real(dp), dimension(:, :) :: bulksurface_resistance1" bulk surface resistance
397 :
398 : ! HISTORY
399 : !> \authors Matthias Zink
400 :
401 : !> \date Apr 2013
402 :
403 : ! Modifications:
404 : ! Matthias Zink Jun 2017 - moved from mo_multi_scale_param_reg.f90 to mo_mpr_pet.f90
405 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
406 :
407 18 : subroutine bulksurface_resistance(LAI0, param, mask0, nodata, cell_id0, nL0_in_L1, Upp_row_L1, Low_row_L1, Lef_col_L1, &
408 12 : Rig_col_L1, bulksurface_resistance1)
409 :
410 1 : use mo_mpr_constants, only : LAI_factor_surfResi, LAI_offset_surfResi, max_surfResist
411 : use mo_upscaling_operators, only : upscale_arithmetic_mean
412 :
413 : implicit none
414 :
415 : ! LAI at level-0
416 : real(dp), dimension(:, :), intent(in) :: LAI0
417 :
418 : ! - global parameter
419 : real(dp), intent(in) :: param
420 :
421 : ! mask at level 0
422 : logical, dimension(:, :), intent(in) :: mask0
423 :
424 : ! - nodata value
425 : real(dp), intent(in) :: nodata
426 :
427 : ! Cell ids of hi res field
428 : integer(i4), dimension(:), intent(in) :: cell_id0
429 :
430 : ! number of l0 cells within a l1 cell
431 : integer(i4), dimension(:), intent(in) :: nL0_in_L1
432 :
433 : ! upper row of a l1 cell in l0 grid
434 : integer(i4), dimension(:), intent(in) :: Upp_row_L1
435 :
436 : ! lower row of a l1 cell in l0 grid
437 : integer(i4), dimension(:), intent(in) :: Low_row_L1
438 :
439 : ! left col of a l1 cell in l0 grid
440 : integer(i4), dimension(:), intent(in) :: Lef_col_L1
441 :
442 : ! right col of a l1 cell in l0 grid
443 : integer(i4), dimension(:), intent(in) :: Rig_col_L1
444 :
445 : ! bulk surface resistance
446 : real(dp), dimension(:, :), intent(out) :: bulksurface_resistance1
447 :
448 : integer(i4) :: tt
449 :
450 : ! dim 1 = number of cells on level 0,
451 : ! dim 2 = number of months in year (12)
452 6 : real(dp), dimension(:, :), allocatable :: bulksurface_resistance0
453 :
454 :
455 : ! initialize some things
456 3351342 : allocate(bulksurface_resistance0 (size(LAI0, 1), size(LAI0, 2))) ; bulksurface_resistance0 = nodata
457 2526 : bulksurface_resistance1 = nodata
458 : !
459 :
460 78 : do tt = 1, size(LAI0, 2)
461 216 : bulksurface_resistance0(:, tt) = param / (LAI0(:, tt) / &
462 3351456 : (LAI_factor_surfResi * LAI0(:, tt) + LAI_offset_surfResi))
463 : ! efeective LAI from McMahon et al ,2013 , HESS supplements
464 :
465 : ! since LAI may be very low, rs becomes very high
466 : ! thus the values are restricted to maximum literaure values (i.e. McMahon et al ,2013 , HESS)
467 : bulksurface_resistance0(:, tt) = merge(max_surfResist, bulksurface_resistance0(:, tt), &
468 3351312 : bulksurface_resistance0(:, tt) .GT. max_surfResist)
469 :
470 : bulksurface_resistance1(:, tt) = upscale_arithmetic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
471 78 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata, bulksurface_resistance0(:, tt))
472 :
473 : end do
474 :
475 6 : end subroutine bulksurface_resistance
476 :
477 : end module mo_mpr_pet
|
