5.13.2-dev0
mHM
The mesoscale Hydrological Model
Loading...
Searching...
No Matches
mo_multi_param_reg.f90
Go to the documentation of this file.
1!> \file mo_multi_param_reg.f90
2!> \brief \copybrief mo_multi_param_reg
3!> \details \copydetails mo_multi_param_reg
4
5!> \brief Multiscale parameter regionalization (MPR).
6!> \details This module provides the routines for multiscale parameter regionalization (MPR).
7!> \authors Stephan Thober, Rohini Kumar
8!> \date Dec 2012
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
12MODULE mo_multi_param_reg
13
14 use mo_kind, only : i4, dp
15 use mo_common_constants, only : nodata_dp, nodata_i4
16 use mo_message, only : message, error_message
17
18 implicit none
19
20 private
21
22 PUBLIC :: mpr ! calculates effective regionalised parameters
23 PUBLIC :: canopy_intercept_param ! estimate effective max. canopy interception
24
25contains
26
27 !> \brief Regionalizing and Upscaling process parameters
28 !> \details calculating process parameters at L0 scale (Regionalization), like:
29 !! - Baseflow recession parameter
30 !! - Soil moisture parameters
31 !! - PET correction for aspect
32 !!
33 !! and upscale these parameters to retrieve effective parameters at scale L1.
34 !! Further parameter regionalizations are done for:
35 !! - snow accumulation and melting parameters
36 !! - threshold parameter for runoff generation on impervious layer
37 !! - karstic percolation loss
38 !! - setting up the Regionalized Routing Parameters
39 !!
40 !> \changelog
41 !! - Stephan Thober Jan 2013
42 !! - updated calling sequence for upscaling operators
43 !! - Luis Samaniego Feb 2013
44 !! - calling sequence, initial CHECK, call mpr_runoff
45 !! - Stephan Thober Feb 2013
46 !! - added subroutine for karstic percolation loss removed L1_, L0_ in variable names
47 !! - Stephan Thober Aug 2015
48 !! - moved regionalization of routing to mRM
49 !! - Rohini Kumar Mar 2016
50 !! - changes for handling multiple soil database options
51 !! - Zink M. & Demirel M.C. Mar 2017
52 !! - Added Jarvis soil water stress function at SM process(3)
53 !! - Demirel M.C. & S. Stisen Apr 2017
54 !! - Added FC dependency on root fraction coefficient at SM process(3)
55 !! - Robert Schweppe Dec 2017
56 !! - added loop over LCscenes inside MPR, renamed variables rewrite
57 !! - Robert Schweppe Jun 2018
58 !! - refactoring and reformatting
59 !! - Demirel M.C. & S. Stisen Jun 2020
60 !! - Added Feddes and global FC dependency on root fraction coefficient at SM process(3)=4
61 !! - Rohini Kumar Oct 2021
62 !! - Added Neutron count module to mHM integrate into develop branch (5.11.2)
63 !! - Sebastian Müller Mar 2023
64 !! - made L1_alpha, L1_kSlowFlow, L1_kBaseFlow and L1_kPerco land cover dependent
65 !> \authors Stephan Thober, Rohini Kumar
66 !> \date Dec 2012
67subroutine mpr(mask0, geoUnit0, soilId0, Asp0, gridded_LAI0, LCover0, slope_emp0, y0, Id0, upper_bound1, lower_bound1, &
68 left_bound1, right_bound1, n_subcells1, fSealed1, alpha1, degDayInc1, degDayMax1, degDayNoPre1, fAsp1, &
69 HarSamCoeff1, PrieTayAlpha1, aeroResist1, surfResist1, fRoots1, kFastFlow1, kSlowFlow1, kBaseFlow1, &
70 kPerco1, karstLoss1, soilMoistFC1, soilMoistSat1, soilMoistExp1, jarvis_thresh_c1, tempThresh1, &
71 unsatThresh1, sealedThresh1, wiltingPoint1, maxInter1, petLAIcorFactor, &
72 No_Count1, bulkDens1, latticeWater1, COSMICL31, &
73 parameterset )
74
75 use mo_common_variables, only : global_parameters, processmatrix
76 use mo_mpr_smhorizons, only : mpr_smhorizons
77 use mo_mpr_global_variables, only : horizondepth_mhm, fracsealed_cityarea, iflag_soildb, nsoilhorizons_mhm, &
78 soildb
79 use mo_mpr_pet, only : bulksurface_resistance, pet_correctbyasp, pet_correctbylai, priestley_taylor_alpha
80 use mo_mpr_runoff, only : mpr_runoff
81 use mo_mpr_soilmoist, only : mpr_sm
82 use mo_upscaling_operators, only : l0_fractionalcover_in_lx, &
83 upscale_arithmetic_mean
84 use mo_mpr_neutrons, only: mpr_neutrons
85 implicit none
86
87 !> mask at level 0 field
88 logical, dimension(:, :), intent(in) :: mask0
89 !> L0 geological units
90 integer(i4), dimension(:), intent(in) :: geounit0
91 !> soil Ids at level 0
92 integer(i4), dimension(:, :), intent(in) :: soilid0
93 !> [degree] Aspect at Level 0
94 real(dp), dimension(:), intent(in) :: asp0
95 !> LAI grid at level 0, with dim2 = time
96 real(dp), dimension(:, :), intent(in) :: gridded_lai0
97 !> land cover at level 0
98 integer(i4), dimension(:, :), intent(in) :: lcover0
99 !> Empirical quantiles of slope
100 real(dp), dimension(:), intent(in) :: slope_emp0
101 !> Cell ids at level 0
102 integer(i4), dimension(:), intent(in) :: id0
103 !> Upper row of hi res block
104 integer(i4), dimension(:), intent(in) :: upper_bound1
105 !> Lower row of hi res block
106 integer(i4), dimension(:), intent(in) :: lower_bound1
107 !> Left column of hi res block
108 integer(i4), dimension(:), intent(in) :: left_bound1
109 !> Right column of hi res block
110 integer(i4), dimension(:), intent(in) :: right_bound1
111 !> Number of L0 cells within a L1 cell
112 integer(i4), dimension(:), intent(in) :: n_subcells1
113 !> y0 at level 0
114 real(dp), dimension(:), intent(in) :: y0
115 !> [1] fraction of sealed area
116 real(dp), dimension(:, :, :), intent(inout) :: fsealed1
117 !> Parameter that determines the rel. contribution to SM, upscal. Bulk den.
118 real(dp), dimension(:, :, :), intent(inout) :: soilmoistexp1
119 !> [1] jarvis critical value for norm SWC
120 real(dp), dimension(:, :, :), intent(inout) :: jarvis_thresh_c1
121 !> [10^-3 m] depth of saturated SM
122 real(dp), dimension(:, :, :), intent(inout) :: soilmoistsat1
123 !> [10^-3 m] field capacity
124 real(dp), dimension(:, :, :), intent(inout) :: soilmoistfc1
125 !> [10^-3 m] permanent wilting point
126 real(dp), dimension(:, :, :), intent(inout) :: wiltingpoint1
127 !> fraction of roots in soil horizon
128 real(dp), dimension(:, :, :), intent(inout) :: froots1
129 !> [degreeC] threshold temperature for snow rain
130 real(dp), dimension(:, :, :), intent(inout) :: tempthresh1
131 !> [mm-1 degreeC-1] Degree-day factor with no precipitation
132 real(dp), dimension(:, :, :), intent(inout) :: degdaynopre1
133 !> [mm-1 degreeC-1] Maximum Degree-day factor
134 real(dp), dimension(:, :, :), intent(inout) :: degdaymax1
135 !> [d-1 degreeC-1] Increase of the Degree-day factor per mm of increase in precipitation
136 real(dp), dimension(:, :, :), intent(inout) :: degdayinc1
137 !> [1] PET correction for Aspect at level 1
138 real(dp), dimension(:, :, :), intent(inout) :: fasp1
139 !> [1] PET Hargreaves Samani coeff. at level 1
140 real(dp), dimension(:, :, :), intent(inout) :: harsamcoeff1
141 !> [1] PET Priestley Taylor coeff. at level 1
142 real(dp), dimension(:, :, :), intent(inout) :: prietayalpha1
143 !> [s m-1] PET aerodynamical resitance at level 1
144 real(dp), dimension(:, :, :), intent(inout) :: aeroresist1
145 !> [s m-1] PET bulk surface resitance at level 1
146 real(dp), dimension(:, :, :), intent(inout) :: surfresist1
147 !> threshold parameter
148 real(dp), dimension(:, :, :), intent(inout) :: sealedthresh1
149 !> [10^-3 m] Threshhold water depth in upper reservoir (for Runoff contribution)
150 real(dp), dimension(:, :, :), intent(inout) :: unsatthresh1
151 !> [10^-3 m] Recession coefficient of the upper reservoir, upper outlet
152 real(dp), dimension(:, :, :), intent(inout) :: kfastflow1
153 !> [10^-3 m] Recession coefficient of the upper reservoir, lower outlet
154 real(dp), dimension(:, :, :), intent(inout) :: kslowflow1
155 !> Level 1 baseflow recession
156 real(dp), dimension(:, :, :), intent(inout) :: kbaseflow1
157 !> [1] Exponent for the upper reservoir
158 real(dp), dimension(:, :, :), intent(inout) :: alpha1
159 !> [d-1] percolation coefficient
160 real(dp), dimension(:, :, :), intent(inout) :: kperco1
161 !> karstic percolation loss
162 real(dp), dimension(:, :, :), intent(inout) :: karstloss1
163 !> max interception
164 real(dp), dimension(:, :, :), intent(inout) :: maxinter1
165 !> pet cor factor at level-1
166 real(dp), dimension(:, :, :), intent(inout) :: petlaicorfactor
167 !> [-] inital neutron count
168 real(dp), dimension(:, :, :), intent(inout) :: no_count1
169 !> [gcm-3] bulk density
170 real(dp), dimension(:, :, :), intent(inout) :: bulkdens1
171 !> [mm/mm] lattice water content
172 real(dp), dimension(:, :, :), intent(inout) :: latticewater1
173 !> [-] cosmic L3 parameter
174 real(dp), dimension(:, :, :), intent(inout) :: cosmicl31
175
176 !> array of global parameters
177 real(dp), dimension(:), intent(in), optional, target :: parameterset
178
179 ! array of global parameters
180 real(dp), dimension(:), pointer :: param
181
182 real(dp), dimension(:, :, :), allocatable :: thetas_till
183
184 real(dp), dimension(:, :, :), allocatable :: thetafc_till
185
186 real(dp), dimension(:, :, :), allocatable :: thetapw_till
187
188 ! saturated hydraulic conductivity
189 real(dp), dimension(:, :, :), allocatable :: ks
190
191 ! Bulk density
192 real(dp), dimension(:, :, :), allocatable :: db
193 real(dp), dimension(:, :), allocatable :: thetas
194 real(dp), dimension(:, :), allocatable :: thetafc
195 real(dp), dimension(:, :), allocatable :: thetapw
196
197 ! neutron count
198 real(dp), dimension(:,:,:), allocatable :: latwat_till
199 real(dp), dimension(:,:,:), allocatable :: cosmic_l3_till
200 real(dp), dimension(:,:), allocatable :: latwat ! lattice water
201 real(dp), dimension(:,:), allocatable :: cosmic_l3 ! COSMIC parameter L3
202
203
204 ! relative variability of saturated
205 ! hydraulic cound. for Horizantal flow
206 real(dp), dimension(:), allocatable :: ksvar_h0
207
208 ! relative variability of saturated
209 ! hydraulic cound. for vertical flow
210 real(dp), dimension(:), allocatable :: ksvar_v0
211
212 ! soil mositure deficit from
213 ! field cap. w.r.t to saturation
214 real(dp), dimension(:), allocatable :: sms_fc0
215
216 ! L0 baseflow parameter
217 real(dp), dimension(size(Id0, 1)) :: k2_0
218
219 ! L1 baseflow parameter
220 real(dp), dimension(:), allocatable :: k2_1
221
222 ! L0 Aspect
223 real(dp), dimension(size(Id0, 1)) :: fasp0
224
225 ! number of soil classes
226 integer(i4) :: msoil
227
228 ! maximum of number of Tillage horizons
229 integer(i4) :: mtill
230
231 ! maximum number of horizons
232 integer(i4) :: mhor
233
234 ! number of Landcover classes
235 integer(i4) :: mlc
236
237 ! indexing of parameter vector - start
238 integer(i4) :: istart
239
240 ! indexing of parameter vector - end
241 integer(i4) :: iend
242
243 ! 2nd indexing of parameter vector - start
244 integer(i4) :: istart2
245
246 ! 2nd indexing of parameter vector - end
247 integer(i4) :: iend2
248
249 ! counter for looping over LCscenes
250 integer(i4) :: iilc
251
252 ! [1] Fraction of forest cover
253 real(dp), dimension(size(fSealed1, dim = 1)) :: fforest1
254
255 ! [1] Fraction of permeable cover
256 real(dp), dimension(size(fSealed1, dim = 1)) :: fperm1
257
258
259 if (present(parameterset)) then
260 param => parameterset
261 else
262 param => global_parameters(:, 3)
263 end if
264
265
266 ! loop over all LCover scenes
267 do iilc = 1, size(lcover0, 2)
268
269 ! estimate land cover fractions for dominant landcover class
270 ! fSealed is intent inout, the rest only intent in
271 fforest1(:) = l0_fractionalcover_in_lx(lcover0(:, iilc), 1, mask0, &
272 upper_bound1, &
273 lower_bound1, &
274 left_bound1, &
275 right_bound1, &
276 n_subcells1)
277 fsealed1(:, 1, iilc) = l0_fractionalcover_in_lx(lcover0(:, iilc), 2, mask0, &
278 upper_bound1, &
279 lower_bound1, &
280 left_bound1, &
281 right_bound1, &
282 n_subcells1)
283 fperm1(:) = l0_fractionalcover_in_lx(lcover0(:, iilc), 3, mask0, &
284 upper_bound1, &
285 lower_bound1, &
286 left_bound1, &
287 right_bound1, &
288 n_subcells1)
289 !---------------------------------------------------------
290 ! Update fractions of sealed area fractions
291 ! based on the sealing fraction[0-1] in cities
292 !---------------------------------------------------------
293 ! a factor is applied to the sealed area, effectively reducing it
294 fsealed1(:, 1, iilc) = fracsealed_cityarea * fsealed1(:, 1, iilc)
295 ! the forest area is kept constant, but the permeable area is increased so that the
296 ! sum off all fractions equals 1 again
297 fperm1(:) = 1.0_dp - fsealed1(:, 1, iilc) - fforest1(:)
298
299 ! ------------------------------------------------------------------
300 ! snow parameters
301 ! ------------------------------------------------------------------
302 select case(processmatrix(2,1))
303 case(1)
304
305 istart = processmatrix(2, 3) - processmatrix(2, 2) + 1
306 iend = processmatrix(2, 3)
307
308 call snow_acc_melt_param(param(istart : iend), & ! intent(in)
309 fforest1, fsealed1(:, 1, iilc), fperm1, & ! intent(in)
310 tempthresh1(:, 1, iilc), degdaynopre1(:, 1, iilc), & ! intent(out)
311 degdayinc1(:, 1, iilc), degdaymax1(:, 1, iilc) & ! intent(out)
312 )
313 case DEFAULT
314 call error_message('***ERROR: Process description for process "snow pack" does not exist! mo_multi_param_reg')
315 end select
316
317 ! ------------------------------------------------------------------
318 ! Soil moisture parametrization
319 ! ------------------------------------------------------------------
320 msoil = size(soildb%is_present, 1)
321 mlc = maxval(lcover0(:, iilc), (lcover0(:, iilc) .ne. nodata_i4))
322
323 ! depending on which kind of soil database processing is to be performed
324 if(iflag_soildb .eq. 0)then
325 mtill = maxval(soildb%nTillHorizons, (soildb%nTillHorizons .ne. nodata_i4))
326 mhor = maxval(soildb%nHorizons, (soildb%nHorizons .ne. nodata_i4))
327 else if(iflag_soildb .eq. 1) then
328 ! here for each soil type both till and non-till soil hydraulic properties are to be estimated
329 ! since a given soil type can lie in any horizon (till or non-till ones)
330 ! adopt it in a way that it do not break the consistency of iFlag_soilDB = 0
331 ! ** NOTE: SDB_nTillHorizons and SDB_nHorizons are also assigned in
332 ! this flag option (see mo_soildatabase.f90 file - read_soil_LUT).
333 ! But we are not using those variables here since in this case we have not
334 ! varying number of soil horizons or either tillage horizons.
335 ! So assigning them with a value = 1 is more than enough.
336 mtill = 1
337 mhor = 1
338 end if
339
340 allocate(thetas_till(msoil, mtill, mlc))
341 allocate(thetafc_till(msoil, mtill, mlc))
342 allocate(thetapw_till(msoil, mtill, mlc))
343 allocate(thetas(msoil, mhor))
344 allocate(thetafc(msoil, mhor))
345 allocate(thetapw(msoil, mhor))
346 allocate(ks(msoil, mhor, mlc))
347 allocate(db(msoil, mhor, mlc))
348
349 ! neutron count related ones
350 ! allocate and initalize here
351 allocate( latwat_till(msoil, mtill, mlc))
352 allocate(cosmic_l3_till(msoil, mtill, mlc))
353 allocate( latwat(msoil, mhor ))
354 allocate( cosmic_l3(msoil, mhor ))
355 latwat_till = 0.000001_dp
356 cosmic_l3_till = 0.000001_dp
357 cosmic_l3 = 0.000001_dp
358 latwat = 0.000001_dp
359
360
361
362 ! earlier these variables were allocated with size(soilId0,1)
363 ! in which the variable "soilId0" changes according to the iFlag_soilDB
364 ! so better to use other variable which is common to both soilDB (0 AND 1) flags
365 allocate(ksvar_h0(size(id0, 1)))
366 allocate(ksvar_v0(size(id0, 1)))
367 allocate( sms_fc0(size(id0, 1)))
368
369 select case(processmatrix(3,1))
370 case(1)
371 ! first thirteen parameters go to this routine
372 istart = processmatrix(3, 3) - processmatrix(3, 2) + 1
373 iend = processmatrix(3, 3) - 4
374
375 ! next four parameters go here
376 ! (the first three for the fRoots and the fourth one for the beta)
377 istart2 = processmatrix(3, 3) - 4 + 1
378 iend2 = processmatrix(3, 3)
379
380 case(2)
381 ! first thirteen parameters go to this routine
382 istart = processmatrix(3, 3) - processmatrix(3, 2) + 1
383 iend = processmatrix(3, 3) - 5
384
385 ! next four parameters go here
386 ! (the first three for the fRoots and the fourth one for the beta)
387 istart2 = processmatrix(3, 3) - 5 + 1
388 iend2 = processmatrix(3, 3) - 1
389 ! last parameter is jarvis parameter - no need to be regionalized
390 jarvis_thresh_c1 = param(processmatrix(3, 3))
391 !write(*,*) 'jarvis_thresh_c1 = ', jarvis_thresh_c1
392 !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
393
394 case(3)
395 ! first thirteen parameters go to this routine
396 istart = processmatrix(3, 3) - processmatrix(3, 2) + 1
397 iend = processmatrix(3, 3) - 9
398
399 ! next four parameters go here
400 ! (the first three for the fRoots and the fourth one for the beta)
401 istart2 = processmatrix(3, 3) - 8
402 iend2 = processmatrix(3, 3) - 1
403
404 ! last parameter is jarvis parameter - no need to be regionalized
405 jarvis_thresh_c1 = param(processmatrix(3, 3))
406
407 !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
408 !write(*,*) 'jarvis_thresh_c1 = ', jarvis_thresh_c1
409
410 case(4)
411 ! first thirteen parameters go to this routine
412 istart = processmatrix(3, 3) - processmatrix(3, 2) + 1
413 iend = processmatrix(3, 3) - 8
414
415 ! next four parameters go here
416 ! (the first three for the fRoots and the fourth one for the beta)
417 istart2 = processmatrix(3, 3) - 7
418 iend2 = processmatrix(3, 3)
419 !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
420
421 case DEFAULT
422 call error_message('***ERROR: Process description for process "soil moisture parametrization"', &
423 'does not exist! mo_multi_param_reg')
424 end select
425
426 call mpr_sm(param(istart : iend), processmatrix, &
427 soildb%is_present, soildb%nHorizons, soildb%nTillHorizons, &
428 soildb%sand, soildb%clay, soildb%DbM, &
429 id0, soilid0, lcover0(:, iilc), &
430 thetas_till, thetafc_till, thetapw_till, thetas, &
431 thetafc, thetapw, ks, db, ksvar_h0, ksvar_v0, sms_fc0)
432
433 ! >> neutron count related parameters
434 if ( processmatrix(10,1) .GT. 0 ) &
435 call mpr_neutrons( processmatrix(10,1), & ! IN: processmatrix case
436 param( processmatrix(10,3)-processmatrix(10,2)+1:processmatrix(10,3) ) , & ! IN: global parameter set
437 soildb%is_present , & ! IN: flag indicating presence of soil
438 soildb%nHorizons , & ! IN: Number of Horizons of Soiltype
439 soildb%nTillHorizons , & ! IN: Number of tillage Horizons
440 lcover0(:, iilc) , & ! IN: land cover ids at level 0
441 soildb%clay , & ! IN: clay content
442 soildb%DbM , & ! IN: mineral Bulk density
443 db , & ! IN: Bulk density
444 cosmic_l3_till , & ! OUT: COSMIC parameter L3 tillage layer
445 latwat_till , & ! OUT: COSMIC parameter Lattice Water tillage layer
446 cosmic_l3 , & ! OUT: COSMIC parameter L3
447 latwat & ! OUT: COSMIC parameter Lattice Water
448 )
449
450 call mpr_smhorizons(param(istart2:iend2), processmatrix, &
451 iflag_soildb, nsoilhorizons_mhm, horizondepth_mhm, &
452 lcover0(:, iilc), soilid0, &
453 soildb%nHorizons, soildb%nTillHorizons, &
454 thetas_till, thetafc_till, thetapw_till, &
455 thetas, thetafc, thetapw, &
456 soildb%Wd, db, soildb%DbM, soildb%RZdepth, &
457 mask0, id0, &
458 upper_bound1, lower_bound1, left_bound1, right_bound1, n_subcells1, &
459 soilmoistexp1(:, :, iilc), soilmoistsat1(:, :, iilc), soilmoistfc1(:, :, iilc), &
460 wiltingpoint1(:, :, iilc), froots1(:, :, iilc), &
461 ! >>>>>> neutron count
462 latwat_till, cosmic_l3_till, latwat, cosmic_l3, &
463 bulkdens1(:,:,iilc), latticewater1(:,:,iilc), cosmicl31(:,:,iilc) &
464 )
465
466 deallocate(thetas_till)
467 deallocate(thetafc_till)
468 deallocate(thetapw_till)
469 deallocate(thetas)
470 deallocate(thetafc)
471 deallocate(thetapw)
472 deallocate(ks)
473 deallocate(db)
474
475 ! neutron count
476 deallocate( latwat_till )
477 deallocate( cosmic_l3_till )
478 deallocate( latwat )
479 deallocate( cosmic_l3 )
480
481 ! ------------------------------------------------------------------
482 ! potential evapotranspiration (PET)
483 ! ------------------------------------------------------------------
484 ! Penman-Monteith method is only method that is LCscene dependent
485 if (processmatrix(5, 1) == 3) then
486 istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
487 iend = processmatrix(5, 3)
488 call aerodynamical_resistance(gridded_lai0, lcover0(:, iilc), param(istart : iend - 1), mask0, &
489 id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
490 aeroresist1(:, :, iilc))
491 else if (processmatrix(5, 1) == -1) then
492 istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
493 iend = processmatrix(5, 3)
494
495 call pet_correctbylai(param(istart : iend), nodata_dp, &
496 lcover0(:, iilc), gridded_lai0, mask0, id0, &
497 upper_bound1, lower_bound1, left_bound1, &
498 right_bound1, n_subcells1, petlaicorfactor(:, :, iilc))
499 end if
500
501 ! ------------------------------------------------------------------
502 ! interflow
503 ! ------------------------------------------------------------------
504 select case(processmatrix(6, 1))
505 case (1)
506 !
507 istart = processmatrix(6, 3) - processmatrix(6, 2) + 1
508 iend = processmatrix(6, 3)
509 ! TODO: this subroutine should be split into each param (or at least extract kFastFlow1)
510 ! because it is in the loop unnecessarily
511 call mpr_runoff(lcover0(:, iilc), mask0, sms_fc0, slope_emp0, &
512 ksvar_h0, param(istart : iend), id0, upper_bound1, lower_bound1, &
513 left_bound1, right_bound1, n_subcells1, unsatthresh1(:, 1, 1), kfastflow1(:, 1, iilc), &
514 kslowflow1(:, 1, iilc), alpha1(:, 1, iilc))
515 case DEFAULT
516 call error_message('***ERROR: Process description for process "interflow" does not exist! mo_multi_param_reg')
517 END select
518
519 ! ------------------------------------------------------------------
520 ! percolation cofficient, karstic percolation loss
521 ! ------------------------------------------------------------------
522 select case(processmatrix(7, 1))
523 case(1)
524
525 istart = processmatrix(7, 3) - processmatrix(7, 2) + 1
526 iend = processmatrix(7, 3)
527 call karstic_layer(& ! In
528 param(istart : iend), & ! In
529 geounit0, mask0, & ! In
530 sms_fc0, ksvar_v0, id0, & ! In
531 n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, & ! In
532 karstloss1(:, 1, 1), kperco1(:, 1, iilc) & ! Out
533 )
534
535 case DEFAULT
536 call error_message('***ERROR: Process description for process "percolation" does not exist! mo_multi_param_reg')
537 end select
538
539 deallocate(ksvar_h0)
540 deallocate(ksvar_v0)
541 deallocate(sms_fc0)
542
543 end do !! >>>>>>> LAND COVER SCENE LOOP
544
545
546 ! ------------------------------------------------------------------
547 ! sealed area threshold for runoff generation
548 ! ------------------------------------------------------------------
549 select case(processmatrix(4, 1))
550 case (1)
551 istart = processmatrix(4, 3) - processmatrix(4, 2) + 1
552 iend = processmatrix(4, 3)
553 call iper_thres_runoff(param(istart : iend), sealedthresh1)
554 case DEFAULT
555 call error_message('***ERROR: Process description for process "runoff_generation" does not exist! mo_multi_param_reg')
556 end select
557
558 ! ------------------------------------------------------------------
559 ! potential evapotranspiration (PET)
560 ! ------------------------------------------------------------------
561 select case(processmatrix(5, 1))
562 case(-1) ! LAI correction of input PET
563 iend = -9999 ! dummy statement
564 case(0) ! aspect correction of input PET
565 istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
566 iend = processmatrix(5, 3)
567 call pet_correctbyasp(id0, y0, asp0, param(istart : iend), nodata_dp, fasp0)
568 fasp1(:, 1, 1) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
569 left_bound1, right_bound1, id0, mask0, nodata_dp, fasp0)
570 case(1) ! Hargreaves-Samani method
571 istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
572 iend = processmatrix(5, 3)
573 call pet_correctbyasp(id0, y0, asp0, param(istart : iend - 1), nodata_dp, fasp0)
574 fasp1(:, 1, 1) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
575 left_bound1, right_bound1, id0, mask0, nodata_dp, fasp0)
576 harsamcoeff1 = param(iend)
577 case(2) ! Priestley-Taylor Method
578 istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
579 iend = processmatrix(5, 3)
580 call priestley_taylor_alpha(gridded_lai0, param(istart : iend), &
581 mask0, nodata_dp, id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
582 prietayalpha1(:, :, 1))
583 case(3) ! Penman-Monteith method
584 ! aerodynamic resistance is calculated inside LCscene loop
585 istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
586 iend = processmatrix(5, 3)
587 call bulksurface_resistance(gridded_lai0, param(iend), mask0, &
588 nodata_dp, id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
589 surfresist1(:, :, 1))
590 case default
591 call error_message('***ERROR: Process description for process "pet correction" does not exist! mo_multi_param_reg')
592 end select
593 ! ------------------------------------------------------------------
594 ! baseflow recession parameter
595 ! ------------------------------------------------------------------
596 select case(processmatrix(9, 1))
597 case(1)
598 ! the number of process parameters, so the number in processMatrix(9,2) has
599 ! to be equal to the size of geo_unit_list
600 istart = processmatrix(9, 3) - processmatrix(9, 2) + 1
601 iend = processmatrix(9, 3)
602
603 call baseflow_param(param(istart : iend), geounit0, k2_0)
604
605 ! Upscale by arithmetic mean
606 allocate(k2_1(size(kbaseflow1, 1)))
607 k2_1 = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
608 left_bound1, right_bound1, id0, mask0, nodata_dp, k2_0)
609 ! loop over all LCover scenes
610 do iilc = 1, size(lcover0, 2)
611 kbaseflow1(:, 1, iilc) = k2_1
612 end do
613 deallocate(k2_1)
614
615 ! correction and unit conversion
616 ! if percolation is ON: correct K2 such that it is at least k1
617 ! since kSlowFlow1 is LCover dependent, kBaseFlow1 is too
618 if (processmatrix(7, 1) .gt. 0) kbaseflow1 = merge(kslowflow1, kbaseflow1, kbaseflow1 .lt. kslowflow1)
619
620 case DEFAULT
621 call error_message('***ERROR: Process description for process "baseflow Recession" does not exist! mo_multi_param_reg')
622 end select
623
624 ! ------------------------------------------------------------------
625 ! Neutron count related parameters
626 ! >> only N0 parameter - others are defined above in soil parameters
627 ! ------------------------------------------------------------------
628 select case(processmatrix(10, 1))
629 case(0)
630 ! do nothing
631 case(1)
632 ! the number of process parameters, so the number in processMatrix(9,2) has
633 istart = processmatrix(10, 3) - processmatrix(10, 2) + 1
634 iend = processmatrix(10, 3)
635 no_count1 = param(istart) ! >> 1st parameter --> N0 parameter
636 case(2)
637 ! the number of process parameters, so the number in processMatrix(9,2) has
638 istart = processmatrix(10, 3) - processmatrix(10, 2) + 1
639 iend = processmatrix(10, 3)
640 no_count1 = param(istart) ! >> 1st parameter --> N0 parameter
641 case DEFAULT
642 call error_message('***ERROR: Process description for process "Neutron count" does not exist! mo_multi_param_reg')
643 end select
644
645
646 !-------------------------------------------------------------------
647 ! call regionalization of parameters related to LAI
648 ! it is now outside of mHM since LAI is now dynamic variable
649 !-------------------------------------------------------------------
650 call canopy_intercept_param(processmatrix, param(:), &
651 gridded_lai0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, id0, mask0, &
652 nodata_dp, maxinter1(:, :, 1))
653
654 end subroutine mpr
655
656 ! ----------------------------------------------------------------------------
657
658 ! NAME
659 ! baseflow_param
660
661 ! PURPOSE
662 !> \brief baseflow recession parameter
663
664 !> \details This subroutine calculates the baseflow recession parameter
665 !> based on the geological units at the Level 0 scale. For each level 0
666 !> cell, it assigns the value specified in the parameter array param for the
667 !> geological unit in this cell.
668 !> Global parameters needed (see mhm_parameter.nml):
669 !> - param(1) = GeoParam(1,:)
670 !> - param(2) = GeoParam(2,:)
671 !> - ...
672
673 ! INTENT(IN)
674 !> \param[in] "real(dp), dimension(:) :: param" list of required parameters
675 !> \param[in] "integer(i4), dimension(:) :: geoUnit0" ids of geological units at L0
676
677 ! INTENT(OUT)
678 !> \param[out] "real(dp), dimension(:) :: k2_0" - baseflow recession parameter at Level 0
679
680 ! HISTORY
681 !> \authors Stephan Thober, Rohini Kumar
682
683 !> \date Dec 2012
684
685 ! Modifications:
686 ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
687 ! Robert Schweppe Jun 2018 - refactoring and reformatting
688
689 subroutine baseflow_param(param, geoUnit0, k2_0)
690
691 use mo_mpr_global_variables, only : geounitlist
692 !$ use omp_lib
693
694 implicit none
695
696 ! list of required parameters
697 real(dp), dimension(:), intent(in) :: param
698
699 ! ids of geological units at L0
700 integer(i4), dimension(:), intent(in) :: geoUnit0
701
702 ! - baseflow recession parameter at Level 0
703 real(dp), dimension(:), intent(out) :: k2_0
704
705 ! loop variable
706 integer(i4) :: ii
707
708 ! geo unit
709 integer(i4), dimension(1) :: gg
710
711
712 if (size(param) .ne. size(geounitlist)) &
713 call error_message(' mo_multi_param_reg: baseflow_param: size mismatch, subroutine baseflow parameters ')
714
715 k2_0 = nodata_dp
716
717 !$OMP PARALLEL
718 !$OMP DO PRIVATE(gg) SCHEDULE(STATIC)
719 do ii = 1, size(k2_0)
720 ! get parameter index in geoUnitList
721 gg = minloc(abs(geounitlist - geounit0(ii)))
722 k2_0(ii) = param(gg(1))
723 end do
724 !$OMP END DO
725 !$OMP END PARALLEL
726
727 end subroutine baseflow_param
728
729 ! ----------------------------------------------------------------------------
730
731 ! NAME
732 ! snow_acc_melt_param
733
734 ! PURPOSE
735 !> \brief Calculates the snow parameters.
736
737 !> \details This subroutine calculates the snow parameters
738 !> threshold temperature (TT), degree-day factor without precipitation (DD)
739 !> and maximum degree-day factor (DDmax) as well as increase of degree-day
740 !> factor per mm of increase in precipitation (IDDP).
741
742 !> Global parameters needed (see mhm_parameter.nml):
743 !> - param(1) = snowTreshholdTemperature
744 !> - param(2) = degreeDayFactor_forest
745 !> - param(3) = degreeDayFactor_impervious
746 !> - param(4) = degreeDayFactor_pervious
747 !> - param(5) = increaseDegreeDayFactorByPrecip
748 !> - param(6) = maxDegreeDayFactor_forest
749 !> - param(7) = maxDegreeDayFactor_impervious
750 !> - param(8) = maxDegreeDayFactor_pervious
751 !> INTENT(IN)
752
753 ! INTENT(IN)
754 !> \param[in] "real(dp), dimension(8) :: param" eight global parameters
755 !> \param[in] "real(dp), dimension(:) :: fForest1" [1] fraction of forest cover
756 !> \param[in] "real(dp), dimension(:) :: fIperm1" [1] fraction of sealed area
757 !> \param[in] "real(dp), dimension(:) :: fPerm1" [1] fraction of permeable area
758
759 ! INTENT(OUT)
760 !> \param[out] "real(dp), dimension(:) :: tempThresh1" [degreeC] threshold temperature for snow rain
761 !> \param[out] "real(dp), dimension(:) :: degDayNoPre1" [mm-1 degreeC-1] Degree-day factor with
762 !> \param[out] "real(dp), dimension(:) :: degDayInc1" [d-1 degreeC-1] Increase of the Degree-day
763 !> \param[out] "real(dp), dimension(:) :: degDayMax1" [mm-1 degreeC-1] Maximum Degree-day factor
764
765 ! HISTORY
766 !> \authors Stephan Thober, Rohini Kumar
767
768 !> \date Dec 2012
769
770 ! Modifications:
771 ! Juliane Mai Oct 2013 - OLD parametrization --> param(1) = snowTreshholdTemperature
772 ! --> param(2) = degreeDayFactor_forest
773 ! --> param(3) = degreeDayFactor_impervious
774 ! --> param(4) = degreeDayFactor_pervious
775 ! --> param(5) = increaseDegreeDayFactorByPrecip
776 ! --> param(6) = maxDegreeDayFactor_forest
777 ! --> param(7) = maxDegreeDayFactor_impervious
778 ! --> param(8) = maxDegreeDayFactor_pervious
779 ! -------------------------------
780 ! degreeDayFactor_impervious = degreeDayFactor_forest + delta_1 + delta_2
781 ! degreeDayFactor_pervious = degreeDayFactor_forest + delta_1
782 ! maxDegreeDayFactor_forest = degreeDayFactor_forest + delta_3
783 ! maxDegreeDayFactor_impervious = degreeDayFactor_impervious + delta_5
784 ! = degreeDayFactor_forest + delta_1 + delta_2 + delta_5
785 ! maxDegreeDayFactor_pervious = degreeDayFactor_pervious + delta_4
786 ! = degreeDayFactor_forest + delta_1 + delta_4
787 ! -------------------------------
788 ! NEW parametrization
789 ! --> param(1) = snowTreshholdTemperature
790 ! --> param(2) = degreeDayFactor_forest
791 ! --> param(3) = delta_2
792 ! --> param(4) = delta_1
793 ! --> param(5) = increaseDegreeDayFactorByPrecip
794 ! --> param(6) = delta_3
795 ! --> param(7) = delta_5
796 ! --> param(8) = delta_4
797 ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
798 ! Robert Schweppe Jun 2018 - refactoring and reformatting
799
800 subroutine snow_acc_melt_param(param, fForest1, fIperm1, fPerm1, tempThresh1, degDayNoPre1, degDayInc1, &
801 degDayMax1)
802 implicit none
803
804 ! eight global parameters
805 real(dp), dimension(8), intent(in) :: param
806
807 ! [1] fraction of forest cover
808 real(dp), dimension(:), intent(in) :: fForest1
809
810 ! [1] fraction of sealed area
811 real(dp), dimension(:), intent(in) :: fIperm1
812
813 ! [1] fraction of permeable area
814 real(dp), dimension(:), intent(in) :: fPerm1
815
816 ! [degreeC] threshold temperature for snow rain
817 real(dp), dimension(:), intent(out) :: tempThresh1
818
819 ! [mm-1 degreeC-1] Degree-day factor with
820 real(dp), dimension(:), intent(out) :: degDayNoPre1
821
822 ! [mm-1 degreeC-1] Maximum Degree-day factor
823 real(dp), dimension(:), intent(out) :: degDayMax1
824
825 ! [d-1 degreeC-1] Increase of the Degree-day
826 real(dp), dimension(:), intent(out) :: degDayInc1
827
828 real(dp) :: tmp_degreeDayFactor_forest, tmp_degreeDayFactor_impervious, tmp_degreeDayFactor_pervious
829
830 real(dp) :: tmp_maxDegreeDayFactor_forest, tmp_maxDegreeDayFactor_impervious, tmp_maxDegreeDayFactor_pervious
831
832
833 tmp_degreedayfactor_forest = param(2) ! OLD: param(2)
834 tmp_degreedayfactor_impervious = param(2) + param(4) + param(3) ! OLD: param(3)
835 tmp_degreedayfactor_pervious = param(2) + param(4) ! OLD: param(4)
836 tmp_maxdegreedayfactor_forest = param(2) + param(6) ! OLD: param(6)
837 tmp_maxdegreedayfactor_impervious = param(2) + param(4) + param(3) + param(7) ! OLD: param(7)
838 tmp_maxdegreedayfactor_pervious = param(2) + param(4) + param(8) ! OLD: param(8)
839
840 tempthresh1 = param(1)
841 degdayinc1 = param(5)
842
843 degdaynopre1 = (&
844 tmp_degreedayfactor_forest * fforest1 + &
845 tmp_degreedayfactor_impervious * fiperm1 + &
846 tmp_degreedayfactor_pervious * fperm1)
847 degdaymax1 = (&
848 tmp_maxdegreedayfactor_forest * fforest1 + &
849 tmp_maxdegreedayfactor_impervious * fiperm1 + &
850 tmp_maxdegreedayfactor_pervious * fperm1)
851
852 end subroutine snow_acc_melt_param
853
854 ! ----------------------------------------------------------------------------
855
856 ! NAME
857 ! iper_thres_runoff
858
859 ! PURPOSE
860 !> \brief sets the impervious layer threshold parameter for runoff generation
861
862 !> \details to be done by Kumar
863 !> ....
864
865 !> Global parameters needed (see mhm_parameter.nml):
866 !> - param(1) = imperviousStorageCapacity
867
868 ! INTENT(IN)
869 !> \param[in] "real(dp), dimension(1) :: param" - given threshold parameter
870
871 ! INTENT(OUT)
872 !> \param[out] "real(dp), dimension(:, :, :) :: sealedThresh1"
873
874 ! HISTORY
875 !> \authors Stephan Thober, Rohini Kumar
876
877 !> \date Dec 2012
878
879 ! Modifications:
880 ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
881 ! Robert Schweppe Jun 2018 - refactoring and reformatting
882
883 subroutine iper_thres_runoff(param, sealedThresh1)
884 implicit none
885
886 ! - given threshold parameter
887 real(dp), dimension(1), intent(in) :: param
888
889 real(dp), dimension(:, :, :), intent(out) :: sealedThresh1
890
891
892 sealedthresh1 = param(1)
893
894 end subroutine iper_thres_runoff
895
896 ! ----------------------------------------------------------------------------
897
898 ! NAME
899 ! karstic_layer
900
901 ! PURPOSE
902 !> \brief calculates the Karstic percolation loss
903
904 !> \details This subroutine calls first the karstic_fraction upscaling
905 !> routine for determine the karstic fraction area for every Level 1
906 !> cell. Then, the karstic percolation loss is estimated given two
907 !> shape parameters by
908 !> \f[ karstLoss1 = 1 + ( fKarArea * param(1)) *( (-1)**INT(param(2),i4) ) \f]
909 !> where \f$ karstLoss1 \f$ is the karstic percolation loss and \f$ fKarArea \f$
910 !> is the fraction of karstic area at level 1
911 !> Global parameters needed (see mhm_parameter.nml):
912 !> - param(1) = rechargeCoefficient
913 !> - param(2) = rechargeFactor_karstic
914 !> - param(3) = gain_loss_GWreservoir_karstic
915
916 ! INTENT(IN)
917 !> \param[in] "real(dp), dimension(3) :: param" parameters
918 !> \param[in] "integer(i4), dimension(:) :: geoUnit0" id of the Karstic formation
919 !> \param[in] "logical, dimension(:, :) :: mask0" mask at level 0
920 !> \param[in] "real(dp), dimension(:) :: SMs_FC0" [-] soil mositure deficit from field
921 !> capacity w.r.t to saturation
922 !> \param[in] "real(dp), dimension(:) :: KsVar_V0" [-] relative variability of saturated
923 !> \param[in] "integer(i4), dimension(:) :: Id0" Cell ids of hi res field
924 !> \param[in] "integer(i4), dimension(:) :: n_subcells1" number of l0 cells within a l1 cell
925 !> \param[in] "integer(i4), dimension(:) :: upper_bound1" upper row of a l1 cell in l0 grid
926 !> \param[in] "integer(i4), dimension(:) :: lower_bound1" lower row of a l1 cell in l0 grid
927 !> \param[in] "integer(i4), dimension(:) :: left_bound1" left col of a l1 cell in l0 grid
928 !> \param[in] "integer(i4), dimension(:) :: right_bound1" right col of a l1 cell in l0 grid
929
930 ! INTENT(OUT)
931 !> \param[out] "real(dp), dimension(:) :: karstLoss1" [-] Karstic percolation loss
932 !> \param[out] "real(dp), dimension(:) :: L1_Kp" [d-1] percolation coefficient
933
934 ! HISTORY
935 !> \authors Rohini Kumar, Stephan Thober
936
937 !> \date Feb 2013
938
939 ! Modifications:
940 ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
941 ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
942 ! Robert Schweppe Jun 2018 - refactoring and reformatting
943
944 subroutine karstic_layer(param, geoUnit0, mask0, SMs_FC0, KsVar_V0, Id0, n_subcells1, upper_bound1, lower_bound1, &
945 left_bound1, right_bound1, karstLoss1, L1_Kp)
946
947 use mo_mpr_global_variables, only : geounitlist, geounitkar
948 use mo_upscaling_operators, only : l0_fractionalcover_in_lx, upscale_arithmetic_mean
949 !$ use omp_lib
950
951 implicit none
952
953 ! parameters
954 real(dp), dimension(3), intent(in) :: param
955
956 ! id of the Karstic formation
957 integer(i4), dimension(:), intent(in) :: geoUnit0
958
959 ! mask at level 0
960 logical, dimension(:, :), intent(in) :: mask0
961
962 ! [-] soil mositure deficit from field
963 ! capacity w.r.t to saturation
964 real(dp), dimension(:), intent(in) :: SMs_FC0
965
966 ! [-] relative variability of saturated
967 real(dp), dimension(:), intent(in) :: KsVar_V0
968
969 ! Cell ids of hi res field
970 integer(i4), dimension(:), intent(in) :: Id0
971
972 ! number of l0 cells within a l1 cell
973 integer(i4), dimension(:), intent(in) :: n_subcells1
974
975 ! upper row of a l1 cell in l0 grid
976 integer(i4), dimension(:), intent(in) :: upper_bound1
977
978 ! lower row of a l1 cell in l0 grid
979 integer(i4), dimension(:), intent(in) :: lower_bound1
980
981 ! left col of a l1 cell in l0 grid
982 integer(i4), dimension(:), intent(in) :: left_bound1
983
984 ! right col of a l1 cell in l0 grid
985 integer(i4), dimension(:), intent(in) :: right_bound1
986
987 ! [-] Karstic percolation loss
988 real(dp), dimension(:), intent(out) :: karstLoss1
989
990 ! [d-1] percolation coefficient
991 real(dp), dimension(:), intent(out) :: L1_Kp
992
993 ! fraction of karstic area
994 real(dp), dimension(:), allocatable :: fKarArea
995
996 ! temporal variable
997 real(dp), dimension(size(SMs_FC0, 1)) :: tmp
998
999 integer(i4) :: nGeoUnits
1000
1001 integer(i4) :: i
1002
1003
1004 ! ------------------------------------------------------------------
1005 ! PERCOLATION; 1/Kp = f(Ks)
1006 ! ------------------------------------------------------------------
1007 ! Regionalise Kp with variability of last soil layer property
1008 !$OMP PARALLEL
1009 tmp = merge(param(1) * (1.0_dp + sms_fc0) / (1.0_dp + ksvar_v0), &
1010 nodata_dp, id0 .ne. nodata_i4)
1011 !$OMP END PARALLEL
1012
1013 l1_kp = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
1014 left_bound1, right_bound1, id0, mask0, nodata_dp, tmp)
1015
1016 ! minimum constrains
1017 l1_kp = merge(2.0_dp, l1_kp, l1_kp .lt. 2.0_dp)
1018
1019 ngeounits = size(geounitlist, 1)
1020
1021 ! 1st calculate fraction of Karstic area
1022 allocate(fkararea(size(karstloss1, 1)))
1023 fkararea = 0.0_dp
1024
1025 do i = 1, ngeounits
1026 if(geounitkar(i) .eq. 0) cycle
1027 fkararea(:) = l0_fractionalcover_in_lx(geounit0, geounitlist(i), mask0, &
1028 upper_bound1, lower_bound1, left_bound1, right_bound1, n_subcells1)
1029 end do
1030
1031 ! 2nd calculate karstLoss1
1032 karstloss1 = 1.0_dp - (fkararea * param(2))
1033
1034 deallocate(fkararea)
1035
1036 end subroutine karstic_layer
1037
1038 ! ----------------------------------------------------------------------------
1039
1040 ! NAME
1041 ! canopy_intercept_param
1042
1043 ! PURPOSE
1044 !> \brief estimate effective maximum interception capacity at L1
1045
1046 !> \details estimate effective maximum interception capacity at L1 for a given
1047 !> Leaf Area Index field.
1048 !> Global parameters needed (see mhm_parameter.nml):
1049 !> Process Case 1:
1050 !> - param(1) = canopyInterceptionFactor
1051
1052 ! INTENT(IN)
1053 !> \param[in] "integer(i4), dimension(:, :) :: processMatrix" indicate processes
1054 !> \param[in] "real(dp), dimension(:) :: param" array of global parameters
1055 !> \param[in] "real(dp), dimension(:, :) :: LAI0" LAI at level-0(nCells0, time)
1056 !> \param[in] "integer(i4), dimension(:) :: n_subcells1" Number of L0 cells within a L1 cell
1057 !> \param[in] "integer(i4), dimension(:) :: upper_bound1" Upper row of high resolution block
1058 !> \param[in] "integer(i4), dimension(:) :: lower_bound1" Lower row of high resolution block
1059 !> \param[in] "integer(i4), dimension(:) :: left_bound1" Left column of high resolution block
1060 !> \param[in] "integer(i4), dimension(:) :: right_bound1" Right column of high resolution block
1061 !> \param[in] "integer(i4), dimension(:) :: Id0" Cell ids at level 0
1062 !> \param[in] "logical, dimension(:, :) :: mask0" mask at level 0 field
1063 !> \param[in] "real(dp) :: nodata" - nodata value
1064
1065 ! INTENT(OUT)
1066 !> \param[out] "real(dp), dimension(:, :) :: max_intercept1" max interception at level-1(nCells1, time)
1067
1068 ! HISTORY
1069 !> \authors Rohini Kumar
1070
1071 !> \date Aug. 2013
1072
1073 ! Modifications:
1074 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1075
1076 subroutine canopy_intercept_param(processMatrix, param, LAI0, n_subcells1, upper_bound1, lower_bound1, left_bound1, &
1077 right_bound1, Id0, mask0, nodata, max_intercept1)
1078
1079 use mo_string_utils, only : num2str
1080 use mo_upscaling_operators, only : upscale_arithmetic_mean
1081
1082 implicit none
1083
1084 ! indicate processes
1085 integer(i4), dimension(:, :), intent(in) :: processmatrix
1086
1087 ! array of global parameters
1088 real(dp), dimension(:), intent(in) :: param
1089
1090 ! LAI at level-0(nCells0, time)
1091 real(dp), dimension(:, :), intent(in) :: lai0
1092
1093 ! Number of L0 cells within a L1 cell
1094 integer(i4), dimension(:), intent(in) :: n_subcells1
1095
1096 ! Upper row of high resolution block
1097 integer(i4), dimension(:), intent(in) :: upper_bound1
1098
1099 ! Lower row of high resolution block
1100 integer(i4), dimension(:), intent(in) :: lower_bound1
1101
1102 ! Left column of high resolution block
1103 integer(i4), dimension(:), intent(in) :: left_bound1
1104
1105 ! Right column of high resolution block
1106 integer(i4), dimension(:), intent(in) :: right_bound1
1107
1108 ! Cell ids at level 0
1109 integer(i4), dimension(:), intent(in) :: id0
1110
1111 ! mask at level 0 field
1112 logical, dimension(:, :), intent(in) :: mask0
1113
1114 ! - nodata value
1115 real(dp), intent(in) :: nodata
1116
1117 ! max interception at level-1(nCells1, time)
1118 real(dp), dimension(:, :), intent(out) :: max_intercept1
1119
1120 integer(i4) :: istart, iend, it
1121
1122 real(dp), dimension(:), allocatable :: max_intercept0
1123
1124 real(dp), dimension(:), allocatable :: gamma_intercept
1125
1126
1127 ! ------------------------------------------------------------------
1128 ! Maximum interception parameter
1129 ! ------------------------------------------------------------------
1130 select case(processmatrix(1, 1))
1131 case(1)
1132 istart = processmatrix(1, 3) - processmatrix(1, 2) + 1
1133 iend = processmatrix(1, 3)
1134
1135 ! allocate space
1136 allocate(gamma_intercept(iend - istart + 1))
1137 allocate(max_intercept0(size(id0, 1)))
1138
1139 ! estimate max. intercept at Level-0
1140 gamma_intercept(:) = param(istart : iend)
1141
1142 do it = 1, size(lai0, 2)
1143 !$OMP PARALLEL
1144 max_intercept0(:) = lai0(:, it) * gamma_intercept(1)
1145 !$OMP END PARALLEL
1146
1147 ! Upscale by arithmetic mean
1148 max_intercept1(:, it) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, left_bound1, &
1149 right_bound1, id0, mask0, nodata, max_intercept0(:))
1150
1151 end do
1152
1153 deallocate(gamma_intercept)
1154 deallocate(max_intercept0)
1155 CASE DEFAULT
1156 call error_message('mo_multi_param_reg: This processMatrix=', num2str(processmatrix(1, 1)), ' is not implemented!')
1157 end select
1158
1159 end subroutine canopy_intercept_param
1160
1161
1162 ! ----------------------------------------------------------------------------
1163
1164 ! NAME
1165 ! aerodynamical_resistance
1166
1167 ! PURPOSE
1168 !> \brief Regionalization of aerodynamic resistance
1169
1170 !> \details estimation of aerodynamical resistance
1171 !> Global parameters needed (see mhm_parameter.nml):
1172 !> - param(1) = canopyheigth_forest
1173 !> - param(2) = canopyheigth_impervious
1174 !> - param(3) = canopyheigth_pervious
1175 !> - param(4) = displacementheight_coeff
1176 !> - param(5) = roughnesslength_momentum_coeff
1177 !> - param(6) = roughnesslength_heat_coeff
1178
1179 ! INTENT(IN)
1180 !> \param[in] "real(dp), dimension(:, :) :: LAI0" LAI at level-0
1181 !> \param[in] "integer(i4), dimension(:) :: LCover0" land cover field
1182 !> \param[in] "real(dp), dimension(6) :: param" input parameter
1183 !> \param[in] "logical, dimension(:, :) :: mask0" mask at level 0
1184 !> \param[in] "integer(i4), dimension(:) :: Id0" Cell ids of hi res field
1185 !> \param[in] "integer(i4), dimension(:) :: n_subcells1" number of l0 cells within a l1 cell
1186 !> \param[in] "integer(i4), dimension(:) :: upper_bound1" upper row of a l1 cell in l0 grid
1187 !> \param[in] "integer(i4), dimension(:) :: lower_bound1" lower row of a l1 cell in l0 grid
1188 !> \param[in] "integer(i4), dimension(:) :: left_bound1" left col of a l1 cell in l0 grid
1189 !> \param[in] "integer(i4), dimension(:) :: right_bound1" right col of a l1 cell in l0 grid
1190
1191 ! INTENT(OUT)
1192 !> \param[out] "real(dp), dimension(:, :) :: aerodyn_resistance1" aerodynmaical resistance
1193
1194 ! HISTORY
1195 !> \authors Matthias Zink
1196
1197 !> \date Apr 2013
1198
1199 ! Modifications:
1200 ! Matthias Zink Jun 2017 - moved from mo_multi_scale_param_reg.f90 to mo_mpr_pet.f90
1201 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1202
1203 subroutine aerodynamical_resistance(LAI0, LCover0, param, mask0, Id0, n_subcells1, upper_bound1, lower_bound1, &
1204 left_bound1, right_bound1, aerodyn_resistance1)
1205
1206 use mo_common_constants, only : eps_dp
1207 use mo_mpr_constants, only : windmeasheight, karman
1208 use mo_upscaling_operators, only : upscale_arithmetic_mean
1209
1210 implicit none
1211
1212 ! LAI at level-0
1213 real(dp), dimension(:, :), intent(in) :: LAI0
1214
1215 ! land cover field
1216 integer(i4), dimension(:), intent(in) :: LCover0
1217
1218 ! input parameter
1219 real(dp), dimension(6), intent(in) :: param
1220
1221 ! mask at level 0
1222 logical, dimension(:, :), intent(in) :: mask0
1223
1224 ! Cell ids of hi res field
1225 integer(i4), dimension(:), intent(in) :: Id0
1226
1227 ! number of l0 cells within a l1 cell
1228 integer(i4), dimension(:), intent(in) :: n_subcells1
1229
1230 ! upper row of a l1 cell in l0 grid
1231 integer(i4), dimension(:), intent(in) :: upper_bound1
1232
1233 ! lower row of a l1 cell in l0 grid
1234 integer(i4), dimension(:), intent(in) :: lower_bound1
1235
1236 ! left col of a l1 cell in l0 grid
1237 integer(i4), dimension(:), intent(in) :: left_bound1
1238
1239 ! right col of a l1 cell in l0 grid
1240 integer(i4), dimension(:), intent(in) :: right_bound1
1241
1242 ! aerodynmaical resistance
1243 real(dp), dimension(:, :), intent(out) :: aerodyn_resistance1
1244
1245 integer(i4) :: tt
1246
1247 real(dp), dimension(:), allocatable :: maxLAI
1248
1249 real(dp), dimension(:), allocatable :: zm
1250
1251 real(dp), dimension(:), allocatable :: canopy_height0
1252
1253 real(dp), dimension(:), allocatable :: zm_zero, zh_zero, displace
1254
1255 ! dim 1 = number of cells on level 0,
1256 ! dim2=month of year
1257 real(dp), dimension(:, :), allocatable :: aerodyn_resistance0
1258
1259
1260 ! initialize some things
1261 allocate(zm(size(lcover0, dim = 1))) ; zm = nodata_dp
1262 allocate(zm_zero(size(lcover0, dim = 1))) ; zm_zero = nodata_dp
1263 allocate(zh_zero(size(lcover0, dim = 1))) ; zh_zero = nodata_dp
1264 allocate(displace(size(lcover0, dim = 1))) ; displace = nodata_dp
1265 allocate(canopy_height0(size(lcover0, dim = 1))) ; canopy_height0 = nodata_dp
1266 allocate(aerodyn_resistance0(size(lcover0, dim = 1), size(lai0, 2))) ; aerodyn_resistance0 = nodata_dp
1267 allocate(maxlai(size(lcover0, dim = 1))) ; maxlai = nodata_dp
1268
1269 ! regionalization of canopy height
1270 ! substitute with canopy height
1271 canopy_height0 = merge(param(1), canopy_height0, lcover0 == 1) ! forest
1272 canopy_height0 = merge(param(2), canopy_height0, lcover0 == 2) ! impervious
1273
1274 ! old implementation used values from LUT statically for all cells (Jan-Dec, 12 values):
1275 ! 7 Intensive-orchards: 2.0, 2.0, 2.0, 2.0, 3.0, 3.5, 4.0, 4.0, 4.0, 2.5, 2.0, 2.0
1276 maxlai = maxval(lai0, dim=2)
1277
1278 do tt = 1, size(lai0, 2)
1279
1280 ! pervious canopy height is scaled with LAI
1281 canopy_height0 = merge((param(3) * lai0(:, tt) / maxlai), canopy_height0, lcover0 == 3) ! pervious
1282
1283 ! estimation of the aerodynamic resistance on the lower level
1284 ! see FAO Irrigation and Drainage Paper No. 56 (p. 19 ff) for more information
1285 zm = windmeasheight
1286 ! correction: if wind measurement height is below canopy height loagarithm becomes negative
1287 zm = merge(canopy_height0 + zm, zm, ((abs(zm - nodata_dp) .GT. eps_dp) .AND. (zm .LT. canopy_height0)))
1288
1289 ! zh = zm
1290 displace = param(4) * canopy_height0
1291 zm_zero = param(5) * canopy_height0
1292 zh_zero = param(6) * zm_zero
1293 !
1294 ! calculate aerodynamic resistance (changes monthly)
1295 aerodyn_resistance0(:, tt) = log((zm - displace) / zm_zero) * log((zm - displace) / zh_zero) / (karman**2.0_dp)
1296 aerodyn_resistance1(:, tt) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
1297 left_bound1, right_bound1, id0, mask0, nodata_dp, aerodyn_resistance0(:, tt))
1298
1299 end do
1300
1301 end subroutine aerodynamical_resistance
1302
1303END MODULE mo_multi_param_reg
mo_common_constants
Provides constants commonly used by mHM, mRM and MPR.
Definition mo_common_constants.f90:15
mo_common_constants::eps_dp
real(dp), parameter, public eps_dp
epsilon(1.0) in double precision
Definition mo_common_constants.f90:25
mo_common_constants::nodata_dp
real(dp), parameter, public nodata_dp
Definition mo_common_constants.f90:31
mo_common_constants::nodata_i4
integer(i4), parameter, public nodata_i4
Definition mo_common_constants.f90:30
mo_common_variables
Provides structures needed by mHM, mRM and/or mpr.
Definition mo_common_variables.F90:27
mo_common_variables::global_parameters
real(dp), dimension(:, :), allocatable, target, public global_parameters
Definition mo_common_variables.F90:135
mo_common_variables::processmatrix
integer(i4), dimension(nprocesses, 3), public processmatrix
Definition mo_common_variables.F90:126
mo_mpr_constants
Provides MPR specific constants.
Definition mo_mpr_constants.f90:12
mo_mpr_constants::windmeasheight
real(dp), parameter, public windmeasheight
assumed meteorol.
Definition mo_mpr_constants.f90:65
mo_mpr_constants::karman
real(dp), parameter, public karman
von karman constant
Definition mo_mpr_constants.f90:67
mo_mpr_global_variables
Global variables for mpr only.
Definition mo_mpr_global_variables.f90:20
mo_mpr_global_variables::geounitkar
integer(i4), dimension(:), allocatable, public geounitkar
Definition mo_mpr_global_variables.f90:82
mo_mpr_global_variables::geounitlist
integer(i4), dimension(:), allocatable, public geounitlist
Definition mo_mpr_global_variables.f90:81
mo_mpr_global_variables::nsoilhorizons_mhm
integer(i4), public nsoilhorizons_mhm
Definition mo_mpr_global_variables.f90:36
mo_mpr_global_variables::soildb
type(soiltype), public soildb
Definition mo_mpr_global_variables.f90:75
mo_mpr_global_variables::fracsealed_cityarea
real(dp), public fracsealed_cityarea
Definition mo_mpr_global_variables.f90:98
mo_mpr_global_variables::iflag_soildb
integer(i4), public iflag_soildb
Definition mo_mpr_global_variables.f90:35
mo_mpr_global_variables::horizondepth_mhm
real(dp), dimension(:), allocatable, public horizondepth_mhm
Definition mo_mpr_global_variables.f90:37
mo_mpr_neutrons
Multiscale parameter regionalization (MPR) for neutrons.
Definition mo_mpr_neutrons.f90:12
mo_mpr_neutrons::mpr_neutrons
subroutine, public mpr_neutrons(process_case, param, is_present, nhorizons, ntillhorizons, lcover0, clay, dbm, db, cosmic_l3_till, latwat_till, cosmic_l3, latwat)
multiscale parameter regionalization for neutrons
Definition mo_mpr_neutrons.f90:102
mo_mpr_pet
MPR routine for PET.
Definition mo_mpr_pet.f90:12
mo_mpr_pet::pet_correctbyasp
subroutine, public pet_correctbyasp(id0, latitude_l0, asp0, param, nodata, fasp0)
correction of PET
Definition mo_mpr_pet.f90:215
mo_mpr_pet::bulksurface_resistance
subroutine, public bulksurface_resistance(lai0, param, mask0, nodata, cell_id0, nl0_in_l1, upp_row_l1, low_row_l1, lef_col_l1, rig_col_l1, bulksurface_resistance1)
Regionalization of bulk surface resistance.
Definition mo_mpr_pet.f90:409
mo_mpr_pet::priestley_taylor_alpha
subroutine, public priestley_taylor_alpha(lai0, param, mask0, nodata, cell_id0, nl0_in_l1, upp_row_l1, low_row_l1, lef_col_l1, rig_col_l1, priestley_taylor_alpha1)
Regionalization of priestley taylor alpha.
Definition mo_mpr_pet.f90:313
mo_mpr_pet::pet_correctbylai
subroutine, public pet_correctbylai(param, nodata, lcover0, lai0, mask0, cell_id0, upp_row_l1, low_row_l1, lef_col_l1, rig_col_l1, nl0_in_l1, l1_petlaicorfactor)
estimate PET correction factor based on LAI at L1
Definition mo_mpr_pet.f90:82
mo_mpr_runoff
multiscale parameter regionalization for runoff generation
Definition mo_mpr_runoff.f90:12
mo_mpr_runoff::mpr_runoff
subroutine, public mpr_runoff(lcover0, mask0, sms_fc0, slope_emp0, ksvar_h0, param, cell_id0, upp_row_l1, low_row_l1, lef_col_l1, rig_col_l1, nl0_in_l1, l1_hl1, l1_k0, l1_k1, l1_alpha)
multiscale parameter regionalization for runoff parameters
Definition mo_mpr_runoff.f90:76
mo_mpr_smhorizons
setting up the soil moisture horizons
Definition mo_mpr_smhorizons.f90:12
mo_mpr_smhorizons::mpr_smhorizons
subroutine, public mpr_smhorizons(param, processmatrix, iflag_soil, nhorizons_mhm, horizondepth, lcover0, soilid0, nhorizons, ntillhorizons, thetas_till, thetafc_till, thetapw_till, thetas, thetafc, thetapw, wd, db, dbm, rzdepth, mask0, cell_id0, upp_row_l1, low_row_l1, lef_col_l1, rig_col_l1, nl0_in_l1, l1_beta, l1_sms, l1_fc, l1_pw, l1_froots, latwat_till, cosmic_l3_till, latwat, cosmic_l3, l1_bulkdens, l1_latticewater, l1_cosmicl3)
upscale soil moisture horizons
Definition mo_mpr_smhorizons.f90:132
mo_mpr_soilmoist
Multiscale parameter regionalization (MPR) for soil moisture.
Definition mo_mpr_soilmoist.f90:12
mo_mpr_soilmoist::mpr_sm
subroutine, public mpr_sm(param, processmatrix, is_present, nhorizons, ntillhorizons, sand, clay, dbm, id0, soilid0, lcover0, thetas_till, thetafc_till, thetapw_till, thetas, thetafc, thetapw, ks, db, ksvar_h0, ksvar_v0, sms_fc0)
multiscale parameter regionalization for soil moisture
Definition mo_mpr_soilmoist.f90:102
mo_multi_param_reg
Multiscale parameter regionalization (MPR).
Definition mo_multi_param_reg.f90:12
mo_multi_param_reg::mpr
subroutine, public mpr(mask0, geounit0, soilid0, asp0, gridded_lai0, lcover0, slope_emp0, y0, id0, upper_bound1, lower_bound1, left_bound1, right_bound1, n_subcells1, fsealed1, alpha1, degdayinc1, degdaymax1, degdaynopre1, fasp1, harsamcoeff1, prietayalpha1, aeroresist1, surfresist1, froots1, kfastflow1, kslowflow1, kbaseflow1, kperco1, karstloss1, soilmoistfc1, soilmoistsat1, soilmoistexp1, jarvis_thresh_c1, tempthresh1, unsatthresh1, sealedthresh1, wiltingpoint1, maxinter1, petlaicorfactor, no_count1, bulkdens1, latticewater1, cosmicl31, parameterset)
Regionalizing and Upscaling process parameters.
Definition mo_multi_param_reg.f90:74
mo_multi_param_reg::snow_acc_melt_param
subroutine snow_acc_melt_param(param, fforest1, fiperm1, fperm1, tempthresh1, degdaynopre1, degdayinc1, degdaymax1)
Calculates the snow parameters.
Definition mo_multi_param_reg.f90:802
mo_multi_param_reg::aerodynamical_resistance
subroutine aerodynamical_resistance(lai0, lcover0, param, mask0, id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, aerodyn_resistance1)
Regionalization of aerodynamic resistance.
Definition mo_multi_param_reg.f90:1205
mo_multi_param_reg::karstic_layer
subroutine karstic_layer(param, geounit0, mask0, sms_fc0, ksvar_v0, id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, karstloss1, l1_kp)
calculates the Karstic percolation loss
Definition mo_multi_param_reg.f90:946
mo_multi_param_reg::baseflow_param
subroutine baseflow_param(param, geounit0, k2_0)
baseflow recession parameter
Definition mo_multi_param_reg.f90:690
mo_multi_param_reg::canopy_intercept_param
subroutine, public canopy_intercept_param(processmatrix, param, lai0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, id0, mask0, nodata, max_intercept1)
estimate effective maximum interception capacity at L1
Definition mo_multi_param_reg.f90:1078
mo_multi_param_reg::iper_thres_runoff
subroutine iper_thres_runoff(param, sealedthresh1)
sets the impervious layer threshold parameter for runoff generation
Definition mo_multi_param_reg.f90:884
mo_upscaling_operators
Module containing upscaling operators.
Definition mo_upscaling_operators.f90:12
mo_upscaling_operators::upscale_arithmetic_mean
real(dp) function, dimension(size(nl0_cells_in_l1_cell, 1)), public upscale_arithmetic_mean(nl0_cells_in_l1_cell, l1_upper_rowid_cell, l1_lower_rowid_cell, l1_left_colonid_cell, l1_right_colonid_cell, l0_cellid, mask0, nodata_value, l0_finescale_data)
aritmetic mean
Definition mo_upscaling_operators.f90:268
mo_upscaling_operators::l0_fractionalcover_in_lx
real(dp) function, dimension(size(l0upbound_inlx, 1)), public l0_fractionalcover_in_lx(datain0, classid, mask0, l0upbound_inlx, l0downbound_inlx, l0leftbound_inlx, l0rightbound_inlx, ntcells0_inlx)
fractional coverage of a given class of L0 fields in Lx field (Lx = L1 or L11)
Definition mo_upscaling_operators.f90:154