Line data Source code
1 : !> \file mo_mpr_startup.f90
2 : !> \brief \copybrief mo_mpr_startup
3 : !> \details \copydetails mo_mpr_startup
4 :
5 : !> \brief Startup procedures for mHM.
6 : !> \details This module initializes all variables required to run mHM. This
7 : !! module needs to be run only one time at the beginning of a simulation if re-starting files do not exist.
8 : !> \authors Luis Samaniego, Rohini Kumar
9 : !> \date Dec 2012
10 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
11 : !! mHM is released under the LGPLv3+ license \license_note
12 : !> \ingroup f_mpr
13 : MODULE mo_mpr_startup
14 :
15 : USE mo_kind, ONLY : i4, dp
16 : use mo_common_constants, only : nodata_i4, nodata_dp ! global nodata values (i4, dp)
17 :
18 : IMPLICIT NONE
19 :
20 : PRIVATE
21 :
22 : PUBLIC :: mpr_initialize, init_eff_params ! initialization sequence
23 :
24 : CONTAINS
25 :
26 : !> \brief Initialize main mHM variables
27 : !> \details Initialize main mHM variables for a given domain.
28 : !! Calls the following procedures in this order:
29 : !! - Constant initialization.
30 : !! - Generate soil database.
31 : !! - Checking inconsistencies input fields.
32 : !! - Variable initialization at level-0.
33 : !! - Variable initialization at level-1.
34 : !! - Variable initialization at level-11.
35 : !! - Space allocation of remaining variable/parameters.
36 : !! Global variables will be used at this stage.
37 : !> \changelog
38 : !! - Luis Samaniego Mar 2008
39 : !! - fully distributed multilayer
40 : !! - Rohini Kumar Oct 2010
41 : !! - matrix to vector version
42 : !! - openmp parallelization
43 : !! - routing level 11
44 : !! - Luis Samaniego Jul 2012
45 : !! - removal of IMSL dependencies
46 : !! - Luis Samaniego Dec 2012
47 : !! - modular version
48 : !! - Rohini Kumar May 2013
49 : !! - code cleaned and error checks
50 : !! - Rohini Kumar Nov 2013
51 : !! - updated documentation
52 : !! - Stephan Thober Jun 2014
53 : !! - copied L2 initialization from mo_meteo_forcings
54 : !! - Stephan Thober Jun 2014
55 : !! - updated flag for read_restart
56 : !! - Stephan Thober Aug 2015
57 : !! - removed initialisation of routing
58 : !! - Rohini Kumar Mar 2016
59 : !! - changes for handling multiple soil database options
60 : !! - Robert Schweppe Jun 2018
61 : !! - refactoring and reformatting
62 : !> \authors Luis Samaniego, Rohini Kumar
63 : !> \date Dec 2012
64 13 : subroutine mpr_initialize
65 :
66 : use mo_common_variables, only : l0_l1_remap, level0, level1, domainMeta, resolutionHydrology
67 : use mo_grid, only : init_lowres_level, set_domain_indices
68 : use mo_kind, only : i4
69 : use mo_read_latlon, only : read_latlon
70 : use mo_soil_database, only : generate_soil_database
71 :
72 : implicit none
73 :
74 : integer(i4) :: iDomain
75 :
76 :
77 : ! soilDB common for all domains
78 13 : call generate_soil_database()
79 :
80 64 : allocate(level1(domainMeta%nDomains))
81 64 : allocate(l0_l1_remap(domainMeta%nDomains))
82 :
83 : ! L0 and L1 initialization
84 38 : do iDomain = 1, domainMeta%nDomains
85 25 : if (iDomain .eq. 1) then
86 13 : call L0_check_input(domainMeta%L0DataFrom(iDomain))
87 13 : call L0_variable_init(domainMeta%L0DataFrom(iDomain))
88 : ! ToDo: adopt to parallel
89 : ! ToDo: check change
90 12 : else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
91 : ! this needs only be done if there is new input
92 : call L0_check_input(domainMeta%L0DataFrom(iDomain))
93 9 : call L0_variable_init(domainMeta%L0DataFrom(iDomain))
94 : end if
95 :
96 0 : call init_lowres_level(level0(domainMeta%L0DataFrom(iDomain)), resolutionHydrology(iDomain), &
97 25 : level1(iDomain), l0_l1_remap(iDomain))
98 : ! read lat lon coordinates for level 1
99 25 : call read_latlon(iDomain, "lon", "lat", "level1", level1(iDomain))
100 :
101 : ! Parameter fields have to be allocated in any case
102 63 : call init_eff_params(level1(iDomain)%nCells)
103 :
104 : end do
105 :
106 : call set_domain_indices(level1)
107 :
108 13 : end subroutine mpr_initialize
109 :
110 :
111 : !> \brief Check for errors in L0 input data
112 : !> \details Check for possible errors in input data (morphological and land cover) at level-0
113 : !> \changelog
114 : !! - Rohini Kumar Aug 2013
115 : !! - added iFlag_LAI_data_format to handle LAI options, and changed within the code made accordingly
116 : !! - Rohini Kumar Sep 2013
117 : !! - read input data for routing processes according & Stephan Thober, to process_matrix flag
118 : !! - Stephan Thober Aug 2015
119 : !! - moved check of L0 routing variables to mRM
120 : !! - Rohini Kumar Mar 2016
121 : !! - changes for handling multiple soil database options
122 : !! - Robert Schweppe Jun 2018
123 : !! - refactoring and reformatting
124 : !> \authors Rohini Kumar
125 : !> \date Jan 2013
126 22 : subroutine L0_check_input(iDomain)
127 :
128 13 : use mo_common_constants, only : eps_dp
129 : use mo_common_variables, only : L0_LCover, L0_elev, level0, nLCoverScene
130 : use mo_message, only : error_message
131 : use mo_mpr_global_variables, only : L0_asp, L0_geoUnit, L0_gridded_LAI, &
132 : L0_slope, L0_soilId, iFlag_soilDB, nSoilHorizons_mHM, timeStep_LAI_input
133 : use mo_string_utils, only : num2str
134 : use mo_utils, only : eq
135 :
136 : implicit none
137 :
138 : !> domain id
139 : integer(i4), intent(in) :: iDomain
140 :
141 : integer(i4) :: k, n, nH
142 :
143 : CHARACTER(len=1024) :: message_text = ''
144 :
145 : ! START CHECKING VARIABLES
146 972972 : do k = level0(iDomain)%iStart, level0(iDomain)%iEnd
147 :
148 : ! elevation [m]
149 972950 : if (abs(L0_elev(k) - nodata_dp) .lt. eps_dp) then
150 0 : message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(iDomain, '(I5)'))
151 : call error_message(' Error: elevation has missing value within the valid masked area at cell in domain ', &
152 0 : trim(message_text))
153 : end if
154 :
155 : ! slope [%]
156 972950 : if (abs(L0_slope(k) - nodata_dp) .lt. eps_dp) then
157 0 : message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(iDomain, '(I5)'))
158 : call error_message(' Error: slope has missing value within the valid masked area at cell in domain ', &
159 0 : trim(message_text))
160 : end if
161 :
162 : ! aspect [degree]
163 972950 : if (abs(L0_asp(k) - nodata_dp) .lt. eps_dp) then
164 0 : message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(iDomain, '(I5)'))
165 : call error_message(' Error: aspect has missing values within the valid masked area at cell in domain ', &
166 0 : trim(message_text))
167 : end if
168 :
169 : ! soil-Id [-]
170 972950 : nH = 1 !> by default; when iFlag_soilDB = 0
171 972950 : if (iFlag_soilDB .eq. 1) nH = nSoilHorizons_mHM
172 : ! another option to handle multiple soil horizons properties
173 1945900 : do n = 1, nH
174 1945900 : if (L0_soilId(k, n) .eq. nodata_i4) then
175 0 : message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(iDomain, '(I5)')) // ',' // trim(num2str(n, '(I5)'))
176 : call error_message(' Error: soil id has missing values within the valid masked area at cell in domain and horizon ', &
177 0 : trim(message_text))
178 : end if
179 : end do
180 :
181 : ! geological-Id [-]
182 972950 : if (L0_geoUnit(k) .eq. nodata_i4) then
183 0 : message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(iDomain, '(I5)'))
184 : call error_message(' Error: geological formation id has missing values within the valid masked area at cell in domain ', &
185 0 : trim(message_text))
186 : end if
187 :
188 : ! landcover scenes
189 2918850 : do n = 1, nLCoverScene
190 2918850 : if (L0_LCover(k, n) .eq. nodata_i4) then
191 0 : message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(iDomain, '(I5)')) // ',' // trim(num2str(n, '(I5)'))
192 : call error_message(' Error: land cover id has missing values within the valid masked area at cell in domain and scene ', &
193 0 : trim(message_text))
194 : end if
195 : end do
196 :
197 : ! land cover scenes related to LAI
198 972972 : if(timeStep_LAI_input .EQ. 0) then
199 926405 : if (eq(L0_gridded_LAI(k, 1), nodata_dp)) then
200 0 : message_text = trim(num2str(k, '(G5.3)')) // ',' // trim(num2str(iDomain, '(I5)'))
201 : call error_message(' Error: gridded LAI has missing values within the valid masked area at cell in domain ', &
202 0 : trim(message_text))
203 : end if
204 : end if
205 :
206 : end do
207 :
208 22 : end subroutine L0_check_input
209 :
210 :
211 : !> \brief level 0 variable initialization
212 : !> \details following tasks are performed for L0 data sets
213 : !! - cell id & numbering
214 : !! - storage of cell cordinates (row and coloum id)
215 : !! - empirical dist. of terrain slope
216 : !! - flag to determine the presence of a particular soil id
217 : !! in this configuration of the model run
218 : !! If a variable is added or removed here, then it also has to
219 : !! be added or removed in the subroutine config_variables_set in
220 : !! module mo_restart and in the subroutine set_config in module
221 : !! mo_set_netcdf_restart
222 : !> \changelog
223 : !! - Rohini Kumar & Matthias Cuntz May 2014
224 : !! - cell area calulation based on a regular lat-lon grid or on a regular X-Y coordinate system
225 : !! - Matthias Cuntz May 2014
226 : !! - changed empirical distribution function so that doubles get the same value
227 : !! - Matthias Zink & Matthias Cuntz Feb 2016
228 : !! - code speed up due to reformulation of CDF calculation
229 : !! - Rohini Kumar Mar 2016
230 : !! - changes for handling multiple soil database options
231 : !! - Maren Kaluza Feb 2018
232 : !! - removed slope_val, temp, only sort the index to speed up finding the empirical distribution slope_emp
233 : !! - Robert Schweppe Jun 2018
234 : !! - refactoring and reformatting
235 : !> \authors Rohini Kumar
236 : !> \date Jan 2013
237 22 : subroutine L0_variable_init(iDomain)
238 :
239 22 : use mo_append, only : append
240 : use mo_common_variables, only : level0
241 : use mo_grid, only : L0_grid_setup
242 : use mo_mpr_global_variables, only : L0_slope, L0_slope_emp, L0_soilId, iFlag_soilDB, nSoilHorizons_mHM, &
243 : nSoilTypes, soilDB
244 : use mo_orderpack, only : sort_index
245 : use mo_utils, only : eq
246 :
247 : implicit none
248 :
249 : !> domain id
250 : integer(i4), intent(in) :: iDomain
251 :
252 22 : real(dp), dimension(:), allocatable :: slope_emp
253 :
254 22 : integer(i4), dimension(:), allocatable :: slope_sorted_index
255 :
256 : integer(i4) :: i, j, k, nH, i_sort, i_sortpost
257 :
258 : ! STEPS ::
259 :
260 :
261 : !--------------------------------------------------------
262 : ! 1) Estimate each variable locally for a given domain
263 : ! 2) Pad each variable to its corresponding global one
264 : !--------------------------------------------------------
265 : !------------------------------------------------------
266 : ! Assign whether a given soil type is present or not
267 : !------------------------------------------------------
268 22 : if (iDomain .eq. 1) then
269 39 : allocate(soilDB%is_present(nSoilTypes))
270 19188 : soilDB%is_present(:) = 0_i4
271 : end if
272 :
273 22 : call L0_grid_setup(level0(iDomain))
274 :
275 : !---------------------------------------------------
276 : ! Estimate empirical distribution of slope
277 : !---------------------------------------------------
278 110 : allocate(slope_emp(level0(iDomain)%nCells), slope_sorted_index(level0(iDomain)%nCells))
279 :
280 : ! get sorted data and sorted indexes to remap later
281 22 : slope_sorted_index = sort_index(L0_slope(level0(iDomain)%iStart : level0(iDomain)%iEnd))
282 :
283 : ! empirical distribution of slopes = cumulated number points with slopes that are <= the slope at this point
284 : !
285 : ! sorted data emp. CDF
286 : ! 9 | x x 7/8 | x x
287 : ! | |
288 : ! 8 | x 5/8 | x
289 : ! | |
290 : ! 5 | x x x 4/8 | x x x
291 : ! | |
292 : ! 2 | x 1/8 | x
293 : ! |__________________ |__________________
294 : !
295 : ! highest slope value = highest rank or No. of data points / (data points + 1)
296 0 : slope_emp(slope_sorted_index(level0(iDomain)%nCells)) = real(level0(iDomain)%nCells, dp) / &
297 22 : real(level0(iDomain)%nCells + 1_i4, dp)
298 :
299 : ! backward loop to check if the preceding data point has the same slope value
300 972950 : do i = level0(iDomain)%nCells - 1, 1, -1
301 972928 : i_sort=slope_sorted_index(i)
302 972928 : i_sortpost=slope_sorted_index(i+1)
303 972950 : if (eq(L0_slope(level0(iDomain)%iStart-1_i4+i_sort), L0_slope(level0(iDomain)%iStart-1_i4+i_sortpost))) then
304 : ! if yes: assign the same probabitity
305 971948 : slope_emp(i_sort) = slope_emp(i_sortpost)
306 : else
307 : ! if not: assign rank / (data points + 1)
308 980 : slope_emp(i_sort) = real(i, dp) / real(level0(iDomain)%nCells + 1_i4, dp)
309 : end if
310 : end do
311 :
312 : ! EXAMPLE
313 : ! in = [ 7, 20, 31, 31, 12, 31, 42 ]
314 : ! sorted = [ 7, 12, 20, 31, 31, 31, 42 ]
315 : ! index = [ 1, 5, 2, 3, 4, 6, 7 ]
316 : ! temp = [ 1, 2, 3, 6, 6, 6, 7 ]
317 : ! out = [ 1, 3, 6, 6, 2, 6, 7 ] / (len(out) + 1 )
318 :
319 : !--------------------------------------------------------
320 : ! Start padding up local variables to global variables
321 : !--------------------------------------------------------
322 22 : call append(L0_slope_emp, slope_emp)
323 :
324 22 : nH = 1_i4 !> by default; when iFlag_soilDB = 0
325 22 : if (iFlag_soilDB .eq. 1) nH = nSoilHorizons_mHM
326 44 : do i = 1, nH
327 972994 : do k = level0(iDomain)%iStart, level0(iDomain)%iEnd
328 972950 : j = L0_soilId(k, i)
329 972972 : soilDB%is_present(j) = 1_i4
330 : end do
331 : end do
332 :
333 : ! free space
334 22 : deallocate(slope_emp, slope_sorted_index)
335 :
336 :
337 44 : end subroutine L0_variable_init
338 :
339 :
340 : !> \brief Allocation of space for mHM related L1 and L11 variables.
341 : !> \details Allocation of space for mHM related L1 and L11 variables (e.g., states,
342 : !! fluxes, and parameters) for a given domain. Variables allocated here is
343 : !! defined in them mo_global_variables.f90 file. After allocating any variable
344 : !! in this routine, initalize them in the following variables_default_init
345 : !! subroutine.
346 : !> \changelog
347 : !! - R. Kumar Sep 2013
348 : !! - documentation added according to the template
349 : !! - S. Thober Aug 2015
350 : !! - removed routing related variables
351 : !! - Zink M. Demirel C. Mar 2017
352 : !! - Init Jarvis soil water stress variable at SM process(3)
353 : !! - Robert Schweppe Dec 2017
354 : !! - restructured allocation in variables_alloc, expanded dimensions of effective parameters
355 : !! - Robert Schweppe Jun 2018
356 : !! - refactoring and reformatting
357 : !! - Rohini Kumar Oct 2021
358 : !! - Added Neutron count module to mHM integrate into develop branch (5.11.2)
359 : !! - Sebastian Müller Mar 2023
360 : !! - made L1_alpha, L1_kSlowFlow, L1_kBaseFlow and L1_kPerco land cover dependent
361 : !> \authors Rohini Kumar
362 : !> \date Jan 2013
363 26 : subroutine init_eff_params(ncells1)
364 :
365 22 : use mo_append, only : append
366 : use mo_constants, only : YearMonths
367 : use mo_common_constants, only : P1_InitStateFluxes
368 : use mo_common_variables, only : nLCoverScene
369 : use mo_mpr_global_variables, only : L1_HarSamCoeff, L1_PrieTayAlpha, L1_aeroResist, L1_alpha, L1_degDay, &
370 : L1_degDayInc, L1_degDayMax, L1_degDayNoPre, L1_fAsp, L1_fRoots, L1_fSealed, &
371 : L1_jarvis_thresh_c1, L1_kBaseFlow, L1_kPerco, L1_kSlowFlow, L1_karstLoss, &
372 : L1_kfastFlow, L1_maxInter, L1_petLAIcorFactor, L1_sealedThresh, L1_soilMoistExp, &
373 : L1_soilMoistFC, L1_soilMoistSat, L1_surfResist, L1_tempThresh, L1_unsatThresh, &
374 : L1_wiltingPoint, nLAI, nSoilHorizons_mHM, &
375 : L1_No_Count, L1_bulkDens, L1_latticeWater, L1_COSMICL3
376 :
377 : implicit none
378 :
379 : !> number of L1 cells
380 : integer(i4), intent(in) :: ncells1
381 :
382 26 : real(dp), dimension(:, :, :), allocatable :: dummy_3D
383 :
384 : integer(i4) :: max_extent
385 :
386 :
387 : ! get maximum extent of one dimension 2 or 3
388 26 : max_extent = max(nSoilHorizons_mHM, int(YearMonths, i4), nLCoverScene, nLAI)
389 :
390 : ! for appending and intialization
391 130 : allocate(dummy_3D(nCells1, max_extent, nLCoverScene))
392 :
393 27884 : dummy_3D = P1_InitStateFluxes
394 :
395 : !-------------------------------------------
396 : ! EFFECTIVE PARAMETERS
397 : !-------------------------------------------
398 26 : call append(L1_fSealed, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
399 : ! exponent for the upper reservoir
400 26 : call append(L1_alpha, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
401 : ! increase of the Degree-day factor per mm of increase in precipitation
402 26 : call append(L1_degDayInc, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
403 : ! maximum degree-day factor
404 26 : call append(L1_degDayMax, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
405 : ! degree-day factor with no precipitation
406 26 : call append(L1_degDayNoPre, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
407 : ! degree-day factor
408 26 : call append(L1_degDay, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
409 : ! Karstic percolation loss
410 26 : call append(L1_karstLoss, dummy_3D(:, 1 : 1, 1 : 1))
411 : ! PET correction factor due to terrain aspect
412 26 : call append(L1_fAsp, dummy_3D(:, 1 : 1, 1 : 1))
413 : ! PET correction factor due to LAI
414 26 : call append(L1_petLAIcorFactor, dummy_3D(:, 1 : nLAI, 1 : nLCoverScene))
415 : ! PET Hargreaves Samani coefficient
416 26 : call append(L1_HarSamCoeff, dummy_3D(:, 1 : 1, 1 : 1))
417 : ! PET Prietley Taylor coefficient
418 26 : call append(L1_PrieTayAlpha, dummy_3D(:, 1 : nLAI, 1 : 1))
419 : ! PET aerodynamical resistance
420 26 : call append(L1_aeroResist, dummy_3D(:, 1 : nLAI, 1 : nLCoverScene))
421 : ! PET bulk surface resistance
422 26 : call append(L1_surfResist, dummy_3D(:, 1 : nLAI, 1 : 1))
423 : ! Fraction of roots in soil horizons
424 26 : call append(L1_fRoots, dummy_3D(:, 1 : nSoilHorizons_mHM, 1 : nLCoverScene))
425 : ! Maximum interception
426 26 : call append(L1_maxInter, dummy_3D(:, 1 : nLAI, 1 : 1))
427 : ! fast interflow recession coefficient
428 26 : call append(L1_kfastFlow, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
429 : ! slow interflow recession coefficient
430 26 : call append(L1_kSlowFlow, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
431 : ! baseflow recession coefficient
432 26 : call append(L1_kBaseFlow, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
433 : ! percolation coefficient
434 26 : call append(L1_kPerco, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
435 : ! Soil moisture below which actual ET is reduced linearly till PWP
436 26 : call append(L1_soilMoistFC, dummy_3D(:, 1 : nSoilHorizons_mHM, 1 : nLCoverScene))
437 : ! Saturation soil moisture for each horizon [mm]
438 26 : call append(L1_soilMoistSat, dummy_3D(:, 1 : nSoilHorizons_mHM, 1 : nLCoverScene))
439 : ! jarvis critical value for normalized soil water content
440 26 : call append(L1_jarvis_thresh_c1, dummy_3D(:, 1 : 1, 1 : 1))
441 : ! Exponential parameter to how non-linear is the soil water retention
442 26 : call append(L1_soilMoistExp, dummy_3D(:, 1 : nSoilHorizons_mHM, 1 : nLCoverScene))
443 : ! Threshold temperature for snow/rain
444 26 : call append(L1_tempThresh, dummy_3D(:, 1 : 1, 1 : nLCoverScene))
445 : ! Threshold water depth controlling fast interflow
446 26 : call append(L1_unsatThresh, dummy_3D(:, 1 : 1, 1 : 1))
447 : ! Threshold water depth for surface runoff in sealed surfaces
448 26 : call append(L1_sealedThresh, dummy_3D(:, 1 : 1, 1 : 1))
449 : ! Permanent wilting point
450 26 : call append(L1_wiltingPoint, dummy_3D(:, 1 : nSoilHorizons_mHM, 1 : nLCoverScene))
451 : ! neutron count related parameters
452 26 : call append(L1_No_Count, dummy_3D(:, 1:1, 1:1)) ! N0 count
453 26 : call append(L1_bulkDens, dummy_3D(:, 1:nSoilHorizons_mHM, 1:nLCoverScene)) ! bulk density
454 26 : call append(L1_latticeWater, dummy_3D(:, 1:nSoilHorizons_mHM, 1:nLCoverScene)) ! lattice water
455 26 : call append(L1_COSMICL3, dummy_3D(:, 1:nSoilHorizons_mHM, 1:nLCoverScene)) ! cosmic L3 parameter
456 :
457 : ! free space
458 26 : if (allocated(dummy_3D)) deallocate(dummy_3D)
459 :
460 26 : end subroutine init_eff_params
461 :
462 : END MODULE mo_mpr_startup
|