Line data Source code
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
12 : MODULE 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 :
25 : contains
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
67 621 : subroutine mpr(mask0, geoUnit0, soilId0, Asp0, gridded_LAI0, LCover0, slope_emp0, y0, Id0, upper_bound1, lower_bound1, &
68 1035 : left_bound1, right_bound1, n_subcells1, fSealed1, alpha1, degDayInc1, degDayMax1, degDayNoPre1, fAsp1, &
69 621 : HarSamCoeff1, PrieTayAlpha1, aeroResist1, surfResist1, fRoots1, kFastFlow1, kSlowFlow1, kBaseFlow1, &
70 621 : kPerco1, karstLoss1, soilMoistFC1, soilMoistSat1, soilMoistExp1, jarvis_thresh_c1, tempThresh1, &
71 621 : unsatThresh1, sealedThresh1, wiltingPoint1, maxInter1, petLAIcorFactor, &
72 414 : No_Count1, bulkDens1, latticeWater1, COSMICL31, &
73 207 : 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 207 : real(dp), dimension(:), pointer :: param
181 :
182 207 : real(dp), dimension(:, :, :), allocatable :: thetaS_till
183 :
184 207 : real(dp), dimension(:, :, :), allocatable :: thetaFC_till
185 :
186 207 : real(dp), dimension(:, :, :), allocatable :: thetaPW_till
187 :
188 : ! saturated hydraulic conductivity
189 207 : real(dp), dimension(:, :, :), allocatable :: Ks
190 :
191 : ! Bulk density
192 207 : real(dp), dimension(:, :, :), allocatable :: Db
193 207 : real(dp), dimension(:, :), allocatable :: thetaS
194 207 : real(dp), dimension(:, :), allocatable :: thetaFC
195 207 : real(dp), dimension(:, :), allocatable :: thetaPW
196 :
197 : ! neutron count
198 207 : real(dp), dimension(:,:,:), allocatable :: latWat_till
199 207 : real(dp), dimension(:,:,:), allocatable :: COSMIC_L3_till
200 207 : real(dp), dimension(:,:), allocatable :: latWat ! lattice water
201 207 : real(dp), dimension(:,:), allocatable :: COSMIC_L3 ! COSMIC parameter L3
202 :
203 :
204 : ! relative variability of saturated
205 : ! hydraulic cound. for Horizantal flow
206 207 : real(dp), dimension(:), allocatable :: KsVar_H0
207 :
208 : ! relative variability of saturated
209 : ! hydraulic cound. for vertical flow
210 207 : real(dp), dimension(:), allocatable :: KsVar_V0
211 :
212 : ! soil mositure deficit from
213 : ! field cap. w.r.t to saturation
214 207 : real(dp), dimension(:), allocatable :: SMs_FC0
215 :
216 : ! L0 baseflow parameter
217 9583982 : real(dp), dimension(size(Id0, 1)) :: k2_0
218 :
219 : ! L1 baseflow parameter
220 207 : real(dp), dimension(:), allocatable :: k2_1
221 :
222 : ! L0 Aspect
223 9584189 : 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 8022 : real(dp), dimension(size(fSealed1, dim = 1)) :: fForest1
254 :
255 : ! [1] Fraction of permeable cover
256 8022 : real(dp), dimension(size(fSealed1, dim = 1)) :: fPerm1
257 :
258 :
259 207 : if (present(parameterset)) then
260 207 : param => parameterset
261 : else
262 0 : param => global_parameters(:, 3)
263 : end if
264 :
265 :
266 : ! loop over all LCover scenes
267 621 : 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 414 : n_subcells1)
277 828 : fSealed1(:, 1, iiLC) = L0_fractionalCover_in_Lx(LCover0(:, iiLC), 2, mask0, &
278 : upper_bound1, &
279 : lower_bound1, &
280 : left_bound1, &
281 : right_bound1, &
282 16458 : 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 414 : 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 15630 : 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 16044 : fPerm1(:) = 1.0_dp - fSealed1(:, 1, iiLC) - fForest1(:)
298 :
299 : ! ------------------------------------------------------------------
300 : ! snow parameters
301 : ! ------------------------------------------------------------------
302 828 : select case(processMatrix(2,1))
303 : case(1)
304 :
305 414 : iStart = processMatrix(2, 3) - processMatrix(2, 2) + 1
306 414 : iEnd = processMatrix(2, 3)
307 :
308 0 : 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 414 : )
313 : case DEFAULT
314 414 : 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 414 : msoil = size(soilDB%is_present, 1)
321 19167964 : mLC = maxval(LCover0(:, iiLC), (LCover0(:, iiLC) .ne. nodata_i4))
322 :
323 : ! depending on which kind of soil database processing is to be performed
324 414 : if(iFlag_soilDB .eq. 0)then
325 611064 : mtill = maxval(soilDB%nTillHorizons, (soilDB%nTillHorizons .ne. nodata_i4))
326 611064 : mHor = maxval(soilDB%nHorizons, (soilDB%nHorizons .ne. nodata_i4))
327 0 : 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 0 : mtill = 1
337 0 : mHor = 1
338 : end if
339 :
340 2070 : allocate(thetaS_till(msoil, mtill, mLC))
341 1656 : allocate(thetaFC_till(msoil, mtill, mLC))
342 1656 : allocate(thetaPW_till(msoil, mtill, mLC))
343 1656 : allocate(thetaS(msoil, mHor))
344 1242 : allocate(thetaFC(msoil, mHor))
345 1242 : allocate(thetaPW(msoil, mHor))
346 2070 : allocate(Ks(msoil, mHor, mLC))
347 1656 : allocate(Db(msoil, mHor, mLC))
348 :
349 : ! neutron count related ones
350 : ! allocate and initalize here
351 1656 : allocate( latWat_till(msoil, mtill, mLC))
352 1656 : allocate(COSMIC_L3_till(msoil, mtill, mLC))
353 1242 : allocate( latWat(msoil, mHor ))
354 1242 : allocate( COSMIC_L3(msoil, mHor ))
355 1835262 : latWat_till = 0.000001_dp
356 1835262 : COSMIC_L3_till = 0.000001_dp
357 1834020 : COSMIC_L3 = 0.000001_dp
358 1834020 : 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 1656 : allocate(KsVar_H0(size(Id0, 1)))
366 828 : allocate(KsVar_V0(size(Id0, 1)))
367 828 : allocate( SMs_FC0(size(Id0, 1)))
368 :
369 822 : select case(processMatrix(3,1))
370 : case(1)
371 : ! first thirteen parameters go to this routine
372 408 : iStart = processMatrix(3, 3) - processMatrix(3, 2) + 1
373 408 : 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 408 : iStart2 = processMatrix(3, 3) - 4 + 1
378 408 : iEnd2 = processMatrix(3, 3)
379 :
380 : case(2)
381 : ! first thirteen parameters go to this routine
382 2 : iStart = processMatrix(3, 3) - processMatrix(3, 2) + 1
383 2 : 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 2 : iStart2 = processMatrix(3, 3) - 5 + 1
388 2 : iEnd2 = processMatrix(3, 3) - 1
389 : ! last parameter is jarvis parameter - no need to be regionalized
390 74 : 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 2 : iStart = processMatrix(3, 3) - processMatrix(3, 2) + 1
397 2 : 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 2 : iStart2 = processMatrix(3, 3) - 8
402 2 : iEnd2 = processMatrix(3, 3) - 1
403 :
404 : ! last parameter is jarvis parameter - no need to be regionalized
405 74 : 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 2 : iStart = processMatrix(3, 3) - processMatrix(3, 2) + 1
413 2 : 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 2 : iStart2 = processMatrix(3, 3) - 7
418 2 : 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 414 : 'does not exist! mo_multi_param_reg')
424 : end select
425 :
426 0 : 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 414 : thetaFC, thetaPW, Ks, Db, KsVar_H0, KsVar_V0, SMs_FC0)
432 :
433 : ! >> neutron count related parameters
434 414 : if ( processMatrix(10,1) .GT. 0 ) &
435 : call mpr_neutrons( processMatrix(10,1), & ! IN: processmatrix case
436 0 : 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 12 : )
449 :
450 0 : 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 414 : )
465 :
466 414 : deallocate(thetaS_till)
467 414 : deallocate(thetaFC_till)
468 414 : deallocate(thetaPW_till)
469 414 : deallocate(thetaS)
470 414 : deallocate(thetaFC)
471 414 : deallocate(thetaPW)
472 414 : deallocate(Ks)
473 414 : deallocate(Db)
474 :
475 : ! neutron count
476 414 : deallocate( latWat_till )
477 414 : deallocate( COSMIC_L3_till )
478 414 : deallocate( latWat )
479 414 : deallocate( COSMIC_L3 )
480 :
481 : ! ------------------------------------------------------------------
482 : ! potential evapotranspiration (PET)
483 : ! ------------------------------------------------------------------
484 : ! Penman-Monteith method is only method that is LCscene dependent
485 414 : if (processMatrix(5, 1) == 3) then
486 12 : iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
487 12 : iEnd = processMatrix(5, 3)
488 0 : 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 12 : aeroResist1(:, :, iiLC))
491 402 : else if (processMatrix(5, 1) == -1) then
492 4 : iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
493 4 : iEnd = processMatrix(5, 3)
494 :
495 0 : call pet_correctbyLAI(param(iStart : iEnd), nodata_dp, &
496 : LCover0(:, iiLC), gridded_LAI0, mask0, Id0, &
497 : upper_bound1, lower_bound1, left_bound1, &
498 4 : right_bound1, n_subcells1, petLAIcorFactor(:, :, iiLC))
499 : end if
500 :
501 : ! ------------------------------------------------------------------
502 : ! interflow
503 : ! ------------------------------------------------------------------
504 828 : select case(processMatrix(6, 1))
505 : case (1)
506 : !
507 414 : iStart = processMatrix(6, 3) - processMatrix(6, 2) + 1
508 414 : 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 0 : KsVar_H0, param(iStart : iEnd), Id0, upper_bound1, lower_bound1, &
513 : left_bound1, right_bound1, n_subcells1, unsatThresh1(:, 1, 1), kFastFlow1(:, 1, iiLC), &
514 414 : kSlowFlow1(:, 1, iiLC), alpha1(:, 1, iiLC))
515 : case DEFAULT
516 414 : 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 828 : select case(processMatrix(7, 1))
523 : case(1)
524 :
525 414 : iStart = processMatrix(7, 3) - processMatrix(7, 2) + 1
526 414 : iEnd = processMatrix(7, 3)
527 : call karstic_layer(& ! In
528 0 : 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 414 : )
534 :
535 : case DEFAULT
536 414 : call error_message('***ERROR: Process description for process "percolation" does not exist! mo_multi_param_reg')
537 : end select
538 :
539 414 : deallocate(KsVar_H0)
540 414 : deallocate(KsVar_V0)
541 621 : deallocate(SMs_FC0)
542 :
543 : end do !! >>>>>>> LAND COVER SCENE LOOP
544 :
545 :
546 : ! ------------------------------------------------------------------
547 : ! sealed area threshold for runoff generation
548 : ! ------------------------------------------------------------------
549 414 : select case(processMatrix(4, 1))
550 : case (1)
551 207 : iStart = processMatrix(4, 3) - processMatrix(4, 2) + 1
552 207 : iEnd = processMatrix(4, 3)
553 207 : call iper_thres_runoff(param(iStart : iEnd), sealedThresh1)
554 : case DEFAULT
555 207 : 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 207 : select case(processMatrix(5, 1))
562 : case(-1) ! LAI correction of input PET
563 404 : iEnd = -9999 ! dummy statement
564 : case(0) ! aspect correction of input PET
565 197 : iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
566 197 : iEnd = processMatrix(5, 3)
567 197 : call pet_correctbyASP(Id0, y0, Asp0, param(iStart : iEnd), nodata_dp, fAsp0)
568 197 : fAsp1(:, 1, 1) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
569 7662 : left_bound1, right_bound1, Id0, mask0, nodata_dp, fAsp0)
570 : case(1) ! Hargreaves-Samani method
571 1 : iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
572 1 : iEnd = processMatrix(5, 3)
573 1 : call pet_correctbyASP(Id0, y0, Asp0, param(iStart : iEnd - 1), nodata_dp, fAsp0)
574 1 : fAsp1(:, 1, 1) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
575 36 : left_bound1, right_bound1, Id0, mask0, nodata_dp, fAsp0)
576 37 : HarSamCoeff1 = param(iEnd)
577 : case(2) ! Priestley-Taylor Method
578 1 : iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
579 1 : iEnd = processMatrix(5, 3)
580 0 : 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 1 : PrieTayAlpha1(:, :, 1))
583 : case(3) ! Penman-Monteith method
584 : ! aerodynamic resistance is calculated inside LCscene loop
585 6 : iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
586 6 : iEnd = processMatrix(5, 3)
587 0 : call bulksurface_resistance(gridded_LAI0, param(iEnd), mask0, &
588 : nodata_dp, Id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
589 6 : surfResist1(:, :, 1))
590 : case default
591 207 : 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 414 : 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 207 : iStart = processMatrix(9, 3) - processMatrix(9, 2) + 1
601 207 : iEnd = processMatrix(9, 3)
602 :
603 207 : call baseflow_param(param(iStart : iEnd), geoUnit0, k2_0)
604 :
605 : ! Upscale by arithmetic mean
606 621 : allocate(k2_1(size(kBaseFlow1, 1)))
607 : k2_1 = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
608 207 : left_bound1, right_bound1, Id0, mask0, nodata_dp, k2_0)
609 : ! loop over all LCover scenes
610 621 : do iiLC = 1, size(LCover0, 2)
611 15837 : kBaseFlow1(:, 1, iiLC) = k2_1
612 : end do
613 207 : 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 16458 : if (processMatrix(7, 1) .gt. 0) kBaseFlow1 = merge(kSlowFlow1, kBaseFlow1, kBaseFlow1 .lt. kSlowFlow1)
619 :
620 : case DEFAULT
621 207 : 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 213 : 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 6 : iStart = processMatrix(10, 3) - processMatrix(10, 2) + 1
634 6 : iEnd = processMatrix(10, 3)
635 672 : 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 0 : iStart = processMatrix(10, 3) - processMatrix(10, 2) + 1
639 0 : iEnd = processMatrix(10, 3)
640 0 : No_Count1 = param(iStart) ! >> 1st parameter --> N0 parameter
641 : case DEFAULT
642 207 : 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 207 : nodata_dp, maxInter1(:, :, 1))
653 :
654 414 : 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 414 : subroutine baseflow_param(param, geoUnit0, k2_0)
690 :
691 207 : 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 207 : if (size(param) .ne. size(geoUnitList)) &
713 0 : call error_message(' mo_multi_param_reg: baseflow_param: size mismatch, subroutine baseflow parameters ')
714 :
715 9583982 : k2_0 = nodata_dp
716 :
717 : !$OMP PARALLEL
718 : !$OMP DO PRIVATE(gg) SCHEDULE(STATIC)
719 9583982 : do ii = 1, size(k2_0)
720 : ! get parameter index in geoUnitList
721 105421525 : gg = minloc(abs(geoUnitList - geoUnit0(ii)))
722 9583982 : k2_0(ii) = param(gg(1))
723 : end do
724 : !$OMP END DO
725 : !$OMP END PARALLEL
726 :
727 207 : 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 828 : subroutine snow_acc_melt_param(param, fForest1, fIperm1, fPerm1, tempThresh1, degDayNoPre1, degDayInc1, &
801 414 : 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 414 : real(dp) :: tmp_degreeDayFactor_forest, tmp_degreeDayFactor_impervious, tmp_degreeDayFactor_pervious
829 :
830 414 : real(dp) :: tmp_maxDegreeDayFactor_forest, tmp_maxDegreeDayFactor_impervious, tmp_maxDegreeDayFactor_pervious
831 :
832 :
833 414 : tmp_degreeDayFactor_forest = param(2) ! OLD: param(2)
834 414 : tmp_degreeDayFactor_impervious = param(2) + param(4) + param(3) ! OLD: param(3)
835 414 : tmp_degreeDayFactor_pervious = param(2) + param(4) ! OLD: param(4)
836 414 : tmp_maxDegreeDayFactor_forest = param(2) + param(6) ! OLD: param(6)
837 414 : tmp_maxDegreeDayFactor_impervious = param(2) + param(4) + param(3) + param(7) ! OLD: param(7)
838 414 : tmp_maxDegreeDayFactor_pervious = param(2) + param(4) + param(8) ! OLD: param(8)
839 :
840 15630 : tempThresh1 = param(1)
841 15630 : degDayInc1 = param(5)
842 :
843 0 : degDayNoPre1 = (&
844 : tmp_degreeDayFactor_forest * fForest1 + &
845 0 : tmp_degreeDayFactor_impervious * fIperm1 + &
846 16044 : tmp_degreeDayFactor_pervious * fPerm1)
847 0 : degDayMax1 = (&
848 : tmp_maxDegreeDayFactor_forest * fForest1 + &
849 : tmp_maxDegreeDayFactor_impervious * fIperm1 + &
850 16044 : tmp_maxDegreeDayFactor_pervious * fPerm1)
851 :
852 207 : 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 207 : 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 8229 : sealedThresh1 = param(1)
893 :
894 414 : 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 1242 : subroutine karstic_layer(param, geoUnit0, mask0, SMs_FC0, KsVar_V0, Id0, n_subcells1, upper_bound1, lower_bound1, &
945 828 : left_bound1, right_bound1, karstLoss1, L1_Kp)
946 :
947 207 : 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 414 : real(dp), dimension(:), allocatable :: fKarArea
995 :
996 : ! temporal variable
997 19168378 : 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 414 : tmp = merge(param(1) * (1.0_dp + SMs_FC0) / (1.0_dp + KsVar_V0), &
1010 19168378 : nodata_dp, Id0 .ne. nodata_i4)
1011 : !$OMP END PARALLEL
1012 :
1013 : L1_Kp = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
1014 414 : left_bound1, right_bound1, Id0, mask0, nodata_dp, tmp)
1015 :
1016 : ! minimum constrains
1017 15630 : L1_Kp = merge(2.0_dp, L1_Kp, L1_Kp .lt. 2.0_dp)
1018 :
1019 414 : nGeoUnits = size(geoUnitlist, 1)
1020 :
1021 : ! 1st calculate fraction of Karstic area
1022 1242 : allocate(fKarArea(size(karstLoss1, 1)))
1023 16044 : fKarArea = 0.0_dp
1024 :
1025 4554 : do i = 1, nGeoUnits
1026 4140 : if(GeoUnitKar(i) .eq. 0) cycle
1027 0 : fKarArea(:) = L0_fractionalCover_in_Lx(geoUnit0, geoUnitlist(i), mask0, &
1028 4554 : upper_bound1, lower_bound1, left_bound1, right_bound1, n_subcells1)
1029 : end do
1030 :
1031 : ! 2nd calculate karstLoss1
1032 15630 : karstLoss1 = 1.0_dp - (fKarArea * param(2))
1033 :
1034 414 : deallocate(fKarArea)
1035 :
1036 414 : 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 621 : subroutine canopy_intercept_param(processMatrix, param, LAI0, n_subcells1, upper_bound1, lower_bound1, left_bound1, &
1077 621 : right_bound1, Id0, mask0, nodata, max_intercept1)
1078 :
1079 414 : 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 207 : real(dp), dimension(:), allocatable :: max_intercept0
1123 :
1124 207 : real(dp), dimension(:), allocatable :: gamma_intercept
1125 :
1126 :
1127 : ! ------------------------------------------------------------------
1128 : ! Maximum interception parameter
1129 : ! ------------------------------------------------------------------
1130 207 : select case(processMatrix(1, 1))
1131 : case(1)
1132 207 : iStart = processMatrix(1, 3) - processMatrix(1, 2) + 1
1133 207 : iEnd = processMatrix(1, 3)
1134 :
1135 : ! allocate space
1136 621 : allocate(gamma_intercept(iEnd - iStart + 1))
1137 621 : allocate(max_intercept0 (size(Id0, 1)))
1138 :
1139 : ! estimate max. intercept at Level-0
1140 414 : gamma_intercept(:) = param(iStart : iEnd)
1141 :
1142 2709 : do it = 1, size(LAI0, 2)
1143 : !$OMP PARALLEL
1144 115845612 : 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 2709 : right_bound1, Id0, mask0, nodata, max_intercept0(:))
1150 :
1151 : end do
1152 :
1153 207 : deallocate(gamma_intercept)
1154 207 : deallocate(max_intercept0)
1155 : CASE DEFAULT
1156 207 : call error_message('mo_multi_param_reg: This processMatrix=', num2str(processMatrix(1, 1)), ' is not implemented!')
1157 : end select
1158 :
1159 207 : 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 36 : subroutine aerodynamical_resistance(LAI0, LCover0, param, mask0, Id0, n_subcells1, upper_bound1, lower_bound1, &
1204 24 : left_bound1, right_bound1, aerodyn_resistance1)
1205 :
1206 207 : 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 12 : real(dp), dimension(:), allocatable :: maxLAI
1248 :
1249 12 : real(dp), dimension(:), allocatable :: zm
1250 :
1251 12 : real(dp), dimension(:), allocatable :: canopy_height0
1252 :
1253 12 : 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 12 : real(dp), dimension(:, :), allocatable :: aerodyn_resistance0
1258 :
1259 :
1260 : ! initialize some things
1261 558576 : allocate(zm (size(LCover0, dim = 1))) ; zm = nodata_dp
1262 558564 : allocate(zm_zero (size(LCover0, dim = 1))) ; zm_zero = nodata_dp
1263 558564 : allocate(zh_zero (size(LCover0, dim = 1))) ; zh_zero = nodata_dp
1264 558564 : allocate(displace (size(LCover0, dim = 1))) ; displace = nodata_dp
1265 558564 : allocate(canopy_height0 (size(LCover0, dim = 1))) ; canopy_height0 = nodata_dp
1266 6702672 : allocate(aerodyn_resistance0 (size(LCover0, dim = 1), size(LAI0, 2))) ; aerodyn_resistance0 = nodata_dp
1267 558564 : allocate(maxLAI (size(LCover0, dim = 1))) ; maxLAI = nodata_dp
1268 :
1269 : ! regionalization of canopy height
1270 : ! substitute with canopy height
1271 558564 : canopy_height0 = merge(param(1), canopy_height0, LCover0 == 1) ! forest
1272 558564 : 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 12 : maxLAI = MAXVAL(LAI0, dim=2)
1277 :
1278 156 : do tt = 1, size(LAI0, 2)
1279 :
1280 : ! pervious canopy height is scaled with LAI
1281 6702768 : 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 6702624 : zm = WindMeasHeight
1286 : ! correction: if wind measurement height is below canopy height loagarithm becomes negative
1287 6702768 : zm = merge(canopy_height0 + zm, zm, ((abs(zm - nodata_dp) .GT. eps_dp) .AND. (zm .LT. canopy_height0)))
1288 :
1289 : ! zh = zm
1290 6702768 : displace = param(4) * canopy_height0
1291 6702768 : zm_zero = param(5) * canopy_height0
1292 6702768 : zh_zero = param(6) * zm_zero
1293 : !
1294 : ! calculate aerodynamic resistance (changes monthly)
1295 6702624 : 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 156 : left_bound1, right_bound1, Id0, mask0, nodata_dp, aerodyn_resistance0(:, tt))
1298 :
1299 : end do
1300 :
1301 12 : end subroutine aerodynamical_resistance
1302 :
1303 : END MODULE mo_multi_param_reg
|
