Line data Source code
1 : !> \file mo_soil_database.f90
2 : !> \brief \copybrief mo_soil_database
3 : !> \details \copydetails mo_soil_database
4 :
5 : !> \brief Generating soil database from input file.
6 : !> \details This module provides the routines for generating the soil database for mHM from an ASCII input file.
7 : !! One routine \e read_soil_LUT reads a soil LookUpTable, performs some consistency checks and returns an initial
8 : !! soil database.
9 : !! The second routine \e generate_soil_database calculates based on the initial one the proper soil database.
10 : !> \authors Juliane Mai
11 : !> \date Dec 2012
12 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
13 : !! mHM is released under the LGPLv3+ license \license_note
14 : !> \ingroup f_mpr
15 : MODULE mo_soil_database
16 :
17 : ! This module to provide a soil database for mHM.
18 :
19 : ! Written Juliane Mai, Dec 2012
20 :
21 : use mo_kind, only : i4, dp
22 : use mo_message, only: message, error_message
23 : use mo_string_utils, only : num2str
24 : use mo_os, only : check_path_isfile
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : PUBLIC :: read_soil_LUT ! Reads the soil LUT file
31 : PUBLIC :: generate_soil_database ! Generates the soil database
32 :
33 : ! ------------------------------------------------------------------
34 :
35 : CONTAINS
36 :
37 : ! ------------------------------------------------------------------
38 :
39 : ! NAME
40 : ! read_soil_LUT
41 :
42 : ! PURPOSE
43 : !> \brief Reads the soil LUT file.
44 :
45 : !> \details Reads the soil LookUpTable file and checks for consistency.
46 :
47 : ! INTENT(IN)
48 : !> \param[in] "character(len = *) :: filename" filename of the soil LUT
49 :
50 : ! HISTORY
51 : !> \authors Juliane Mai
52 :
53 : !> \date Dec 2012
54 :
55 : ! Modifications:
56 : ! Luis Samaniego Nov 2013 - transform relation op. == -> .eq. etc
57 : ! Rohini Kumar Mar 2016 - new variables for handling different soil databases
58 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
59 :
60 13 : subroutine read_soil_LUT(filename)
61 :
62 : use mo_common_constants, only : eps_dp, nodata_dp, nodata_i4
63 : use mo_mpr_constants, only : nLCover_class
64 : use mo_mpr_file, only : usoil_database
65 : use mo_mpr_global_variables, only : HorizonDepth_mHM, iFlag_soilDB, nSoilHorizons_mHM, nSoilTypes, soilDB, &
66 : tillageDepth
67 :
68 : implicit none
69 :
70 : ! filename of the soil LUT
71 : character(len = *), intent(in) :: filename
72 :
73 : character(len = 10) :: dummy
74 :
75 : integer(i4) :: ios
76 :
77 : integer(i4) :: ii, jj, kk
78 :
79 : integer(i4) :: nR, nH, max_i
80 :
81 13 : real(dp) :: up, down
82 :
83 13 : real(dp) :: cly, snd, bd
84 :
85 13 : real(dp) :: dmin
86 :
87 : integer(i4) :: maxNumberHorizons
88 :
89 : integer(i4) :: minNumberTillHorizons, maxNumberTillHorizons
90 :
91 :
92 13 : SELECT CASE (iFlag_soilDB)
93 : ! classical mHM soil database
94 : CASE(0)
95 13 : ios = 0_i4
96 : !checking whether the file exists
97 13 : call check_path_isfile(path = filename, raise=.true.)
98 13 : open(usoil_database, file = filename, status = 'old', iostat = ios)
99 13 : read(usoil_database, *) dummy, nSoilTypes
100 13 : dummy = dummy // '' ! only to avoid warning
101 :
102 : ! allocate space
103 39 : allocate(soilDB%Id(nSoilTypes))
104 26 : allocate(soilDB%nHorizons(nSoilTypes))
105 26 : allocate(soilDB%nTillHorizons(nSoilTypes))
106 39 : allocate(soilDB%RZdepth(nSoilTypes))
107 : ! initialize
108 19188 : soilDB%Id(:) = nodata_i4
109 19188 : soilDB%nHorizons(:) = nodata_i4
110 19188 : soilDB%nTillHorizons(:) = nodata_i4
111 19188 : soilDB%RZdepth(:) = nodata_dp
112 :
113 : ! check if max soil_type id is equal to nSoilTypes
114 13 : max_i = 0_i4
115 : ! initalise total rows to read
116 13 : nR = 0_i4
117 13 : read(usoil_database, *) dummy
118 38376 : do while (.NOT. (ios .ne. 0))
119 38363 : read(usoil_database, *, IOSTAT = ios) ii, jj, up, down, cly, snd, bd
120 :
121 : ! Checks
122 38363 : if(up .ge. down) then
123 : call error_message('read_soil_LUT: ERROR occurred: Mixed horizon depths in soil type', &
124 0 : num2str(ii, '(I3)'), ' and horizon no.', num2str(jj, '(I3)'))
125 : end if
126 : if(cly .lt. 0.0_dp .OR. cly .gt. 100.0_dp .OR. &
127 : snd .lt. 0.0_dp .OR. snd .gt. 100.0_dp .OR. &
128 38363 : bd .lt. 0.0_dp .OR. bd .gt. 5.0_dp) then
129 : call error_message('read_soil_LUT: ERROR occurred: Inappropriate soil properties in soil type', &
130 0 : num2str(ii, '(I3)'), ' and horizon no.', num2str(jj, '(I3)'))
131 : end if
132 :
133 : ! initalise soil id
134 38363 : if ( ii > size(soilDB%Id) ) call error_message( &
135 : "ERROR: nSoil_Types (", num2str(size(soilDB%Id)), &
136 : ") in soil_classdefinition.txt seems to be to low! Tried to read: ", num2str(ii) &
137 0 : )
138 38363 : soilDB%Id(ii) = ii
139 38363 : soilDB%nHorizons(ii) = jj
140 38363 : if( anint(down, dp) .gt. soilDB%RZdepth(ii) ) soilDB%RZdepth(ii) = anint(down, dp)
141 :
142 38363 : nR = nR + 1_i4
143 38376 : max_i = max(max_i, ii)
144 :
145 : end do
146 :
147 13 : if ( max_i < nSoilTypes ) call error_message( &
148 0 : "ERROR: nSoil_Types (", num2str(size(soilDB%Id)), &
149 : ") in soil_classdefinition.txt seems to be to high! Highest read value: ", num2str(max_i) &
150 0 : )
151 :
152 : ! initalise minimum root zone depth among all soil types
153 19201 : dMin = minval(soilDB%RZdepth(:), soilDB%RZdepth(:) .gt. 0.0_dp)
154 :
155 : ! check the tillage depth...(if possible adjust it..)
156 13 : if(tillageDepth .gt. dMin) then
157 0 : call error_message('read_soil_LUT: ERROR occurred: ', raise=.false.)
158 0 : call error_message(' Tillage depth is greater than overall minimum total soil depth ', raise=.false.)
159 0 : call error_message(' So tillage depth should be at least', num2str(dMin, '(F7.2)'), raise=.false.)
160 0 : call error_message(' Please adjust!')
161 : end if
162 :
163 : ! insert a new tillage soil layer, only in those soil types, in which it is not present
164 13 : rewind(usoil_database)
165 13 : read(usoil_database, *) dummy
166 13 : read(usoil_database, *) dummy
167 :
168 : ! Last row is read twice so, read only upto (nR - 1)
169 38363 : do ii = 1, nR - 1
170 38350 : read(usoil_database, *) jj, nH, up, down, cly, snd, bd
171 38350 : if(tillageDepth .gt. anint(up, dp) .and. tillageDepth .lt. anint(down, dp)) then
172 19175 : soilDB%nHorizons(jj) = soilDB%nHorizons(jj) + 1_i4
173 : end if
174 :
175 : ! identify upto which soil horizon does the tillage depth goes in
176 38363 : if(tillageDepth .gt. anint(up, dp) .and. tillageDepth .le. anint(down, dp)) then
177 19175 : soilDB%nTillHorizons(jj) = nH
178 : end if
179 : end do
180 :
181 19188 : maxNumberHorizons = maxval(soilDB%nHorizons(:))
182 :
183 : ! the variables of soilType are all allocated with maximal number of Horizons although not for each of the
184 : ! nSoilTypes the array will be used
185 : ! loops for array(i,j) should be: j=1, nHorizons(i), otherwise a nodata_dp will appear
186 52 : allocate(soilDB%UD(nSoilTypes, maxNumberHorizons))
187 39 : allocate(soilDB%LD(nSoilTypes, maxNumberHorizons))
188 39 : allocate(soilDB%depth(nSoilTypes, maxNumberHorizons))
189 39 : allocate(soilDB%clay(nSoilTypes, maxNumberHorizons))
190 39 : allocate(soilDB%sand(nSoilTypes, maxNumberHorizons))
191 39 : allocate(soilDB%dbM(nSoilTypes, maxNumberHorizons))
192 :
193 57577 : soilDB%UD(:, :) = nodata_dp
194 57577 : soilDB%LD(:, :) = nodata_dp
195 57577 : soilDB%depth(:, :) = nodata_dp
196 57577 : soilDB%clay(:, :) = nodata_dp
197 57577 : soilDB%sand(:, :) = nodata_dp
198 57577 : soilDB%dbM(:, :) = nodata_dp
199 :
200 : ! allocate space for other derived soil hydraulic properties ...
201 19188 : minNumberTillHorizons = minval(soilDB%nTillHorizons(:))
202 19188 : maxNumberTillHorizons = maxval(soilDB%nTillHorizons(:))
203 :
204 65 : allocate(soilDB%thetaS_till(nSoilTypes, 1 : maxNumberTillHorizons, nLCover_class))
205 52 : allocate(soilDB%thetaS(nSoilTypes, minNumberTillHorizons + 1 : maxNumberHorizons))
206 :
207 39 : allocate(soilDB%thetaFC_till(nSoilTypes, 1 : maxNumberTillHorizons, nLCover_class))
208 39 : allocate(soilDB%thetaFC(nSoilTypes, minNumberTillHorizons + 1 : maxNumberHorizons))
209 :
210 39 : allocate(soilDB%thetaPW_till(nSoilTypes, 1 : maxNumberTillHorizons, nLCover_class))
211 39 : allocate(soilDB%thetaPW(nSoilTypes, minNumberTillHorizons + 1 : maxNumberHorizons))
212 :
213 39 : allocate(soilDB%Db(nSoilTypes, 1 : maxNumberTillHorizons, nLCover_class))
214 :
215 65 : allocate(soilDB%Ks(nSoilTypes, 1 : maxNumberHorizons, nLCover_class))
216 :
217 : ! Initialize with default nodata_dp
218 57616 : soilDB%thetaS_till(:, :, :) = nodata_dp
219 38389 : soilDB%thetaS(:, :) = nodata_dp
220 :
221 57616 : soilDB%thetaFC_till(:, :, :) = nodata_dp
222 38389 : soilDB%thetaFC(:, :) = nodata_dp
223 :
224 57616 : soilDB%thetaPW_till(:, :, :) = nodata_dp
225 38389 : soilDB%thetaPW(:, :) = nodata_dp
226 :
227 57616 : soilDB%Db(:, :, :) = nodata_dp
228 172744 : soilDB%Ks(:, :, :) = nodata_dp
229 :
230 : ! Read again soil properties from the data base
231 13 : rewind(usoil_database)
232 13 : read(usoil_database, *) dummy
233 13 : read(usoil_database, *) dummy
234 :
235 : ! Last row is read twice so, read only upto (nR - 1)
236 13 : kk = 0_i4
237 38363 : do ii = 1, nR - 1_i4
238 38350 : read(usoil_database, *) jj, nH, up, down, cly, snd, bd
239 :
240 : ! to avoid numerical errors in PTF
241 38350 : if(cly .lt. 1.0_dp) cly = 1.0_dp
242 38350 : if(snd .lt. 1.0_dp) snd = 1.0_dp
243 :
244 : ! Physical consistency
245 38350 : if((cly + snd) .gt. 100.0_dp) then
246 0 : cly = cly / (cly + snd)
247 0 : snd = snd / (cly + snd)
248 : end if
249 :
250 : ! check for an extra tillage horizon (if not exists create a layer)
251 38350 : if(tillageDepth .gt. anint(up, dp) .and. tillageDepth .lt. anint(down, dp)) then
252 :
253 19175 : soilDB%UD(jj, nH) = anint(up, dp)
254 19175 : soilDB%LD(jj, nH) = tillageDepth
255 19175 : soilDB%depth(jj, nH) = soilDB%LD(jj, nH) - soilDB%UD(jj, nH)
256 19175 : soilDB%clay(jj, nH) = cly
257 19175 : soilDB%sand(jj, nH) = snd
258 19175 : soilDB%dbM(jj, nH) = bd
259 :
260 : ! initalise the upper depth to new value... and the increment counter
261 19175 : up = tillageDepth
262 : kk = 1_i4
263 : end if
264 :
265 : ! increment nH by one once it encounter the upper loop..
266 38350 : if(kk .eq. 1_i4) nH = nH + 1_i4
267 :
268 38350 : soilDB%UD(jj, nH) = anint(up, dp)
269 38350 : soilDB%LD(jj, nH) = anint(down, dp)
270 38350 : soilDB%depth(jj, nH) = soilDB%LD(jj, nH) - soilDB%UD(jj, nH)
271 38350 : soilDB%clay(jj, nH) = cly
272 38350 : soilDB%sand(jj, nH) = snd
273 38350 : soilDB%dbM(jj, nH) = bd
274 :
275 : ! check for number of soil horizons...
276 38350 : if(nH .gt. soilDB%nHorizons(jj)) then
277 0 : call error_message('read_soil_LUT: ERROR occurred: ', raise=.false.)
278 0 : call error_message(' There is something wrong in allocating horizons in soil data base.', raise=.false.)
279 0 : call error_message(' Please check in code !')
280 : end if
281 :
282 : ! initalise the increment counter to zero
283 38363 : if(nH .eq. soilDB%nHorizons(jj)) kk = 0_i4
284 :
285 : end do
286 13 : close(usoil_database)
287 :
288 : ! soil database for the horizon specific case
289 : CASE(1)
290 :
291 0 : allocate(soilDB%RZdepth(1))
292 0 : allocate(soilDB%UD(1, 1))
293 0 : allocate(soilDB%LD(1, 1))
294 0 : allocate(soilDB%depth(1, 1))
295 0 : allocate(soilDB%thetaS_till(1, 1, 1))
296 0 : allocate(soilDB%thetaS(1, 1))
297 0 : allocate(soilDB%thetaFC_till(1, 1, 1))
298 0 : allocate(soilDB%thetaFC(1, 1))
299 0 : allocate(soilDB%thetaPW_till(1, 1, 1))
300 0 : allocate(soilDB%thetaPW(1, 1))
301 0 : allocate(soilDB%Db(1, 1, 1))
302 0 : allocate(soilDB%Ks(1, 1, 1))
303 :
304 :
305 :
306 : !checking whether the file exists
307 0 : call check_path_isfile(path = filename, raise=.true.)
308 0 : open(usoil_database, file = filename, status = 'old', action = 'read')
309 0 : read(usoil_database, *) dummy, nSoilTypes
310 0 : dummy = dummy // '' ! only to avoid warning
311 0 : allocate(soilDB%Id (nSoilTypes))
312 0 : allocate(soilDB%clay(nSoilTypes, 1))
313 0 : allocate(soilDB%sand(nSoilTypes, 1))
314 0 : allocate(soilDB%dbM (nSoilTypes, 1))
315 0 : soilDB%clay(:, :) = nodata_dp
316 0 : soilDB%sand(:, :) = nodata_dp
317 0 : soilDB%dbM (:, :) = nodata_dp
318 0 : read(usoil_database, *) dummy
319 0 : do kk = 1, nSoilTypes
320 0 : read(usoil_database, *, iostat=ios) jj, cly, snd, bd
321 0 : if ( ios /= 0 ) call error_message( &
322 : "ERROR: nSoil_Types (", num2str(nSoilTypes), ") in soil_classdefinition_iFlag_soilDB_1.txt ", &
323 : "seems to be higher than available soil types!" &
324 0 : )
325 : ! to avoid numerical errors in PTF
326 0 : if(cly .lt. 1.0_dp) cly = 1.0_dp
327 0 : if(snd .lt. 1.0_dp) snd = 1.0_dp
328 0 : soilDB%Id(kk) = jj
329 0 : soilDB%clay(kk, 1) = cly
330 0 : soilDB%sand(kk, 1) = snd
331 0 : soilDB%dbM (kk, 1) = bd
332 : end do
333 0 : close(usoil_database)
334 :
335 : ! assign up to which horizon layer a soil is treated as tillage layer
336 : ! since our horizon information is well defined for modeling too
337 : ! this information is uniform across all soils/modeling cells
338 : ! for compatibility with iFlag_soilDB == 0, assign
339 : ! both nHorizons & tillage horizons
340 0 : allocate(soilDB%nHorizons(1))
341 0 : allocate(soilDB%nTillHorizons(1))
342 0 : soilDB%nHorizons(:) = nSoilHorizons_mHM
343 0 : soilDB%nTillHorizons(:) = -9
344 0 : do kk = 1, nSoilHorizons_mHM
345 0 : if(abs(HorizonDepth_mHM(kk) - tillageDepth) .lt. eps_dp) then
346 0 : soilDB%nTillHorizons(1) = kk
347 : end if
348 : end do
349 :
350 : ! check
351 0 : if(soilDB%nTillHorizons(1) .eq. -9) then
352 : ! rarely could happen *** since this is checked in reading of horizons depths only
353 : ! but is checked here for double confirmation
354 0 : call error_message('***ERROR: specification of tillage depths is not confirming', raise=.false.)
355 0 : call error_message(' with given depths of soil horizons to be modeled.')
356 : else
357 0 : call message()
358 0 : call message('Tillage layers: the tillage horizons are modelled ')
359 0 : call message(' upto mHM layers: ', trim(num2str(soilDB%nTillHorizons(1))))
360 : end if
361 :
362 : CASE DEFAULT
363 0 : call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
364 :
365 : END SELECT
366 :
367 13 : end subroutine read_soil_LUT
368 :
369 : ! ------------------------------------------------------------------
370 :
371 : ! NAME
372 : ! generate_soil_database
373 :
374 : ! PURPOSE
375 : !> \brief Generates soil database.
376 :
377 : !> \details Calculates the proper soil database using the initialized soil database from read_soil_LUT.
378 :
379 : ! HISTORY
380 : !> \authors Juliane Mai
381 :
382 : !> \date Dec 2012
383 :
384 : ! Modifications:
385 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
386 :
387 13 : subroutine generate_soil_database
388 :
389 13 : use mo_common_constants, only : nodata_dp, nodata_i4
390 : use mo_mpr_global_variables, only : HorizonDepth_mHM, iFlag_soilDB, nSoilHorizons_mHM, nSoilTypes, soilDB
391 :
392 : implicit none
393 :
394 : integer(i4) :: ii, jj, kk
395 :
396 13 : real(dp) :: dmin
397 :
398 13 : real(dp) :: dpth_f, dpth_t
399 :
400 : integer(i4) :: layer_f, layer_t
401 :
402 : real(dp), parameter :: small = 0.000001_dp
403 :
404 : ! [mm] soil depth accuracy
405 : real(dp), parameter :: soil_dAccuracy = 0.5_dp
406 :
407 :
408 13 : SELECT CASE (iFlag_soilDB)
409 : ! classical mHM soil database
410 : CASE(0)
411 : ! initalise minimum root zone depth among all soil types
412 19201 : dMin = minval(soilDB%RZdepth(:), soilDB%RZdepth(:) > 0.0_dp)
413 :
414 : ! check
415 13 : if (HorizonDepth_mHM(nSoilHorizons_mHM - 1) .ge. dMin) then
416 0 : call error_message('generate_soil_database: ERROR occurred: ', raise=.false.)
417 0 : call error_message(' The depth of soil Horizons provided for modelling is not appropriate', raise=.false.)
418 0 : call error_message(' The global minimum of total soil horizon depth among all soil type is ', num2str(dMin, '(F7.2)'), raise=.false.)
419 0 : call error_message(' Adjust your modeling soil horizon depth in this range', raise=.false.)
420 0 : call error_message(' OR Increase the soil depth in data base for only those soil types', raise=.false.)
421 0 : call error_message(' whose total depth is smaller than your given modeling depth.')
422 : end if
423 :
424 : ! allocate and initalise depth weight
425 19240 : allocate(soilDB%Wd(nSoilTypes, nSoilHorizons_mHM, maxval(soilDB%nHorizons(:))))
426 115180 : soilDB%Wd(:, :, :) = 0.0_dp
427 :
428 : ! Process further to estimate weight of each horizons
429 : ! weightage according to soil depths
430 19188 : do ii = 1, nSoilTypes
431 28302300 : soilDB%Wd(ii, :, soilDB%nHorizons(ii) + 1_i4 : maxval(soilDB%nHorizons(:))) = nodata_dp
432 :
433 : ! initalise last horizon depth to model w.r.t to surface
434 : ! NOTE:: it is different for each soil
435 19175 : HorizonDepth_mHM(nSoilHorizons_mHM) = soilDB%RZdepth(ii)
436 :
437 : ! Estimate soil properties for each modeling layers
438 57538 : do jj = 1, nSoilHorizons_mHM
439 :
440 : ! modeling depth ( **from --> to ** )
441 : ! take into account the depth accuracy [0.5mm, defined in module..]
442 38350 : dpth_f = 0.0_dp
443 38350 : if(jj .ne. 1_i4) dpth_f = HorizonDepth_mHM(jj - 1)
444 38350 : dpth_t = HorizonDepth_mHM(jj) - soil_dAccuracy
445 :
446 : ! identify to which layer of batabase this mHM horizon lies
447 38350 : layer_f = nodata_i4
448 38350 : layer_t = nodata_i4
449 153400 : do kk = 1, soilDB%nHorizons(ii)
450 115050 : if(dpth_f .ge. soilDB%UD(ii, kk) .and. dpth_f .le. (soilDB%LD(ii, kk) - soil_dAccuracy)) layer_f = kk
451 153400 : if(dpth_t .ge. soilDB%UD(ii, kk) .and. dpth_t .le. (soilDB%LD(ii, kk) - soil_dAccuracy)) layer_t = kk
452 : end do
453 :
454 : ! Check
455 38350 : if(layer_f .le. 0_i4 .or. layer_t .le. 0_i4) then
456 0 : call error_message('generate_soil_database: ERROR occurred: ', raise=.false.)
457 0 : call error_message(' Horizon depths to model do not lie in database for soil type', num2str(ii, '(I3)'), raise=.false.)
458 0 : call error_message(' Please check!')
459 : end if
460 38350 : if(layer_f .gt. layer_t) then
461 0 : call error_message('generate_soil_database: ERROR occurred: ', raise=.false.)
462 0 : call error_message(' Something is wrong in assignment of modeling soil horizons or', raise=.false.)
463 0 : call error_message(' database of soil type ', num2str(ii, '(I3)'), raise=.false.)
464 0 : call error_message(' Please check!')
465 : end if
466 :
467 : ! iF modeling depth of a given horizon falls in a same soil layer
468 57525 : if(layer_f .eq. layer_t) then
469 19175 : soilDB%Wd(ii, jj, layer_f) = 1.0_dp
470 :
471 : ! else estimate depth weightage...
472 : else
473 :
474 : ! for starting layer...
475 19175 : soilDB%Wd(ii, jj, layer_f) = soilDB%LD(ii, layer_f) - dpth_f
476 :
477 : ! for ending layer
478 19175 : soilDB%Wd(ii, jj, layer_t) = (dpth_t + soil_dAccuracy) - soilDB%UD(ii, layer_t)
479 :
480 : ! other intermediate layers weight, if exit
481 19175 : if(layer_t - layer_f .gt. 1_i4) then
482 0 : do kk = layer_f + 1, layer_t - 1
483 0 : soilDB%Wd(ii, jj, kk) = soilDB%LD(ii, kk) - soilDB%UD(ii, kk)
484 : end do
485 : end if
486 :
487 : ! Estimate depth weightage
488 19175 : if(jj .ne. 1_i4) then
489 19175 : soilDB%Wd(ii, jj, 1 : soilDB%nHorizons(ii)) = soilDB%Wd(ii, jj, 1 : soilDB%nHorizons(ii)) / &
490 95875 : (HorizonDepth_mHM(jj) - HorizonDepth_mHM(jj - 1_i4))
491 : else
492 0 : soilDB%Wd(ii, jj, 1 : soilDB%nHorizons(ii)) = soilDB%Wd(ii, jj, 1 : soilDB%nHorizons(ii)) / &
493 0 : HorizonDepth_mHM(jj)
494 : end if
495 :
496 : ! Check (small margin for numerical errors)
497 153400 : if(sum(soilDB%Wd(ii, jj, :), soilDB%Wd(ii, jj, :) .gt. 0.0_dp) .le. 1.0_dp - small .or. &
498 : sum(soilDB%Wd(ii, jj, :), soilDB%Wd(ii, jj, :) .gt. 0.0_dp) .ge. 1.0_dp + small) then
499 0 : call error_message('generate_soil_database: ERROR occurred: ', raise=.false.)
500 0 : call error_message(' Weight assigned for each soil horizons are not correct.', raise=.false.)
501 0 : call error_message(' Please check!')
502 : end if
503 :
504 : end if
505 :
506 : end do
507 : end do
508 : ! soil database for the horizon specific case
509 : CASE(1)
510 : ! right now nothing is done here
511 : ! *** reserved for future changes
512 0 : allocate(soilDB%Wd(1,1,1))
513 :
514 : CASE DEFAULT
515 13 : call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
516 : END SELECT
517 :
518 13 : end subroutine generate_soil_database
519 :
520 : END MODULE mo_soil_database
|
