Line data Source code
1 : !> \file mo_write_fluxes_states.f90
2 : !> \brief \copybrief mo_write_fluxes_states
3 : !> \details \copydetails mo_write_fluxes_states
4 :
5 : !> \brief Creates NetCDF output for different fluxes and state variables of mHM.
6 : !> \details NetCDF is first initialized and later on variables are put to the NetCDF.
7 : !> \changelog
8 : !! - David Schaefer Aug 2015
9 : !! - major rewrite
10 : !! - O. Rakovec, R. Kumar Nov 2017
11 : !! - added project description for the netcdf outputs
12 : !> \authors Matthias Zink
13 : !> \date Apr 2013
14 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
15 : !! mHM is released under the LGPLv3+ license \license_note
16 : !> \ingroup f_mhm
17 : module mo_write_fluxes_states
18 :
19 : use mo_nc_output, only: OutputDataset, OutputVariable, set_attributes, data_dims, data_dtype
20 : use mo_kind, only : i4, dp
21 : use mo_string_utils, only : num2str
22 : use mo_common_constants, only : nodata_dp
23 : use mo_netcdf, only : NcDataset
24 : use mo_global_variables, only : &
25 : output_deflate_level, output_double_precision, output_time_reference, timeStep_model_outputs, outputFlxState
26 : use mo_common_mHM_mRM_variables, only: timeStep
27 : use mo_mpr_global_variables, only : nSoilHorizons_mHM
28 : use mo_common_variables, only : iFlag_cordinate_sys, level1
29 : use mo_file, only : file_mhm_output
30 :
31 : implicit none
32 :
33 : contains
34 :
35 : !> \brief Initialize mHM OutputDataset
36 : !> \details Create and initialize the output file. If new a new output
37 : !! variable needs to be written, this is the first of two
38 : !! procedures to change (second: updateDataset)
39 : !!
40 : !> \changelog
41 : !! - Robert Schweppe Jun 2018
42 : !! - refactoring and reformatting
43 : !! - Pallav Shrestha Mar 2020
44 : !! - iFlag_cordinate_sys based dimensions (dims)
45 : !!
46 : !> \return type(OutputDataset)
47 : !> \authors Matthias Zink
48 : !> \date Apr 2013
49 15 : function mHM_OutputDataset(iDomain, mask1) result(out)
50 : implicit none
51 :
52 : !> domain id
53 : integer(i4), intent(in) :: iDomain
54 : !> L1 mask to reconstruct the data
55 : logical, pointer, intent(in), dimension(:, :) :: mask1
56 :
57 : type(OutputDataset) :: out
58 :
59 : integer(i4) :: ii, nn, nCells
60 : character(3) :: dtype
61 : character(16), dimension(3) :: dims
62 : character(16) :: unit
63 2550 : type(OutputVariable), dimension(size(outputFlxState) * nSoilHorizons_mHM) :: tmpvars
64 :
65 15 : nCells = level1(iDomain)%nCells
66 :
67 : out = OutputDataset( &
68 : iDomain=iDomain, &
69 : level=level1, &
70 : file_name=file_mhm_output, &
71 : double_precision=output_double_precision, &
72 : outputs_frequence=timeStep_model_outputs, &
73 : time_reference=output_time_reference &
74 15 : )
75 :
76 15 : dtype = data_dtype(output_double_precision)
77 15 : unit = fluxesUnit(iDomain)
78 15 : dims = data_dims()
79 :
80 15 : ii = 0
81 :
82 15 : if (outputFlxState(1)) then
83 0 : ii = ii + 1
84 0 : tmpvars(ii) = OutputVariable(&
85 0 : out%nc, "interception", dtype, dims, nCells, mask1, output_deflate_level, .true.)
86 : call set_attributes(&
87 : tmpvars(ii)%nc, "canopy interception storage", "mm", output_double_precision)
88 : end if
89 :
90 15 : if (outputFlxState(2)) then
91 0 : ii = ii + 1
92 0 : tmpvars(ii) = OutputVariable(&
93 0 : out%nc, "snowpack", dtype, dims, nCells, mask1, output_deflate_level, .true.)
94 : call set_attributes(tmpvars(ii)%nc, "depth of snowpack", "mm", output_double_precision)
95 : end if
96 :
97 15 : if (outputFlxState(3)) then
98 45 : do nn = 1, nSoilHorizons_mHM
99 30 : ii = ii + 1
100 30 : tmpvars(ii) = OutputVariable(out%nc, "SWC_L" // trim(num2str(nn, '(i2.2)')), &
101 60 : dtype, dims, nCells, mask1, output_deflate_level, .true.)
102 : call set_attributes(tmpvars(ii)%nc, &
103 45 : 'soil water content of soil layer' // trim(num2str(nn)), "mm", output_double_precision)
104 : end do
105 : end if
106 :
107 15 : if (outputFlxState(4)) then
108 0 : do nn = 1, nSoilHorizons_mHM
109 0 : ii = ii + 1
110 0 : tmpvars(ii) = OutputVariable(out%nc, "SM_L" // trim(num2str(nn, '(i2.2)')), &
111 0 : dtype, dims, nCells, mask1, output_deflate_level, .true.)
112 : call set_attributes(tmpvars(ii)%nc, &
113 0 : 'volumetric soil moisture of soil layer' // trim(num2str(nn)), "mm mm-1", output_double_precision)
114 : end do
115 : end if
116 :
117 15 : if (outputFlxState(5)) then
118 0 : ii = ii + 1
119 0 : tmpvars(ii) = OutputVariable(&
120 0 : out%nc, "SM_Lall", dtype, dims, nCells, mask1, output_deflate_level, .true.)
121 : call set_attributes(&
122 : tmpvars(ii)%nc, "average soil moisture over all layers", "mm mm-1", output_double_precision)
123 : end if
124 :
125 15 : if (outputFlxState(6)) then
126 0 : ii = ii + 1
127 0 : tmpvars(ii) = OutputVariable(&
128 0 : out%nc, "sealedSTW", dtype, dims, nCells, mask1, output_deflate_level, .true.)
129 : call set_attributes(&
130 : tmpvars(ii)%nc, "reservoir of sealed areas (sealedSTW)", "mm", output_double_precision)
131 : end if
132 :
133 15 : if (outputFlxState(7)) then
134 0 : ii = ii + 1
135 0 : tmpvars(ii) = OutputVariable(&
136 0 : out%nc, "unsatSTW", dtype, dims, nCells, mask1, output_deflate_level, .true.)
137 : call set_attributes(&
138 : tmpvars(ii)%nc, "reservoir of unsaturated zone", "mm", output_double_precision)
139 : end if
140 :
141 15 : if (outputFlxState(8)) then
142 0 : ii = ii + 1
143 0 : tmpvars(ii) = OutputVariable(&
144 0 : out%nc, "satSTW", dtype, dims, nCells, mask1, output_deflate_level, .true.)
145 : call set_attributes(&
146 : tmpvars(ii)%nc, "water level in groundwater reservoir", "mm", output_double_precision)
147 : end if
148 :
149 15 : if (outputFlxState(18)) then
150 0 : ii = ii + 1
151 0 : tmpvars(ii) = OutputVariable(&
152 0 : out%nc, "neutrons", dtype, dims, nCells, mask1, output_deflate_level, .true.)
153 : call set_attributes(&
154 : tmpvars(ii)%nc, "ground albedo neutrons", "cph", output_double_precision)
155 : end if
156 :
157 15 : if (outputFlxState(9)) then
158 2 : ii = ii + 1
159 2 : tmpvars(ii) = OutputVariable(&
160 4 : out%nc, "PET", dtype, dims, nCells, mask1, output_deflate_level)
161 : call set_attributes(&
162 2 : tmpvars(ii)%nc, "potential Evapotranspiration", trim(unit), output_double_precision)
163 : end if
164 :
165 15 : if (outputFlxState(10)) then
166 5 : ii = ii + 1
167 5 : tmpvars(ii) = OutputVariable(&
168 10 : out%nc, "aET", dtype, dims, nCells, mask1, output_deflate_level)
169 : call set_attributes(&
170 5 : tmpvars(ii)%nc, "actual Evapotranspiration", trim(unit), output_double_precision)
171 : end if
172 :
173 15 : if (outputFlxState(11)) then
174 15 : ii = ii + 1
175 15 : tmpvars(ii) = OutputVariable(&
176 30 : out%nc, "Q", dtype, dims, nCells, mask1, output_deflate_level)
177 : call set_attributes(&
178 15 : tmpvars(ii)%nc, "total runoff generated by every cell", trim(unit), output_double_precision)
179 : end if
180 :
181 15 : if (outputFlxState(12)) then
182 0 : ii = ii + 1
183 0 : tmpvars(ii) = OutputVariable(&
184 0 : out%nc, "QD", dtype, dims, nCells, mask1, output_deflate_level)
185 : call set_attributes(tmpvars(ii)%nc, &
186 0 : "direct runoff generated by every cell (runoffSeal)", trim(unit), output_double_precision)
187 : end if
188 :
189 15 : if (outputFlxState(13)) then
190 0 : ii = ii + 1
191 0 : tmpvars(ii) = OutputVariable(&
192 0 : out%nc, "QIf", dtype, dims, nCells, mask1, output_deflate_level)
193 : call set_attributes(tmpvars(ii)%nc, &
194 0 : "fast interflow generated by every cell (fastRunoff)", trim(unit), output_double_precision)
195 : end if
196 :
197 15 : if (outputFlxState(14)) then
198 0 : ii = ii + 1
199 0 : tmpvars(ii) = OutputVariable(&
200 0 : out%nc, "QIs", dtype, dims, nCells, mask1, output_deflate_level)
201 : call set_attributes(tmpvars(ii)%nc, &
202 0 : "slow interflow generated by every cell (slowRunoff)", trim(unit), output_double_precision)
203 : end if
204 :
205 15 : if (outputFlxState(15)) then
206 0 : ii = ii + 1
207 0 : tmpvars(ii) = OutputVariable(&
208 0 : out%nc, "QB", dtype, dims, nCells, mask1, output_deflate_level)
209 : call set_attributes(&
210 0 : tmpvars(ii)%nc, "baseflow generated by every cell", trim(unit), output_double_precision)
211 : end if
212 :
213 15 : if (outputFlxState(16)) then
214 2 : ii = ii + 1
215 2 : tmpvars(ii) = OutputVariable(&
216 4 : out%nc, "recharge", dtype, dims, nCells, mask1, output_deflate_level)
217 : call set_attributes(&
218 2 : tmpvars(ii)%nc, "groundwater recharge", trim(unit), output_double_precision)
219 : end if
220 :
221 15 : if (outputFlxState(17)) then
222 0 : do nn = 1, nSoilHorizons_mHM
223 0 : ii = ii + 1
224 0 : tmpvars(ii) = OutputVariable(&
225 : out%nc, "soil_infil_L" // trim(num2str(nn, '(i2.2)')), &
226 0 : dtype, dims, nCells, mask1, output_deflate_level)
227 : call set_attributes(tmpvars(ii)%nc, &
228 0 : "infiltration flux from soil layer" // trim(num2str(nn)), unit, output_double_precision)
229 : end do
230 : end if
231 :
232 15 : if (outputFlxState(19)) then
233 0 : do nn = 1, nSoilHorizons_mHM
234 0 : ii = ii + 1
235 0 : tmpvars(ii) = OutputVariable(&
236 : out%nc, "aET_L" // trim(num2str(nn, '(i2.2)')), &
237 0 : dtype, dims, nCells, mask1, output_deflate_level)
238 : call set_attributes(tmpvars(ii)%nc, &
239 : 'actual Evapotranspiration from soil layer' // trim(num2str(nn)), &
240 0 : "mm " // trim(unit), output_double_precision)
241 : end do
242 : end if
243 :
244 15 : if (outputFlxState(20)) then
245 0 : ii = ii + 1
246 0 : tmpvars(ii) = OutputVariable(&
247 0 : out%nc, "preEffect", dtype, dims, nCells, mask1, output_deflate_level)
248 : call set_attributes(&
249 0 : tmpvars(ii)%nc, "effective precipitation", trim(unit), output_double_precision)
250 : end if
251 :
252 15 : if (outputFlxState(21)) then
253 0 : ii = ii + 1
254 0 : tmpvars(ii) = OutputVariable(&
255 0 : out%nc, "Qsm", dtype, dims, nCells, mask1, output_deflate_level)
256 : call set_attributes(&
257 : tmpvars(ii)%nc, "Average liquid water generated from solid to liquid phase change in the snow", &
258 0 : trim(unit), output_double_precision)
259 : end if
260 :
261 99 : allocate(out%vars(ii))
262 84 : out%vars = tmpvars(1 : ii)
263 :
264 15 : end function mHM_OutputDataset
265 :
266 : !> \brief Update all variables.
267 : !> \details Call the type bound procedure updateVariable for
268 : !> all output variables. If a new output
269 : !> variable needs to be written, this is the second
270 : !> of two procedures to change (first: newOutputDataset)
271 : !!
272 : !> \changelog
273 : !! - L. Samaniego et al. Dec 2013
274 : !! - nullify pointer Matthias Zink, Feb. 2014
275 : !! - added aditional output: pet V. Prykhodk, J. Mai, Nov. 2014
276 : !! - adding new variable infilSoil
277 : !! - case 16 David Schaefer , Jun. 2015
278 : !! - major rewrite
279 : !! - Robert Schweppe Jun 2018
280 : !! - refactoring and reformatting
281 : !> \authors Matthias Zink
282 : !> \date Apr 2013
283 263040 : subroutine mHM_updateDataset(nc_mhm, sidx, eidx, L1_fSealed, L1_fNotSealed, L1_inter, L1_snowPack, L1_soilMoist, &
284 131520 : L1_soilMoistSat, L1_sealSTW, L1_unsatSTW, L1_satSTW, L1_neutrons, L1_pet, L1_aETSoil, &
285 131520 : L1_aETCanopy, L1_aETSealed, L1_total_runoff, L1_runoffSeal, L1_fastRunoff, L1_slowRunoff, &
286 131520 : L1_baseflow, L1_percol, L1_infilSoil, L1_preEffect, L1_melt)
287 :
288 : implicit none
289 :
290 : class(OutputDataset), intent(inout), target :: nc_mhm
291 :
292 : integer(i4), intent(in) :: sidx !< start index
293 : integer(i4), intent(in) :: eidx !< end index
294 : real(dp), intent(in), dimension(:) :: L1_fSealed
295 : real(dp), intent(in), dimension(:) :: L1_fNotSealed
296 : real(dp), intent(in), dimension(:) :: L1_inter
297 : real(dp), intent(in), dimension(:) :: L1_snowPack
298 : real(dp), intent(in), dimension(:, :) :: L1_soilMoist
299 : real(dp), intent(in), dimension(:, :) :: L1_soilMoistSat
300 : real(dp), intent(in), dimension(:) :: L1_sealSTW
301 : real(dp), intent(in), dimension(:) :: L1_unsatSTW
302 : real(dp), intent(in), dimension(:) :: L1_satSTW
303 : real(dp), intent(in), dimension(:) :: L1_neutrons
304 : real(dp), intent(in), dimension(:) :: L1_pet
305 : real(dp), intent(in), dimension(:, :) :: L1_aETSoil
306 : real(dp), intent(in), dimension(:) :: L1_aETCanopy
307 : real(dp), intent(in), dimension(:) :: L1_aETSealed
308 : real(dp), intent(in), dimension(:) :: L1_total_runoff
309 : real(dp), intent(in), dimension(:) :: L1_runoffSeal
310 : real(dp), intent(in), dimension(:) :: L1_fastRunoff
311 : real(dp), intent(in), dimension(:) :: L1_slowRunoff
312 : real(dp), intent(in), dimension(:) :: L1_baseflow
313 : real(dp), intent(in), dimension(:) :: L1_percol
314 : real(dp), intent(in), dimension(:, :) :: L1_infilSoil
315 : real(dp), intent(in), dimension(:) :: L1_preEffect
316 : real(dp), intent(in), dimension(:) :: L1_melt
317 :
318 131520 : type(OutputVariable), pointer, dimension(:) :: vars
319 : integer(i4) :: ii, nn
320 :
321 131520 : ii = 0
322 131520 : vars => nc_mhm%vars
323 :
324 131520 : if (outputFlxState(1)) then
325 0 : ii = ii + 1
326 0 : call vars(ii)%updateVariable(L1_inter(sidx : eidx))
327 : end if
328 :
329 131520 : if (outputFlxState(2)) then
330 0 : ii = ii + 1
331 0 : call vars(ii)%updateVariable(L1_snowPack(sidx : eidx))
332 : end if
333 :
334 131520 : if (outputFlxState(3)) then
335 394560 : do nn = 1, nSoilHorizons_mHM
336 263040 : ii = ii + 1
337 394560 : call vars(ii)%updateVariable(L1_soilMoist(sidx : eidx, nn))
338 : end do
339 : end if
340 :
341 131520 : if (outputFlxState(4)) then
342 0 : do nn = 1, nSoilHorizons_mHM
343 0 : ii = ii + 1
344 0 : call vars(ii)%updateVariable(L1_soilMoist(sidx : eidx, nn) / L1_soilMoistSat(sidx : eidx, nn))
345 : end do
346 : end if
347 :
348 131520 : if (outputFlxState(5)) then
349 0 : ii = ii + 1
350 0 : call vars(ii)%updateVariable(sum(L1_soilMoist(sidx : eidx, :), dim = 2) / sum(L1_soilMoistSat(sidx : eidx, :), dim = 2))
351 : end if
352 :
353 131520 : if (outputFlxState(6)) then
354 0 : ii = ii + 1
355 0 : call vars(ii)%updateVariable(L1_sealSTW(sidx : eidx))
356 : end if
357 :
358 131520 : if (outputFlxState(7)) then
359 0 : ii = ii + 1
360 0 : call vars(ii)%updateVariable(L1_unsatSTW(sidx : eidx))
361 : end if
362 :
363 131520 : if (outputFlxState(8)) then
364 0 : ii = ii + 1
365 0 : call vars(ii)%updateVariable(L1_satSTW(sidx : eidx))
366 : end if
367 :
368 131520 : if (outputFlxState(18)) then
369 0 : ii = ii + 1
370 0 : call vars(ii)%updateVariable(L1_neutrons(sidx : eidx))
371 : end if
372 :
373 131520 : if (outputFlxState(9)) then
374 17520 : ii = ii + 1
375 17520 : call vars(ii)%updateVariable(L1_pet(sidx : eidx))
376 : end if
377 :
378 131520 : if (outputFlxState(10)) then
379 43800 : ii = ii + 1
380 0 : call vars(ii)%updateVariable(sum(L1_aETSoil(sidx : eidx, :), dim = 2) * L1_fNotSealed(sidx : eidx) &
381 4511400 : + L1_aETCanopy(sidx : eidx) + L1_aETSealed(sidx : eidx) * L1_fSealed(sidx : eidx))
382 : end if
383 :
384 131520 : if (outputFlxState(11)) then
385 131520 : ii = ii + 1
386 131520 : call vars(ii)%updateVariable(L1_total_runoff(sidx : eidx))
387 : end if
388 :
389 131520 : if (outputFlxState(12)) then
390 0 : ii = ii + 1
391 0 : call vars(ii)%updateVariable(L1_runoffSeal(sidx : eidx) * L1_fSealed(sidx : eidx))
392 : end if
393 :
394 131520 : if (outputFlxState(13)) then
395 0 : ii = ii + 1
396 0 : call vars(ii)%updateVariable(L1_fastRunoff(sidx : eidx) * L1_fNotSealed(sidx : eidx))
397 : end if
398 :
399 131520 : if (outputFlxState(14)) then
400 0 : ii = ii + 1
401 0 : call vars(ii)%updateVariable(L1_slowRunoff(sidx : eidx) * L1_fNotSealed(sidx : eidx))
402 : end if
403 :
404 131520 : if (outputFlxState(15)) then
405 0 : ii = ii + 1
406 0 : call vars(ii)%updateVariable(L1_baseflow(sidx : eidx) * L1_fNotSealed(sidx : eidx))
407 : end if
408 :
409 131520 : if (outputFlxState(16)) then
410 17520 : ii = ii + 1
411 613200 : call vars(ii)%updateVariable(L1_percol(sidx : eidx) * L1_fNotSealed(sidx : eidx))
412 : end if
413 :
414 131520 : if (outputFlxState(17)) then
415 0 : do nn = 1, nSoilHorizons_mHM
416 0 : ii = ii + 1
417 0 : call vars(ii)%updateVariable(L1_infilSoil(sidx : eidx, nn) * L1_fNotSealed(sidx : eidx))
418 : end do
419 : end if
420 :
421 131520 : if (outputFlxState(19)) then
422 0 : do nn = 1, nSoilHorizons_mHM
423 0 : ii = ii + 1
424 0 : call vars(ii)%updateVariable(L1_aETSoil(sidx : eidx, nn) * L1_fNotSealed(sidx : eidx))
425 : end do
426 : end if
427 :
428 131520 : if (outputFlxState(20)) then
429 0 : ii = ii + 1
430 0 : call vars(ii)%updateVariable(L1_preEffect(sidx : eidx))
431 : end if
432 :
433 131520 : if (outputFlxState(21)) then
434 0 : ii = ii + 1
435 0 : call vars(ii)%updateVariable(L1_melt(sidx : eidx))
436 : end if
437 :
438 131535 : end subroutine mHM_updateDataset
439 :
440 : !> \brief Generate a unit string
441 : !> \details Generate the unit string for the output variable netcdf attribute based on modeling timestep
442 : !!
443 : !> \changelog
444 : !! - Robert Schweppe Jun 2018
445 : !! - refactoring and reformatting
446 : !!
447 : !> \return character(16)
448 : !> \authors David Schaefer
449 : !> \date June 2015
450 15 : function fluxesUnit(iDomain)
451 :
452 131520 : use mo_common_mhm_mrm_variables, only : nTstepDay, simPer
453 :
454 : implicit none
455 :
456 : integer(i4), intent(in) :: iDomain !< selected domain
457 :
458 : character(16) :: fluxesUnit
459 15 : real(dp) :: ntsteps
460 :
461 15 : if (timestep * timestep_model_outputs .eq. 1) then
462 0 : fluxesUnit = 'mm h-1'
463 15 : else if (timestep_model_outputs > 1) then
464 3 : fluxesUnit = 'mm ' // trim(adjustl(num2str(timestep))) // 'h-1'
465 12 : else if (timestep_model_outputs .eq. 0) then
466 12 : ntsteps = (simPer(iDomain)%julEnd - simPer(iDomain)%julStart + 1) * nTstepDay
467 12 : fluxesUnit = 'mm ' // trim(adjustl(num2str(nint(ntsteps)))) // 'h-1'
468 0 : else if (timestep_model_outputs .eq. -1) then
469 0 : fluxesUnit = 'mm d-1'
470 0 : else if (timestep_model_outputs .eq. -2) then
471 0 : fluxesUnit = 'mm month-1'
472 0 : else if (timestep_model_outputs .eq. -3) then
473 0 : fluxesUnit = 'mm a-1'
474 : else
475 0 : fluxesUnit = ''
476 : end if
477 :
478 15 : end function fluxesUnit
479 :
480 : end module mo_write_fluxes_states
|
