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 : file_namelist_param &
162 : )
163 :
164 14 : use mo_common_constants, only : maxNoDomains, nodata_i4
165 : use mo_common_variables, only : domainMeta
166 : use mo_check, only : check_dir
167 : USE mo_string_utils, ONLY : num2str
168 : use mo_namelists, only : nml_config_riv_temp
169 :
170 : implicit none
171 :
172 : class(riv_temp_type), intent(inout) :: self
173 : character(*), intent(in) :: file_namelist !< mhm namelist file
174 : character(*), intent(in) :: file_namelist_param !< mhm parameter namelist file
175 :
176 : integer(i4) :: iDomain, domainID
177 :
178 : ! allocate the directory arrays
179 3 : allocate(self%dir_riv_widths(domainMeta%nDomains))
180 :
181 1 : call nml_config_riv_temp%read(file_namelist)
182 1 : self%albedo_water = nml_config_riv_temp%albedo_water
183 1 : self%pt_a_water = nml_config_riv_temp%pt_a_water
184 1 : self%emissivity_water = nml_config_riv_temp%emissivity_water
185 1 : self%turb_heat_ex_coeff = nml_config_riv_temp%turb_heat_ex_coeff
186 1 : self%delta_iter = nml_config_riv_temp%delta_iter
187 1 : self%max_iter = nml_config_riv_temp%max_iter
188 1 : self%step_iter = nml_config_riv_temp%step_iter
189 1 : self%riv_widths_file = nml_config_riv_temp%riv_widths_file
190 1 : self%riv_widths_name = nml_config_riv_temp%riv_widths_name
191 3 : self%dir_riv_widths = nml_config_riv_temp%dir_riv_widths(1:domainMeta%nDomains)
192 :
193 2 : do iDomain = 1, domainMeta%nDomains
194 1 : domainID = domainMeta%indices(iDomain)
195 : call check_dir( &
196 0 : path=self%dir_riv_widths(iDomain), &
197 : text="(domain "//trim(num2str(domainID,'(I3)'))//") River widths directory:", &
198 : raise=.true., &
199 : tab=4, &
200 : text_length=40 &
201 2 : )
202 : end do
203 :
204 1 : end subroutine config
205 :
206 : !> \brief initalize the \ref riv_temp_type class for the current domain
207 1 : subroutine init( &
208 : self, &
209 : nCells &
210 : )
211 1 : use mo_append, only : append
212 : use mo_mrm_constants, only : nRoutingStates
213 : ! use mo_common_variables, only : level0, domainMeta
214 :
215 : implicit none
216 :
217 : class(riv_temp_type), intent(inout) :: self
218 : integer(i4), intent(in) :: nCells !< number of level-11 cells for the current domain
219 :
220 1 : real(dp), dimension(:), allocatable :: dummy_Vector11
221 1 : real(dp), dimension(:, :), allocatable :: dummy_Matrix11_IT
222 :
223 : ! dummy vector and matrix
224 3 : allocate(dummy_Vector11 (nCells))
225 3 : allocate(dummy_Matrix11_IT(nCells, nRoutingStates))
226 35 : dummy_Vector11(:) = 0.0_dp
227 71 : dummy_Matrix11_IT(:, :) = 0.0_dp
228 :
229 : ! simulated energy flux at each node
230 1 : call append(self%netNode_E_mod, dummy_Vector11)
231 : ! simulated river temperature at each node
232 1 : call append(self%river_temp, dummy_Vector11)
233 : ! Total outflow from cells L11 at time tt
234 1 : call append(self%netNode_E_out, dummy_Vector11)
235 : ! Total discharge inputs at t-1 and t
236 1 : call append(self%netNode_E_IN, dummy_Matrix11_IT)
237 : ! Routed outflow leaving a node
238 1 : call append(self%netNode_E_R, dummy_Matrix11_IT)
239 :
240 : ! free space
241 1 : deallocate(dummy_Vector11)
242 1 : deallocate(dummy_Matrix11_IT)
243 :
244 1 : end subroutine init
245 :
246 : !> \brief initialize the river area of \ref riv_temp_type class for the current domain
247 1 : subroutine init_area( &
248 : self, &
249 : iDomain, &
250 2 : L11_netPerm, &
251 1 : L11_fromN, &
252 1 : L11_length, &
253 : nLinks, &
254 : nCells, &
255 : nrows, &
256 : ncols, &
257 1 : L11_mask &
258 : )
259 1 : use mo_append, only : append
260 : use mo_read_nc, only: read_const_nc
261 :
262 : implicit none
263 :
264 : class(riv_temp_type), intent(inout) :: self
265 : integer(i4), intent(in) :: iDomain !< Domain ID
266 : integer(i4), dimension(:), intent(in) :: L11_netPerm !< L11 routing order
267 : integer(i4), dimension(:), intent(in) :: L11_fromN !< L11 source grid cell order
268 : real(dp), dimension(:), intent(in) :: L11_length !< L11 link length
269 : integer(i4), intent(in) :: nLinks !< number of L11 links in the current domain
270 : integer(i4), intent(in) :: nCells !< number of L11 cells of the current domain
271 : integer(i4), intent(in) :: ncols !< Number of columns
272 : integer(i4), intent(in) :: nrows !< Number of rows
273 : logical, dimension(:, :), intent(in) :: L11_mask !< the mask for valid cells in the original grid (nrows, ncols)
274 :
275 1 : real(dp), dimension(:,:), allocatable :: L11_data ! read data from file
276 1 : real(dp), dimension(:), allocatable :: L11_riv_widths, L11_riv_areas
277 :
278 : integer(i4) :: i, k, iNode
279 :
280 : call read_const_nc(&
281 0 : trim(self%dir_riv_widths(iDomain)), &
282 : nrows, &
283 : ncols, &
284 : self%riv_widths_name, &
285 : L11_data, &
286 : self%riv_widths_file &
287 1 : )
288 :
289 3 : allocate(L11_riv_widths(nCells))
290 1 : L11_riv_widths(:) = pack(L11_data(:,:), mask=L11_mask)
291 1 : call append(self%L11_riv_widths, L11_riv_widths)
292 :
293 3 : allocate(L11_riv_areas(nCells))
294 : ! at area at Outlet is 0
295 36 : L11_riv_areas = 0.0_dp
296 34 : do k = 1, nLinks
297 : ! get LINK routing order -> i
298 33 : i = L11_netPerm(k)
299 33 : iNode = L11_fromN(i)
300 34 : L11_riv_areas(iNode) = L11_length(i) * L11_riv_widths(iNode)
301 : end do
302 1 : call append(self%L11_riv_areas, L11_riv_areas)
303 :
304 1 : deallocate(L11_data)
305 1 : deallocate(L11_riv_widths)
306 1 : deallocate(L11_riv_areas)
307 :
308 1 : end subroutine init_area
309 :
310 : !> \brief initialize the river temperature of \ref riv_temp_type class for the current domain
311 1 : subroutine init_riv_temp( &
312 : self, &
313 2 : temp_air, &
314 1 : efecarea, &
315 1 : L1_L11_Id, &
316 1 : L11_areacell, &
317 1 : L11_L1_Id, &
318 : map_flag &
319 : )
320 :
321 1 : use mo_mrm_pre_routing, only : L11_meteo_acc
322 :
323 : implicit none
324 :
325 : class(riv_temp_type), intent(inout) :: self
326 : real(dp), dimension(:), intent(in) :: temp_air !< air temperature [degC] for current timestep
327 : real(dp), intent(in), dimension(:) :: efecarea !< effective area in [km2] at Level 1
328 : integer(i4), intent(in), dimension(:) :: L1_L11_Id !< L11 Ids mapped on L1
329 : real(dp), intent(in), dimension(:) :: L11_areacell !< effective area in [km2] at Level 11
330 : integer(i4), intent(in), dimension(:) :: L11_L1_Id !< L1 Ids mapped on L11
331 : logical, intent(in) :: map_flag !< Flag indicating whether routing resolution is higher than hydrologic one
332 :
333 : ! map temperature from L1 to L11
334 : call L11_meteo_acc( &
335 1 : temp_air, efecarea, L1_L11_Id, L11_areacell, L11_L1_Id, map_flag, self%river_temp(self%s11 : self%e11))
336 :
337 : ! assure positive temperature
338 35 : self%river_temp(self%s11 : self%e11) = max(self%delta_T, self%river_temp(self%s11 : self%e11))
339 :
340 1 : end subroutine init_riv_temp
341 :
342 : !> \brief reset \ref riv_temp_type class for next timestep
343 6553 : subroutine reset_timestep(self)
344 : implicit none
345 :
346 : class(riv_temp_type), intent(inout) :: self
347 :
348 : ! set sub time-step counter to 0
349 6553 : self%ts_cnt = 0_i4
350 :
351 229355 : self%L1_runoff_E = 0.0_dp
352 235908 : self%L1_acc_strd = 0.0_dp
353 235908 : self%L1_acc_ssrd = 0.0_dp
354 235908 : self%L1_acc_temp = 0.0_dp
355 :
356 1 : end subroutine reset_timestep
357 :
358 : !> \brief allocate lateral temp components of \ref riv_temp_type class for current domain
359 1 : subroutine alloc_lateral( &
360 : self, &
361 : nCells &
362 : )
363 : implicit none
364 :
365 : class(riv_temp_type), intent(inout) :: self
366 : integer(i4), intent(in) :: nCells !< number of level-1 cells for the current domain
367 :
368 3 : allocate(self%L1_runoff_E(nCells))
369 3 : allocate(self%L1_acc_strd(nCells))
370 3 : allocate(self%L1_acc_ssrd(nCells))
371 3 : allocate(self%L1_acc_temp(nCells))
372 : ! meteo arrays
373 3 : allocate(self%L1_ssrd_calc(nCells))
374 3 : allocate(self%L1_strd_calc(nCells))
375 3 : allocate(self%L1_tann_calc(nCells))
376 : ! init these arrays to 0
377 1 : call self%reset_timestep()
378 :
379 6553 : end subroutine alloc_lateral
380 :
381 : !> \brief deallocate lateral temp components of \ref riv_temp_type
382 1 : subroutine dealloc_lateral( &
383 : self &
384 : )
385 : implicit none
386 :
387 : class(riv_temp_type), intent(inout) :: self
388 :
389 0 : deallocate(self%L1_runoff_E)
390 1 : deallocate(self%L1_acc_strd)
391 1 : deallocate(self%L1_acc_ssrd)
392 1 : deallocate(self%L1_acc_temp)
393 1 : deallocate(self%L1_ssrd_calc)
394 1 : deallocate(self%L1_strd_calc)
395 1 : deallocate(self%L1_tann_calc)
396 :
397 1 : end subroutine dealloc_lateral
398 :
399 : !> \brief accumulate energy sources of \ref riv_temp_type
400 13104 : subroutine acc_source_e( &
401 : self, &
402 26208 : fSealed_area_fraction, &
403 13104 : fast_interflow, &
404 13104 : slow_interflow, &
405 13104 : baseflow, &
406 13104 : direct_runoff, &
407 13104 : temp_air &
408 : )
409 :
410 1 : use mo_mrm_pre_routing, only : calc_L1_runoff_E
411 :
412 : implicit none
413 :
414 : class(riv_temp_type), intent(inout) :: self
415 : real(dp), dimension(:), intent(in) :: fSealed_area_fraction !< sealed area fraction [1]
416 : real(dp), dimension(:), intent(in) :: fast_interflow !< \f$ q_0 \f$ Fast runoff component [mm TS-1]
417 : real(dp), dimension(:), intent(in) :: slow_interflow !< \f$ q_1 \f$ Slow runoff component [mm TS-1]
418 : real(dp), dimension(:), intent(in) :: baseflow !< \f$ q_2 \f$ Baseflow [mm TS-1]
419 : real(dp), dimension(:), intent(in) :: direct_runoff !< \f$ q_D \f$ Direct runoff from impervious areas [mm TS-1]
420 : real(dp), dimension(:), intent(in) :: temp_air !< air temperature [K]
421 :
422 : ! increase the sub time-step counter
423 13104 : self%ts_cnt = self%ts_cnt + 1_i4
424 :
425 : ! caclucate the temperature energy of the runoffs at L1 in [K mm]
426 : ! automatically accumulate them
427 : call calc_L1_runoff_E( &
428 : fSealed_area_fraction, &
429 : fast_interflow, slow_interflow, baseflow, direct_runoff, &
430 0 : temp_air, self%L1_tann_calc, &
431 0 : self%L1_runoff_E & ! will be added here
432 13104 : )
433 : ! accumulate meteo forcings (will be averaged with sub time-step counter later)
434 471744 : self%L1_acc_ssrd = self%L1_acc_ssrd + self%L1_ssrd_calc
435 471744 : self%L1_acc_strd = self%L1_acc_strd + self%L1_strd_calc
436 471744 : self%L1_acc_temp = self%L1_acc_temp + temp_air
437 :
438 13104 : end subroutine acc_source_e
439 :
440 : !> \brief finalize energy sources of \ref riv_temp_type
441 6552 : subroutine finalize_source_E( &
442 : self, &
443 13104 : efecarea, &
444 6552 : L1_L11_Id, &
445 6552 : L11_areacell, &
446 6552 : L11_L1_Id, &
447 : timestep, &
448 : map_flag &
449 : )
450 :
451 13104 : use mo_mrm_pre_routing, only : L11_meteo_acc, L11_runoff_acc
452 :
453 : implicit none
454 :
455 : class(riv_temp_type), intent(inout) :: self
456 : !> effective area in [km2] at Level 1
457 : real(dp), intent(in), dimension(:) :: efecarea
458 : !> L11 Ids mapped on L1
459 : integer(i4), intent(in), dimension(:) :: L1_L11_Id
460 : !> effective area in [km2] at Level 11
461 : real(dp), intent(in), dimension(:) :: L11_areacell
462 : !> L1 Ids mapped on L11
463 : integer(i4), intent(in), dimension(:) :: L11_L1_Id
464 : !> simulation timestep in [h]
465 : integer(i4), intent(in) :: timestep
466 : !> Flag indicating whether routing resolution is higher than hydrologic one
467 : logical, intent(in) :: map_flag
468 :
469 : ! prepare temporal variables for temp-routing
470 6552 : if ( allocated(self%L11_srad_net) ) deallocate(self%L11_srad_net)
471 6552 : if ( allocated(self%L11_lrad_in) ) deallocate(self%L11_lrad_in)
472 6552 : if ( allocated(self%L11_air_temp) ) deallocate(self%L11_air_temp)
473 19656 : allocate(self%L11_srad_net(size(L11_L1_Id, dim = 1)))
474 13104 : allocate(self%L11_lrad_in(size(L11_L1_Id, dim = 1)))
475 13104 : allocate(self%L11_air_temp(size(L11_L1_Id, dim = 1)))
476 :
477 : ! convert the L1 runoff temperature energy from [K mm] to [K m3 s-1] on L11
478 : call L11_runoff_acc( &
479 0 : self%L1_runoff_E, &
480 : efecarea, &
481 : L1_L11_Id, &
482 : L11_areacell, &
483 : L11_L1_Id, &
484 : timestep, &
485 : map_flag, &
486 0 : self%netNode_E_out(self%s11 : self%e11) &
487 6552 : )
488 : ! all meteo forcings are accumulated over sub time-steps and need to be divided by the sub time-step counter
489 : ! short wave radiation on L11
490 229320 : self%L1_acc_ssrd = self%L1_acc_ssrd / real(self%ts_cnt, dp)
491 6552 : call L11_meteo_acc(self%L1_acc_ssrd, efecarea, L1_L11_Id, L11_areacell, L11_L1_Id, map_flag, self%L11_srad_net)
492 : ! calculate the net-radiation from incoming/outging short radiation
493 229320 : self%L11_srad_net = self%L11_srad_net * (1._dp - self%albedo_water)
494 : ! long wave radiation on L11
495 229320 : self%L1_acc_strd = self%L1_acc_strd / real(self%ts_cnt, dp)
496 6552 : call L11_meteo_acc(self%L1_acc_strd, efecarea, L1_L11_Id, L11_areacell, L11_L1_Id, map_flag, self%L11_lrad_in)
497 : ! air temperature on L11
498 229320 : self%L1_acc_temp = self%L1_acc_temp / real(self%ts_cnt, dp)
499 6552 : call L11_meteo_acc(self%L1_acc_temp, efecarea, L1_L11_Id, L11_areacell, L11_L1_Id, map_flag, self%L11_air_temp)
500 :
501 6552 : end subroutine finalize_source_E
502 :
503 : !> \brief get outgoing longwave radiation of \ref riv_temp_type
504 : !> \return outgoing longwave radiation
505 5550108 : real(dp) function get_lrad_out(self, riv_temp) result(lrad_out)
506 :
507 6552 : use mo_constants, only: sigma_dp
508 :
509 : implicit none
510 :
511 : class(riv_temp_type), intent(in) :: self
512 : !> river temperature in K
513 : real(dp), intent(in) :: riv_temp
514 :
515 : ! outgoing longwave radiation from Boltzmann equation
516 5550108 : lrad_out = self%emissivity_water * sigma_dp * riv_temp ** 4_i4
517 :
518 5550108 : end function get_lrad_out
519 :
520 : !> \brief latent heat flux of \ref riv_temp_type
521 : !> \return latent heat flux
522 5550108 : real(dp) function get_lat_heat(self, air_temp, netrad) result(lat_heat)
523 :
524 5550108 : use mo_constants, only : Psychro_dp
525 : use mo_pet, only : slope_satpressure
526 :
527 : implicit none
528 :
529 : class(riv_temp_type), intent(in) :: self
530 : !> air temperature in deg C
531 : real(dp), intent(in) :: air_temp
532 : !> net radiation in W * m-2
533 : real(dp), intent(in) :: netrad
534 :
535 : ! save slope of saturation vapor pressure curve
536 5550108 : real(dp) :: delta
537 :
538 5550108 : delta = slope_satpressure(air_temp) ! slope of saturation vapor pressure curve
539 : ! latent heat for priestley taylor PET formula
540 5550108 : lat_heat = self%pt_a_water * delta / (Psychro_dp + delta) * max(netrad, 0._dp)
541 :
542 5550108 : end function get_lat_heat
543 :
544 : !> \brief sensible heat flux of \ref riv_temp_type
545 : !> \return sensible heat flux
546 5550108 : real(dp) function get_sens_heat(self, air_temp, riv_temp) result(sens_heat)
547 :
548 : implicit none
549 :
550 : class(riv_temp_type), intent(in) :: self
551 : real(dp), intent(in) :: air_temp !< air temperature in [deg C]
552 : real(dp), intent(in) :: riv_temp !< river temperature in [deg C]
553 :
554 : ! sensible heat flux resulting for temp. diff.: river <-> air
555 5550108 : sens_heat = self%turb_heat_ex_coeff * (riv_temp - air_temp)
556 :
557 11100216 : end function get_sens_heat
558 :
559 : !> \brief get complete energy source of \ref riv_temp_type at given cell
560 : !> \return energy IO
561 5550108 : real(dp) function get_E_IO(self, riv_temp, cell) result(E_IO)
562 :
563 5550108 : use mo_constants, only : T0_dp, cp_w_dp
564 : use mo_mhm_constants, only : H2Odens
565 :
566 : implicit none
567 :
568 : class(riv_temp_type), intent(in) :: self
569 : !> given river temperature in K to calculate heat fluxes
570 : real(dp), intent(in) :: riv_temp
571 : !> cell index in the current domain
572 : integer(i4), intent(in) :: cell
573 : !> net radiation calc from short and longwave radiation in/out
574 5550108 : real(dp) :: netrad, sens_heat, lat_heat
575 :
576 : ! net radiation
577 5550108 : netrad = self%L11_srad_net(cell) + self%L11_lrad_in(cell) - self%get_lrad_out(riv_temp)
578 : ! sensible heat flux
579 5550108 : sens_heat = self%get_sens_heat(self%L11_air_temp(cell), riv_temp - T0_dp)
580 : ! latent heat flux
581 5550108 : lat_heat = self%get_lat_heat(self%L11_air_temp(cell), netrad)
582 :
583 : ! convert energy flux [W m-2] to [K m s-1]
584 : ! needs to be multiplied with river-area to result in temp-energy flux [K m3 s-1]
585 5550108 : E_IO = (netrad - sens_heat - lat_heat) / (H2Odens * cp_w_dp)
586 :
587 5550108 : end function get_E_IO
588 :
589 : !> \brief execute the temperature routing of \ref riv_temp_type
590 6552 : subroutine L11_routing_E( &
591 : self, &
592 : nLinks, &
593 13104 : netPerm, &
594 6552 : netLink_fromN, &
595 6552 : netLink_toN, &
596 13104 : netLink_C1, &
597 6552 : netLink_C2, &
598 : nInflowGauges, &
599 13104 : InflowHeadwater, &
600 6552 : InflowNodeList, &
601 6552 : sink_cells, &
602 6552 : L11_qTR, &
603 6552 : L11_Qmod &
604 : )
605 5550108 : use mo_constants, only : T0_dp
606 :
607 : implicit none
608 :
609 : class(riv_temp_type), intent(inout) :: self
610 : !> number of stream segment (reaches)
611 : integer(i4), intent(in) :: nLinks
612 : !> routing order of a given domain (permutation)
613 : integer(i4), dimension(:), intent(in) :: netPerm
614 : !> from node
615 : integer(i4), dimension(:), intent(in) :: netLink_fromN
616 : !> to node
617 : integer(i4), dimension(:), intent(in) :: netLink_toN
618 : !> routing parameter C1 (\cite CMM1988 p. 25-41)
619 : real(dp), dimension(:), intent(in) :: netLink_C1
620 : !> routing parameters C2 (id)
621 : real(dp), dimension(:), intent(in) :: netLink_C2
622 : !> [-] number of inflow points
623 : integer(i4), intent(in) :: nInflowGauges
624 : !> [-] if to consider headwater cells of inflow gauge
625 : logical, dimension(:), intent(in) :: InflowHeadwater
626 : !> [-] L11 ID of inflow points
627 : integer(i4), dimension(:), intent(in) :: InflowNodeList
628 : ! [-] sink nodes
629 : integer(i4), dimension(:), intent(in) :: sink_cells
630 : !> [m3 s-1] Transformed outflow leaving node I at current timestep(Muskingum)
631 : real(dp), intent(in), dimension(:) :: L11_qTR
632 : !> [m3 s-1] Simulated routed discharge
633 : real(dp), intent(in), dimension(:) :: L11_Qmod
634 :
635 : integer(i4) :: i, k, m, iNode, tNode, L11in, L11to
636 6552 : real(dp) :: E_IO_in, E_IO, T_est, T_rout
637 :
638 : ! initialize total input at point time (t+1) in all nodes
639 229320 : self%netNode_E_IN(self%s11 : self%e11, 2) = 0.0_dp
640 :
641 222768 : riv_loop: do k = 1, nLinks
642 : ! get LINK routing order -> i
643 216216 : i = netPerm(k)
644 216216 : iNode = netLink_fromN(i)
645 216216 : tNode = netLink_toN(i)
646 : ! indices for concaternated domains
647 216216 : L11in = iNode + self%s11 - 1_i4
648 216216 : L11to = tNode + self%s11 - 1_i4
649 :
650 : ! accumulate all inputs in iNode
651 216216 : self%netNode_E_IN(L11in, 2) = self%netNode_E_IN(L11in, 2) + self%netNode_E_out(L11in)
652 : ! calculate the river temp at the upstream node
653 216216 : self%river_temp(L11in) = self%netNode_E_IN(L11in, 2) / L11_Qmod(iNode) - T0_dp
654 :
655 : ! first guess for the routed temperature is the upstream temp
656 216216 : T_est = self%river_temp(L11in) + T0_dp
657 : ! calculate the meteo energy source at the upstream node
658 216216 : E_IO_in = self%get_E_IO(T_est, iNode)
659 :
660 : ! initialize iteration
661 216216 : call self%init_iter()
662 : ! perform iteration
663 5333892 : iterloop: do m=1, self%max_iter
664 : ! meteo energy flux depending on the estimated routed temperature (avg. IN- and OUT-node)
665 5333892 : E_IO = (E_IO_in + self%get_E_IO(T_est, tNode)) * 0.5_dp * self%L11_riv_areas(L11in)
666 : ! routing iNode
667 5333892 : self%netNode_E_R(L11in, 2) = self%netNode_E_R(L11in, 1) &
668 5333892 : + netLink_C1(i) * (self%netNode_E_IN(L11in, 1) - self%netNode_E_R(L11in, 1)) &
669 10667784 : + netLink_C2(i) * (self%netNode_E_IN(L11in, 2) + E_IO - self%netNode_E_IN(L11in, 1))
670 : ! calculate the routed temperature
671 5333892 : T_rout = self%netNode_E_R(L11in, 2) / L11_qTR(iNode)
672 : ! check for convergence
673 5333892 : if ( abs(T_rout - T_est) < self%delta_iter ) exit iterloop
674 : ! get new estimate for next iteration
675 5333892 : call self%next_iter(T_est, T_rout)
676 : end do iterloop
677 :
678 : ! add the meteo energy flux to the accumulated incoming energy at the IN-node
679 216216 : self%netNode_E_IN(L11in, 2) = self%netNode_E_IN(L11in, 2) + E_IO
680 :
681 : ! check if the inflow from upstream cells should be deactivated
682 216216 : if (nInflowGauges .GT. 0) then
683 0 : do i = 1, nInflowGauges
684 : ! check if downstream Node (tNode) is inflow gauge and headwaters should be ignored
685 0 : if ((tNode == InflowNodeList(i)) .AND. (.NOT. InflowHeadwater(i))) &
686 0 : self%netNode_E_R(L11in, 2) = 0.0_dp
687 : end do
688 : end if
689 :
690 : ! add routed temp energy to downstream node
691 222768 : self%netNode_E_IN(L11to, 2) = self%netNode_E_IN(L11to, 2) + self%netNode_E_R(L11in, 2)
692 :
693 : end do riv_loop
694 :
695 : ! Accumulate all inputs at sinks
696 : ! calculate riv temp at sinks
697 13104 : do k = 1, size(sink_cells)
698 6552 : tNode = sink_cells(k)
699 6552 : L11to = tNode + self%s11 - 1_i4
700 6552 : self%netNode_E_IN(L11to, 2) = self%netNode_E_IN(L11to, 2) + self%netNode_E_out(L11to)
701 13104 : self%river_temp(L11to) = self%netNode_E_IN(L11to, 2) / L11_Qmod(tNode) - T0_dp
702 : end do
703 :
704 : ! backflow t-> t-1
705 229320 : self%netNode_E_R(self%s11 : self%e11, 1) = self%netNode_E_R(self%s11 : self%e11, 2)
706 229320 : self%netNode_E_IN(self%s11 : self%e11, 1) = self%netNode_E_IN(self%s11 : self%e11, 2)
707 :
708 6552 : end subroutine L11_routing_E
709 :
710 : !> \brief initialize iterative solver of \ref riv_temp_type
711 216216 : subroutine init_iter(self)
712 :
713 : implicit none
714 :
715 : class(riv_temp_type), intent(inout) :: self
716 :
717 : ! first iteration is used to determine the search direction for the interval
718 216216 : self%first_iter = .true.
719 : ! after the interval is found, we will perform a biseciton to nail down the temperature
720 216216 : self%bisect_iter = .false.
721 :
722 6552 : end subroutine init_iter
723 :
724 : !> \brief execute next iteration with iterative solver of \ref riv_temp_type
725 5117676 : subroutine next_iter(self, T_est, T_rout)
726 :
727 : implicit none
728 :
729 : class(riv_temp_type), intent(inout) :: self
730 : real(dp), intent(inout) :: T_est !< estimated river temperature
731 : real(dp), intent(in) :: T_rout !< calculated (routed) river temperature
732 :
733 : ! before performing a bisection we need to search for the interval (with given step-size)
734 5117676 : if ( .not. self%bisect_iter ) then
735 : ! to determine the direction of interval search, we have to wait for the first iteration
736 641761 : if ( self%first_iter ) then
737 216216 : self%first_iter = .false.
738 : ! determine search direction
739 216216 : self%up_iter = ( T_rout > T_est )
740 : end if
741 : ! check if we need to take another step
742 641761 : if ( self%up_iter .eqv. ( T_rout > T_est ) ) then
743 425545 : if ( self%up_iter ) then
744 203382 : T_est = T_est + self%step_iter
745 : else
746 222163 : T_est = T_est - self%step_iter
747 : end if
748 : ! if direction changes, we have found the interval for bisection
749 : else
750 : ! start bisection in next iteration
751 216216 : self%bisect_iter = .true.
752 : ! set interval bounds for bisection
753 216216 : if ( self%up_iter ) then
754 103220 : self%up_bnd_iter = T_est
755 103220 : self%low_bnd_iter = T_est - self%step_iter
756 : else
757 112996 : self%up_bnd_iter = T_est + self%step_iter
758 112996 : self%low_bnd_iter = T_est
759 : end if
760 : ! set estimation to interval center
761 216216 : T_est = (self%up_bnd_iter + self%low_bnd_iter) / 2.0_dp
762 : end if
763 : ! perform bisection
764 : else
765 : ! if routed temp is lower than estimated, use lower interval
766 4475915 : if ( T_rout < T_est ) then
767 2233019 : self%up_bnd_iter = T_est
768 : ! if routed temp is higher than estimated, use upper interval
769 : else
770 2242896 : self%low_bnd_iter = T_est
771 : end if
772 : ! set estimation to interval center
773 4475915 : T_est = (self%up_bnd_iter + self%low_bnd_iter) / 2.0_dp
774 : end if
775 :
776 5333892 : end subroutine next_iter
777 :
778 5117676 : end module mo_mrm_riv_temp_class
|
