Line data Source code
1 : !> \file mo_mrm_riv_temp_class.f90
2 : !> \brief \copybrief mo_mrm_riv_temp_class
3 : !> \details \copydetails mo_mrm_riv_temp_class
4 :
5 : !> \brief Class for the river temperature calculations
6 : !> \details River temperature routing on top of mRM.
7 : !> \warning This feature is still experimental! River freezing is still missing.
8 : !> \version 0.1
9 : !> \authors Sebastian Mueller
10 : !> \date Sep 2020
11 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
12 : !! mHM is released under the LGPLv3+ license \license_note
13 : !> \ingroup f_mrm
14 : module mo_mrm_riv_temp_class
15 :
16 : use mo_kind, only: dp, i4
17 :
18 : implicit none
19 :
20 : private
21 :
22 : !> \class riv_temp_type
23 : !> \brief This is a container to define the river temperature routing in the current time step
24 : !> \details This class provides all procedures to rout river tmperature through the river network.
25 : !> \warning This feature is still experimental!
26 : !! This first version doesn't provide ice covering, which means, that the river temperature
27 : !! can drop below 0 [deg C] in winter.
28 : type, public :: riv_temp_type
29 : ! config settings
30 : logical :: active = .false. !< state if this process is active
31 : integer(i4) :: case = 0_i4 !< the selected process-case option
32 : character(256) :: nml_name = 'config_riv_temp' !< namelist name in mhm namelist
33 : ! riv geometry
34 : character(256), dimension(:), allocatable :: dir_riv_widths !< Directory where river widths are stored
35 : character(256) :: riv_widths_file !< file name for river widths
36 : character(256) :: riv_widths_name !< variable name for river widths
37 : real(dp), public, dimension(:), allocatable :: L11_riv_widths !< river widths in L11
38 : real(dp), public, dimension(:), allocatable :: L11_riv_areas !< river area in L11
39 : ! PET related vars
40 : real(dp) :: albedo_water !< albedo of open water
41 : real(dp) :: pt_a_water !< priestley taylor alpha parameter for PET on open water
42 : ! sensible heat flux related
43 : real(dp) :: emissivity_water !< emissivity of water
44 : real(dp) :: turb_heat_ex_coeff !< lateral heat exchange coefficient water <-> air
45 : !> cutoff value for temperature
46 : real(dp) :: delta_T = 0.1_dp
47 : ! controlling variables for iterative solver
48 : integer(i4) :: max_iter !< input: maximal number of iterations done
49 : real(dp) :: delta_iter !< input: convergence criteria for iterative solver
50 : real(dp) :: step_iter !< input: step-size for linear search
51 : logical :: first_iter !< whether it is at the first iteration (to determine search direction)
52 : logical :: up_iter !< whether the search direction is upwards
53 : logical :: bisect_iter !< whether to do the bisection search part (after the interval is found)
54 : real(dp) :: up_bnd_iter !< upper bound for the current bisection step
55 : real(dp) :: low_bnd_iter !< lower bound for the current bisection step
56 : ! accumulated later fluxes (in current time-step)
57 : real(dp), dimension(:), allocatable :: L1_runoff_E !< runoff energy at L1 level
58 : real(dp), dimension(:), allocatable :: L1_acc_ssrd !< accumulated shortwave radiation at L1 level
59 : real(dp), dimension(:), allocatable :: L1_acc_strd !< accumulated longwave radiation at L1 level
60 : real(dp), dimension(:), allocatable :: L1_acc_temp !< accumulated air temperature at L1 level
61 : real(dp), dimension(:), allocatable :: L1_ssrd_calc !< current shortwave radiation at L1 level
62 : real(dp), dimension(:), allocatable :: L1_strd_calc !< current longwave radiation at L1 level
63 : real(dp), dimension(:), allocatable :: L1_tann_calc !< current mean air temperature at L1 level
64 : integer(i4) :: ts_cnt !< sub time-step counter for accumulation of meteo
65 : ! vars for routing
66 : integer(i4) :: s11 !< starting index for current L11 domain
67 : integer(i4) :: e11 !< ending index for current L11 domain
68 : real(dp), dimension(:,:), allocatable :: netNode_E_IN !< Total energy inputs at t-1 and t
69 : real(dp), dimension(:,:), allocatable :: netNode_E_R !< energy leaving at t-1 and t
70 : real(dp), dimension(:), allocatable :: netNode_E_mod !< Simulated routed energy
71 : real(dp), dimension(:), allocatable :: netNode_E_out !< total energy source from cell in L11
72 : real(dp), dimension(:), allocatable :: L11_srad_net !< net short wave radiation at L11
73 : real(dp), dimension(:), allocatable :: L11_lrad_in !< incoming long wave radiation at L11
74 : real(dp), dimension(:), allocatable :: L11_air_temp !< air temp at L11
75 : ! variable containing river temp for each domain at L11 level
76 : real(dp), dimension(:), allocatable :: river_temp !< resulting river temp at L11 in [deg C]
77 : contains
78 : ! config and inits
79 : !> \copydoc mo_mrm_riv_temp_class::config
80 : procedure :: config !< \see mo_mrm_riv_temp_class::config
81 : !> \copydoc mo_mrm_riv_temp_class::init
82 : procedure :: init !< \see mo_mrm_riv_temp_class::init
83 : !> \copydoc mo_mrm_riv_temp_class::init_area
84 : procedure :: init_area !< \see mo_mrm_riv_temp_class::init_area
85 : !> \copydoc mo_mrm_riv_temp_class::init_riv_temp
86 : procedure :: init_riv_temp !< \see mo_mrm_riv_temp_class::init_riv_temp
87 :
88 : ! source accumulations
89 : !> \copydoc mo_mrm_riv_temp_class::acc_source_e
90 : procedure :: acc_source_E !< \see mo_mrm_riv_temp_class::acc_source_e
91 : !> \copydoc mo_mrm_riv_temp_class::finalize_source_e
92 : procedure :: finalize_source_E !< \see mo_mrm_riv_temp_class::finalize_source_e
93 :
94 : ! temp-energy routing routines
95 : !> \copydoc mo_mrm_riv_temp_class::get_lrad_out
96 : procedure :: get_lrad_out !< \see mo_mrm_riv_temp_class::get_lrad_out
97 : !> \copydoc mo_mrm_riv_temp_class::get_lat_heat
98 : procedure :: get_lat_heat !< \see mo_mrm_riv_temp_class::get_lat_heat
99 : !> \copydoc mo_mrm_riv_temp_class::get_sens_heat
100 : procedure :: get_sens_heat !< \see mo_mrm_riv_temp_class::get_sens_heat
101 : !> \copydoc mo_mrm_riv_temp_class::get_e_io
102 : procedure :: get_E_IO !< \see mo_mrm_riv_temp_class::get_e_io
103 : !> \copydoc mo_mrm_riv_temp_class::l11_routing_e
104 : procedure :: L11_routing_E !< \see mo_mrm_riv_temp_class::l11_routing_e
105 :
106 : ! helper for iterative solver
107 : !> \copydoc mo_mrm_riv_temp_class::init_iter
108 : procedure :: init_iter !< \see mo_mrm_riv_temp_class::init_iter
109 : !> \copydoc mo_mrm_riv_temp_class::next_iter
110 : procedure :: next_iter !< \see mo_mrm_riv_temp_class::next_iter
111 :
112 : ! care taker
113 : !> \copydoc mo_mrm_riv_temp_class::reset_timestep
114 : procedure :: reset_timestep !< \see mo_mrm_riv_temp_class::reset_timestep
115 : !> \copydoc mo_mrm_riv_temp_class::alloc_lateral
116 : procedure :: alloc_lateral !< \see mo_mrm_riv_temp_class::alloc_lateral
117 : !> \copydoc mo_mrm_riv_temp_class::dealloc_lateral
118 : procedure :: dealloc_lateral !< \see mo_mrm_riv_temp_class::dealloc_lateral
119 : !> \copydoc mo_mrm_riv_temp_class::clean_up
120 : procedure :: clean_up !< \see mo_mrm_riv_temp_class::clean_up
121 :
122 : end type riv_temp_type
123 :
124 : contains
125 :
126 :
127 : !> \brief clean up
128 14 : subroutine clean_up( &
129 : self &
130 : )
131 : implicit none
132 :
133 : class(riv_temp_type), intent(inout) :: self
134 :
135 0 : if ( allocated(self%L1_runoff_E) ) deallocate(self%L1_runoff_E)
136 14 : if ( allocated(self%L1_acc_strd) ) deallocate(self%L1_acc_strd)
137 14 : if ( allocated(self%L1_acc_ssrd) ) deallocate(self%L1_acc_ssrd)
138 14 : if ( allocated(self%L1_acc_temp) ) deallocate(self%L1_acc_temp)
139 14 : if ( allocated(self%dir_riv_widths) ) deallocate(self%dir_riv_widths)
140 14 : if ( allocated(self%L11_riv_widths) ) deallocate(self%L11_riv_widths)
141 14 : if ( allocated(self%L11_riv_areas) ) deallocate(self%L11_riv_areas)
142 14 : if ( allocated(self%netNode_E_IN) ) deallocate(self%netNode_E_IN)
143 14 : if ( allocated(self%netNode_E_R) ) deallocate(self%netNode_E_R)
144 14 : if ( allocated(self%netNode_E_mod) ) deallocate(self%netNode_E_mod)
145 14 : if ( allocated(self%netNode_E_out) ) deallocate(self%netNode_E_out)
146 14 : if ( allocated(self%L11_srad_net) ) deallocate(self%L11_srad_net)
147 14 : if ( allocated(self%L11_lrad_in) ) deallocate(self%L11_lrad_in)
148 14 : if ( allocated(self%L11_air_temp) ) deallocate(self%L11_air_temp)
149 14 : if ( allocated(self%river_temp) ) deallocate(self%river_temp)
150 : ! meteo arrays
151 14 : if ( allocated(self%L1_ssrd_calc) ) deallocate(self%L1_ssrd_calc)
152 14 : if ( allocated(self%L1_strd_calc) ) deallocate(self%L1_strd_calc)
153 14 : if ( allocated(self%L1_tann_calc) ) deallocate(self%L1_tann_calc)
154 :
155 14 : end subroutine clean_up
156 :
157 : !> \brief configure the \ref riv_temp_type class from the mhm namelist
158 1 : subroutine config( &
159 : self, &
160 : file_namelist, &
161 : unamelist, &
162 : file_namelist_param, &
163 : unamelist_param &
164 : )
165 :
166 14 : use mo_common_constants, only : maxNoDomains, nodata_i4
167 : use mo_common_variables, only : domainMeta
168 : use mo_nml, only : close_nml, open_nml, position_nml
169 : use mo_check, only : check_dir
170 : USE mo_string_utils, ONLY : num2str
171 :
172 : implicit none
173 :
174 : class(riv_temp_type), intent(inout) :: self
175 : character(*), intent(in) :: file_namelist !< mhm namelist file
176 : character(*), intent(in) :: file_namelist_param !< mhm parameter namelist file
177 : integer, intent(in) :: unamelist
178 : integer, intent(in) :: unamelist_param
179 :
180 : ! parameter to read in
181 1 : real(dp) :: albedo_water ! albedo of open water
182 1 : real(dp) :: pt_a_water ! priestley taylor alpha parameter for PET on open water
183 1 : real(dp) :: emissivity_water ! emissivity of water
184 1 : real(dp) :: turb_heat_ex_coeff ! lateral heat exchange coefficient water <-> air
185 : ! controlling variables for iterative solver
186 : integer(i4) :: max_iter
187 1 : real(dp) :: delta_iter
188 1 : real(dp) :: step_iter
189 : ! files for river widths
190 : character(256), dimension(maxNoDomains) :: dir_riv_widths
191 : character(256) :: riv_widths_file ! file name for river widths
192 : character(256) :: riv_widths_name ! variable name for river widths
193 :
194 : integer(i4) :: iDomain, domainID
195 :
196 : ! namelist for river temperature configuration
197 : namelist /config_riv_temp/ &
198 : albedo_water, &
199 : pt_a_water, &
200 : emissivity_water, &
201 : turb_heat_ex_coeff, &
202 : max_iter, &
203 : delta_iter, &
204 : step_iter, &
205 : riv_widths_file, &
206 : riv_widths_name, &
207 : dir_riv_widths
208 :
209 : ! allocate the directory arrays
210 3 : allocate(self%dir_riv_widths(domainMeta%nDomains))
211 :
212 : ! default values
213 1 : albedo_water = 0.15_dp
214 1 : pt_a_water = 1.26_dp
215 1 : emissivity_water = 0.96_dp
216 1 : turb_heat_ex_coeff = 20.0_dp
217 1 : max_iter = 20_i4
218 1 : delta_iter = 1.0e-02_dp
219 1 : step_iter = 5.0_dp
220 :
221 : ! open the namelist file
222 1 : call open_nml(file_namelist, unamelist, quiet=.true.)
223 : ! find the river-temp config namelist
224 1 : call position_nml(self%nml_name, unamelist)
225 : ! read the river-temp config namelist
226 1 : read(unamelist, nml=config_riv_temp)
227 :
228 1 : self%albedo_water = albedo_water
229 1 : self%pt_a_water = pt_a_water
230 1 : self%emissivity_water = emissivity_water
231 1 : self%turb_heat_ex_coeff = turb_heat_ex_coeff
232 1 : self%delta_iter = delta_iter
233 1 : self%max_iter = max_iter
234 1 : self%step_iter = step_iter
235 1 : self%riv_widths_file = riv_widths_file
236 1 : self%riv_widths_name = riv_widths_name
237 52 : self%dir_riv_widths = dir_riv_widths
238 :
239 : ! closing the namelist file
240 1 : call close_nml(unamelist)
241 :
242 2 : do iDomain = 1, domainMeta%nDomains
243 1 : domainID = domainMeta%indices(iDomain)
244 : call check_dir( &
245 0 : path=self%dir_riv_widths(iDomain), &
246 : text="(domain "//trim(num2str(domainID,'(I3)'))//") River widths directory:", &
247 : raise=.true., &
248 : tab=4, &
249 : text_length=40 &
250 2 : )
251 : end do
252 :
253 1 : end subroutine config
254 :
255 : !> \brief initalize the \ref riv_temp_type class for the current domain
256 1 : subroutine init( &
257 : self, &
258 : nCells &
259 : )
260 1 : use mo_append, only : append
261 : use mo_mrm_constants, only : nRoutingStates
262 : ! use mo_common_variables, only : level0, domainMeta
263 :
264 : implicit none
265 :
266 : class(riv_temp_type), intent(inout) :: self
267 : integer(i4), intent(in) :: nCells !< number of level-11 cells for the current domain
268 :
269 1 : real(dp), dimension(:), allocatable :: dummy_Vector11
270 1 : real(dp), dimension(:, :), allocatable :: dummy_Matrix11_IT
271 :
272 : ! dummy vector and matrix
273 3 : allocate(dummy_Vector11 (nCells))
274 3 : allocate(dummy_Matrix11_IT(nCells, nRoutingStates))
275 35 : dummy_Vector11(:) = 0.0_dp
276 71 : dummy_Matrix11_IT(:, :) = 0.0_dp
277 :
278 : ! simulated energy flux at each node
279 1 : call append(self%netNode_E_mod, dummy_Vector11)
280 : ! simulated river temperature at each node
281 1 : call append(self%river_temp, dummy_Vector11)
282 : ! Total outflow from cells L11 at time tt
283 1 : call append(self%netNode_E_out, dummy_Vector11)
284 : ! Total discharge inputs at t-1 and t
285 1 : call append(self%netNode_E_IN, dummy_Matrix11_IT)
286 : ! Routed outflow leaving a node
287 1 : call append(self%netNode_E_R, dummy_Matrix11_IT)
288 :
289 : ! free space
290 1 : deallocate(dummy_Vector11)
291 1 : deallocate(dummy_Matrix11_IT)
292 :
293 1 : end subroutine init
294 :
295 : !> \brief initialize the river area of \ref riv_temp_type class for the current domain
296 1 : subroutine init_area( &
297 : self, &
298 : iDomain, &
299 2 : L11_netPerm, &
300 1 : L11_fromN, &
301 1 : L11_length, &
302 : nLinks, &
303 : nCells, &
304 : nrows, &
305 : ncols, &
306 1 : L11_mask &
307 : )
308 1 : use mo_append, only : append
309 : use mo_read_nc, only: read_const_nc
310 :
311 : implicit none
312 :
313 : class(riv_temp_type), intent(inout) :: self
314 : integer(i4), intent(in) :: iDomain !< Domain ID
315 : integer(i4), dimension(:), intent(in) :: L11_netPerm !< L11 routing order
316 : integer(i4), dimension(:), intent(in) :: L11_fromN !< L11 source grid cell order
317 : real(dp), dimension(:), intent(in) :: L11_length !< L11 link length
318 : integer(i4), intent(in) :: nLinks !< number of L11 links in the current domain
319 : integer(i4), intent(in) :: nCells !< number of L11 cells of the current domain
320 : integer(i4), intent(in) :: ncols !< Number of columns
321 : integer(i4), intent(in) :: nrows !< Number of rows
322 : logical, dimension(:, :), intent(in) :: L11_mask !< the mask for valid cells in the original grid (nrows, ncols)
323 :
324 1 : real(dp), dimension(:,:), allocatable :: L11_data ! read data from file
325 1 : real(dp), dimension(:), allocatable :: L11_riv_widths, L11_riv_areas
326 :
327 : integer(i4) :: i, k, iNode
328 :
329 : call read_const_nc(&
330 0 : trim(self%dir_riv_widths(iDomain)), &
331 : nrows, &
332 : ncols, &
333 : self%riv_widths_name, &
334 : L11_data, &
335 : self%riv_widths_file &
336 1 : )
337 :
338 3 : allocate(L11_riv_widths(nCells))
339 1 : L11_riv_widths(:) = pack(L11_data(:,:), mask=L11_mask)
340 1 : call append(self%L11_riv_widths, L11_riv_widths)
341 :
342 3 : allocate(L11_riv_areas(nCells))
343 : ! at area at Outlet is 0
344 36 : L11_riv_areas = 0.0_dp
345 34 : do k = 1, nLinks
346 : ! get LINK routing order -> i
347 33 : i = L11_netPerm(k)
348 33 : iNode = L11_fromN(i)
349 34 : L11_riv_areas(iNode) = L11_length(i) * L11_riv_widths(iNode)
350 : end do
351 1 : call append(self%L11_riv_areas, L11_riv_areas)
352 :
353 1 : deallocate(L11_data)
354 1 : deallocate(L11_riv_widths)
355 1 : deallocate(L11_riv_areas)
356 :
357 1 : end subroutine init_area
358 :
359 : !> \brief initialize the river temperature of \ref riv_temp_type class for the current domain
360 1 : subroutine init_riv_temp( &
361 : self, &
362 2 : temp_air, &
363 1 : efecarea, &
364 1 : L1_L11_Id, &
365 1 : L11_areacell, &
366 1 : L11_L1_Id, &
367 : map_flag &
368 : )
369 :
370 1 : use mo_mrm_pre_routing, only : L11_meteo_acc
371 :
372 : implicit none
373 :
374 : class(riv_temp_type), intent(inout) :: self
375 : real(dp), dimension(:), intent(in) :: temp_air !< air temperature [degC] for current timestep
376 : real(dp), intent(in), dimension(:) :: efecarea !< effective area in [km2] at Level 1
377 : integer(i4), intent(in), dimension(:) :: L1_L11_Id !< L11 Ids mapped on L1
378 : real(dp), intent(in), dimension(:) :: L11_areacell !< effective area in [km2] at Level 11
379 : integer(i4), intent(in), dimension(:) :: L11_L1_Id !< L1 Ids mapped on L11
380 : logical, intent(in) :: map_flag !< Flag indicating whether routing resolution is higher than hydrologic one
381 :
382 : ! map temperature from L1 to L11
383 : call L11_meteo_acc( &
384 1 : temp_air, efecarea, L1_L11_Id, L11_areacell, L11_L1_Id, map_flag, self%river_temp(self%s11 : self%e11))
385 :
386 : ! assure positive temperature
387 35 : self%river_temp(self%s11 : self%e11) = max(self%delta_T, self%river_temp(self%s11 : self%e11))
388 :
389 1 : end subroutine init_riv_temp
390 :
391 : !> \brief reset \ref riv_temp_type class for next timestep
392 6553 : subroutine reset_timestep(self)
393 : implicit none
394 :
395 : class(riv_temp_type), intent(inout) :: self
396 :
397 : ! set sub time-step counter to 0
398 6553 : self%ts_cnt = 0_i4
399 :
400 229355 : self%L1_runoff_E = 0.0_dp
401 235908 : self%L1_acc_strd = 0.0_dp
402 235908 : self%L1_acc_ssrd = 0.0_dp
403 235908 : self%L1_acc_temp = 0.0_dp
404 :
405 1 : end subroutine reset_timestep
406 :
407 : !> \brief allocate lateral temp components of \ref riv_temp_type class for current domain
408 1 : subroutine alloc_lateral( &
409 : self, &
410 : nCells &
411 : )
412 : implicit none
413 :
414 : class(riv_temp_type), intent(inout) :: self
415 : integer(i4), intent(in) :: nCells !< number of level-1 cells for the current domain
416 :
417 3 : allocate(self%L1_runoff_E(nCells))
418 3 : allocate(self%L1_acc_strd(nCells))
419 3 : allocate(self%L1_acc_ssrd(nCells))
420 3 : allocate(self%L1_acc_temp(nCells))
421 : ! meteo arrays
422 3 : allocate(self%L1_ssrd_calc(nCells))
423 3 : allocate(self%L1_strd_calc(nCells))
424 3 : allocate(self%L1_tann_calc(nCells))
425 : ! init these arrays to 0
426 1 : call self%reset_timestep()
427 :
428 6553 : end subroutine alloc_lateral
429 :
430 : !> \brief deallocate lateral temp components of \ref riv_temp_type
431 1 : subroutine dealloc_lateral( &
432 : self &
433 : )
434 : implicit none
435 :
436 : class(riv_temp_type), intent(inout) :: self
437 :
438 0 : deallocate(self%L1_runoff_E)
439 1 : deallocate(self%L1_acc_strd)
440 1 : deallocate(self%L1_acc_ssrd)
441 1 : deallocate(self%L1_acc_temp)
442 1 : deallocate(self%L1_ssrd_calc)
443 1 : deallocate(self%L1_strd_calc)
444 1 : deallocate(self%L1_tann_calc)
445 :
446 1 : end subroutine dealloc_lateral
447 :
448 : !> \brief accumulate energy sources of \ref riv_temp_type
449 13104 : subroutine acc_source_e( &
450 : self, &
451 26208 : fSealed_area_fraction, &
452 13104 : fast_interflow, &
453 13104 : slow_interflow, &
454 13104 : baseflow, &
455 13104 : direct_runoff, &
456 13104 : temp_air &
457 : )
458 :
459 1 : use mo_mrm_pre_routing, only : calc_L1_runoff_E
460 :
461 : implicit none
462 :
463 : class(riv_temp_type), intent(inout) :: self
464 : real(dp), dimension(:), intent(in) :: fSealed_area_fraction !< sealed area fraction [1]
465 : real(dp), dimension(:), intent(in) :: fast_interflow !< \f$ q_0 \f$ Fast runoff component [mm TS-1]
466 : real(dp), dimension(:), intent(in) :: slow_interflow !< \f$ q_1 \f$ Slow runoff component [mm TS-1]
467 : real(dp), dimension(:), intent(in) :: baseflow !< \f$ q_2 \f$ Baseflow [mm TS-1]
468 : real(dp), dimension(:), intent(in) :: direct_runoff !< \f$ q_D \f$ Direct runoff from impervious areas [mm TS-1]
469 : real(dp), dimension(:), intent(in) :: temp_air !< air temperature [K]
470 :
471 : ! increase the sub time-step counter
472 13104 : self%ts_cnt = self%ts_cnt + 1_i4
473 :
474 : ! caclucate the temperature energy of the runoffs at L1 in [K mm]
475 : ! automatically accumulate them
476 : call calc_L1_runoff_E( &
477 : fSealed_area_fraction, &
478 : fast_interflow, slow_interflow, baseflow, direct_runoff, &
479 0 : temp_air, self%L1_tann_calc, &
480 0 : self%L1_runoff_E & ! will be added here
481 13104 : )
482 : ! accumulate meteo forcings (will be averaged with sub time-step counter later)
483 471744 : self%L1_acc_ssrd = self%L1_acc_ssrd + self%L1_ssrd_calc
484 471744 : self%L1_acc_strd = self%L1_acc_strd + self%L1_strd_calc
485 471744 : self%L1_acc_temp = self%L1_acc_temp + temp_air
486 :
487 13104 : end subroutine acc_source_e
488 :
489 : !> \brief finalize energy sources of \ref riv_temp_type
490 6552 : subroutine finalize_source_E( &
491 : self, &
492 13104 : efecarea, &
493 6552 : L1_L11_Id, &
494 6552 : L11_areacell, &
495 6552 : L11_L1_Id, &
496 : timestep, &
497 : map_flag &
498 : )
499 :
500 13104 : use mo_mrm_pre_routing, only : L11_meteo_acc, L11_runoff_acc
501 :
502 : implicit none
503 :
504 : class(riv_temp_type), intent(inout) :: self
505 : !> effective area in [km2] at Level 1
506 : real(dp), intent(in), dimension(:) :: efecarea
507 : !> L11 Ids mapped on L1
508 : integer(i4), intent(in), dimension(:) :: L1_L11_Id
509 : !> effective area in [km2] at Level 11
510 : real(dp), intent(in), dimension(:) :: L11_areacell
511 : !> L1 Ids mapped on L11
512 : integer(i4), intent(in), dimension(:) :: L11_L1_Id
513 : !> simulation timestep in [h]
514 : integer(i4), intent(in) :: timestep
515 : !> Flag indicating whether routing resolution is higher than hydrologic one
516 : logical, intent(in) :: map_flag
517 :
518 : ! prepare temporal variables for temp-routing
519 6552 : if ( allocated(self%L11_srad_net) ) deallocate(self%L11_srad_net)
520 6552 : if ( allocated(self%L11_lrad_in) ) deallocate(self%L11_lrad_in)
521 6552 : if ( allocated(self%L11_air_temp) ) deallocate(self%L11_air_temp)
522 19656 : allocate(self%L11_srad_net(size(L11_L1_Id, dim = 1)))
523 13104 : allocate(self%L11_lrad_in(size(L11_L1_Id, dim = 1)))
524 13104 : allocate(self%L11_air_temp(size(L11_L1_Id, dim = 1)))
525 :
526 : ! convert the L1 runoff temperature energy from [K mm] to [K m3 s-1] on L11
527 : call L11_runoff_acc( &
528 0 : self%L1_runoff_E, &
529 : efecarea, &
530 : L1_L11_Id, &
531 : L11_areacell, &
532 : L11_L1_Id, &
533 : timestep, &
534 : map_flag, &
535 0 : self%netNode_E_out(self%s11 : self%e11) &
536 6552 : )
537 : ! all meteo forcings are accumulated over sub time-steps and need to be divided by the sub time-step counter
538 : ! short wave radiation on L11
539 229320 : self%L1_acc_ssrd = self%L1_acc_ssrd / real(self%ts_cnt, dp)
540 6552 : call L11_meteo_acc(self%L1_acc_ssrd, efecarea, L1_L11_Id, L11_areacell, L11_L1_Id, map_flag, self%L11_srad_net)
541 : ! calculate the net-radiation from incoming/outging short radiation
542 229320 : self%L11_srad_net = self%L11_srad_net * (1._dp - self%albedo_water)
543 : ! long wave radiation on L11
544 229320 : self%L1_acc_strd = self%L1_acc_strd / real(self%ts_cnt, dp)
545 6552 : call L11_meteo_acc(self%L1_acc_strd, efecarea, L1_L11_Id, L11_areacell, L11_L1_Id, map_flag, self%L11_lrad_in)
546 : ! air temperature on L11
547 229320 : self%L1_acc_temp = self%L1_acc_temp / real(self%ts_cnt, dp)
548 6552 : call L11_meteo_acc(self%L1_acc_temp, efecarea, L1_L11_Id, L11_areacell, L11_L1_Id, map_flag, self%L11_air_temp)
549 :
550 6552 : end subroutine finalize_source_E
551 :
552 : !> \brief get outgoing longwave radiation of \ref riv_temp_type
553 : !> \return outgoing longwave radiation
554 5550108 : real(dp) function get_lrad_out(self, riv_temp) result(lrad_out)
555 :
556 6552 : use mo_constants, only: sigma_dp
557 :
558 : implicit none
559 :
560 : class(riv_temp_type), intent(in) :: self
561 : !> river temperature in K
562 : real(dp), intent(in) :: riv_temp
563 :
564 : ! outgoing longwave radiation from Boltzmann equation
565 5550108 : lrad_out = self%emissivity_water * sigma_dp * riv_temp ** 4_i4
566 :
567 5550108 : end function get_lrad_out
568 :
569 : !> \brief latent heat flux of \ref riv_temp_type
570 : !> \return latent heat flux
571 5550108 : real(dp) function get_lat_heat(self, air_temp, netrad) result(lat_heat)
572 :
573 5550108 : use mo_constants, only : Psychro_dp
574 : use mo_pet, only : slope_satpressure
575 :
576 : implicit none
577 :
578 : class(riv_temp_type), intent(in) :: self
579 : !> air temperature in deg C
580 : real(dp), intent(in) :: air_temp
581 : !> net radiation in W * m-2
582 : real(dp), intent(in) :: netrad
583 :
584 : ! save slope of saturation vapor pressure curve
585 5550108 : real(dp) :: delta
586 :
587 5550108 : delta = slope_satpressure(air_temp) ! slope of saturation vapor pressure curve
588 : ! latent heat for priestley taylor PET formula
589 5550108 : lat_heat = self%pt_a_water * delta / (Psychro_dp + delta) * max(netrad, 0._dp)
590 :
591 5550108 : end function get_lat_heat
592 :
593 : !> \brief sensible heat flux of \ref riv_temp_type
594 : !> \return sensible heat flux
595 5550108 : real(dp) function get_sens_heat(self, air_temp, riv_temp) result(sens_heat)
596 :
597 : implicit none
598 :
599 : class(riv_temp_type), intent(in) :: self
600 : real(dp), intent(in) :: air_temp !< air temperature in [deg C]
601 : real(dp), intent(in) :: riv_temp !< river temperature in [deg C]
602 :
603 : ! sensible heat flux resulting for temp. diff.: river <-> air
604 5550108 : sens_heat = self%turb_heat_ex_coeff * (riv_temp - air_temp)
605 :
606 11100216 : end function get_sens_heat
607 :
608 : !> \brief get complete energy source of \ref riv_temp_type at given cell
609 : !> \return energy IO
610 5550108 : real(dp) function get_E_IO(self, riv_temp, cell) result(E_IO)
611 :
612 5550108 : use mo_constants, only : T0_dp, cp_w_dp
613 : use mo_mhm_constants, only : H2Odens
614 :
615 : implicit none
616 :
617 : class(riv_temp_type), intent(in) :: self
618 : !> given river temperature in K to calculate heat fluxes
619 : real(dp), intent(in) :: riv_temp
620 : !> cell index in the current domain
621 : integer(i4), intent(in) :: cell
622 : !> net radiation calc from short and longwave radiation in/out
623 5550108 : real(dp) :: netrad, sens_heat, lat_heat
624 :
625 : ! net radiation
626 5550108 : netrad = self%L11_srad_net(cell) + self%L11_lrad_in(cell) - self%get_lrad_out(riv_temp)
627 : ! sensible heat flux
628 5550108 : sens_heat = self%get_sens_heat(self%L11_air_temp(cell), riv_temp - T0_dp)
629 : ! latent heat flux
630 5550108 : lat_heat = self%get_lat_heat(self%L11_air_temp(cell), netrad)
631 :
632 : ! convert energy flux [W m-2] to [K m s-1]
633 : ! needs to be multiplied with river-area to result in temp-energy flux [K m3 s-1]
634 5550108 : E_IO = (netrad - sens_heat - lat_heat) / (H2Odens * cp_w_dp)
635 :
636 5550108 : end function get_E_IO
637 :
638 : !> \brief execute the temperature routing of \ref riv_temp_type
639 6552 : subroutine L11_routing_E( &
640 : self, &
641 : nLinks, &
642 13104 : netPerm, &
643 6552 : netLink_fromN, &
644 6552 : netLink_toN, &
645 13104 : netLink_C1, &
646 6552 : netLink_C2, &
647 : nInflowGauges, &
648 13104 : InflowHeadwater, &
649 6552 : InflowNodeList, &
650 6552 : L11_qTR, &
651 6552 : L11_Qmod &
652 : )
653 5550108 : use mo_constants, only : T0_dp
654 :
655 : implicit none
656 :
657 : class(riv_temp_type), intent(inout) :: self
658 : !> number of stream segment (reaches)
659 : integer(i4), intent(in) :: nLinks
660 : !> routing order of a given domain (permutation)
661 : integer(i4), dimension(:), intent(in) :: netPerm
662 : !> from node
663 : integer(i4), dimension(:), intent(in) :: netLink_fromN
664 : !> to node
665 : integer(i4), dimension(:), intent(in) :: netLink_toN
666 : !> routing parameter C1 (\cite CMM1988 p. 25-41)
667 : real(dp), dimension(:), intent(in) :: netLink_C1
668 : !> routing parameters C2 (id)
669 : real(dp), dimension(:), intent(in) :: netLink_C2
670 : !> [-] number of inflow points
671 : integer(i4), intent(in) :: nInflowGauges
672 : !> [-] if to consider headwater cells of inflow gauge
673 : logical, dimension(:), intent(in) :: InflowHeadwater
674 : !> [-] L11 ID of inflow points
675 : integer(i4), dimension(:), intent(in) :: InflowNodeList
676 : !> [m3 s-1] Transformed outflow leaving node I at current timestep(Muskingum)
677 : real(dp), intent(in), dimension(:) :: L11_qTR
678 : !> [m3 s-1] Simulated routed discharge
679 : real(dp), intent(in), dimension(:) :: L11_Qmod
680 :
681 : integer(i4) :: i, k, m, iNode, tNode, L11in, L11to
682 6552 : real(dp) :: E_IO_in, E_IO, T_est, T_rout
683 :
684 : ! initialize total input at point time (t+1) in all nodes
685 229320 : self%netNode_E_IN(self%s11 : self%e11, 2) = 0.0_dp
686 :
687 222768 : riv_loop: do k = 1, nLinks
688 : ! get LINK routing order -> i
689 216216 : i = netPerm(k)
690 216216 : iNode = netLink_fromN(i)
691 216216 : tNode = netLink_toN(i)
692 : ! indices for concaternated domains
693 216216 : L11in = iNode + self%s11 - 1_i4
694 216216 : L11to = tNode + self%s11 - 1_i4
695 :
696 : ! accumulate all inputs in iNode
697 216216 : self%netNode_E_IN(L11in, 2) = self%netNode_E_IN(L11in, 2) + self%netNode_E_out(L11in)
698 : ! calculate the river temp at the upstream node
699 216216 : self%river_temp(L11in) = self%netNode_E_IN(L11in, 2) / L11_Qmod(iNode) - T0_dp
700 :
701 : ! first guess for the routed temperature is the upstream temp
702 216216 : T_est = self%river_temp(L11in) + T0_dp
703 : ! calculate the meteo energy source at the upstream node
704 216216 : E_IO_in = self%get_E_IO(T_est, iNode)
705 :
706 : ! initialize iteration
707 216216 : call self%init_iter()
708 : ! perform iteration
709 5333892 : iterloop: do m=1, self%max_iter
710 : ! meteo energy flux depending on the estimated routed temperature (avg. IN- and OUT-node)
711 5333892 : E_IO = (E_IO_in + self%get_E_IO(T_est, tNode)) * 0.5_dp * self%L11_riv_areas(L11in)
712 : ! routing iNode
713 5333892 : self%netNode_E_R(L11in, 2) = self%netNode_E_R(L11in, 1) &
714 5333892 : + netLink_C1(i) * (self%netNode_E_IN(L11in, 1) - self%netNode_E_R(L11in, 1)) &
715 10667784 : + netLink_C2(i) * (self%netNode_E_IN(L11in, 2) + E_IO - self%netNode_E_IN(L11in, 1))
716 : ! calculate the routed temperature
717 5333892 : T_rout = self%netNode_E_R(L11in, 2) / L11_qTR(iNode)
718 : ! check for convergence
719 5333892 : if ( abs(T_rout - T_est) < self%delta_iter ) exit iterloop
720 : ! get new estimate for next iteration
721 5333892 : call self%next_iter(T_est, T_rout)
722 : end do iterloop
723 :
724 : ! add the meteo energy flux to the accumulated incoming energy at the IN-node
725 216216 : self%netNode_E_IN(L11in, 2) = self%netNode_E_IN(L11in, 2) + E_IO
726 :
727 : ! check if the inflow from upstream cells should be deactivated
728 216216 : if (nInflowGauges .GT. 0) then
729 0 : do i = 1, nInflowGauges
730 : ! check if downstream Node (tNode) is inflow gauge and headwaters should be ignored
731 0 : if ((tNode == InflowNodeList(i)) .AND. (.NOT. InflowHeadwater(i))) &
732 0 : self%netNode_E_R(L11in, 2) = 0.0_dp
733 : end do
734 : end if
735 :
736 : ! add routed temp energy to downstream node
737 222768 : self%netNode_E_IN(L11to, 2) = self%netNode_E_IN(L11to, 2) + self%netNode_E_R(L11in, 2)
738 :
739 : end do riv_loop
740 :
741 : ! Accumulate all inputs in tNode (netNode_E_out) ONLY for last link
742 6552 : tNode = netLink_toN(netPerm(nLinks))
743 6552 : L11to = tNode + self%s11 - 1_i4
744 6552 : self%netNode_E_IN(L11to, 2) = self%netNode_E_IN(L11to, 2) + self%netNode_E_out(L11to)
745 : ! calculate riv temp at last link
746 6552 : self%river_temp(L11to) = self%netNode_E_IN(L11to, 2) / L11_Qmod(tNode) - T0_dp
747 :
748 : ! backflow t-> t-1
749 229320 : self%netNode_E_R(self%s11 : self%e11, 1) = self%netNode_E_R(self%s11 : self%e11, 2)
750 229320 : self%netNode_E_IN(self%s11 : self%e11, 1) = self%netNode_E_IN(self%s11 : self%e11, 2)
751 :
752 6552 : end subroutine L11_routing_E
753 :
754 : !> \brief initialize iterative solver of \ref riv_temp_type
755 216216 : subroutine init_iter(self)
756 :
757 : implicit none
758 :
759 : class(riv_temp_type), intent(inout) :: self
760 :
761 : ! first iteration is used to determine the search direction for the interval
762 216216 : self%first_iter = .true.
763 : ! after the interval is found, we will perform a biseciton to nail down the temperature
764 216216 : self%bisect_iter = .false.
765 :
766 6552 : end subroutine init_iter
767 :
768 : !> \brief execute next iteration with iterative solver of \ref riv_temp_type
769 5117676 : subroutine next_iter(self, T_est, T_rout)
770 :
771 : implicit none
772 :
773 : class(riv_temp_type), intent(inout) :: self
774 : real(dp), intent(inout) :: T_est !< estimated river temperature
775 : real(dp), intent(in) :: T_rout !< calculated (routed) river temperature
776 :
777 : ! before performing a bisection we need to search for the interval (with given step-size)
778 5117676 : if ( .not. self%bisect_iter ) then
779 : ! to determine the direction of interval search, we have to wait for the first iteration
780 641761 : if ( self%first_iter ) then
781 216216 : self%first_iter = .false.
782 : ! determine search direction
783 216216 : self%up_iter = ( T_rout > T_est )
784 : end if
785 : ! check if we need to take another step
786 641761 : if ( self%up_iter .eqv. ( T_rout > T_est ) ) then
787 425545 : if ( self%up_iter ) then
788 203382 : T_est = T_est + self%step_iter
789 : else
790 222163 : T_est = T_est - self%step_iter
791 : end if
792 : ! if direction changes, we have found the interval for bisection
793 : else
794 : ! start bisection in next iteration
795 216216 : self%bisect_iter = .true.
796 : ! set interval bounds for bisection
797 216216 : if ( self%up_iter ) then
798 103220 : self%up_bnd_iter = T_est
799 103220 : self%low_bnd_iter = T_est - self%step_iter
800 : else
801 112996 : self%up_bnd_iter = T_est + self%step_iter
802 112996 : self%low_bnd_iter = T_est
803 : end if
804 : ! set estimation to interval center
805 216216 : T_est = (self%up_bnd_iter + self%low_bnd_iter) / 2.0_dp
806 : end if
807 : ! perform bisection
808 : else
809 : ! if routed temp is lower than estimated, use lower interval
810 4475915 : if ( T_rout < T_est ) then
811 2233019 : self%up_bnd_iter = T_est
812 : ! if routed temp is higher than estimated, use upper interval
813 : else
814 2242896 : self%low_bnd_iter = T_est
815 : end if
816 : ! set estimation to interval center
817 4475915 : T_est = (self%up_bnd_iter + self%low_bnd_iter) / 2.0_dp
818 : end if
819 :
820 5333892 : end subroutine next_iter
821 :
822 5117676 : end module mo_mrm_riv_temp_class
|
