Line data Source code
1 : !> \file mo_mhm_interface_run.f90
2 : !> \brief \copybrief mo_mhm_interface_run
3 : !> \details \copydetails mo_mhm_interface_run
4 :
5 : !> \brief Module providing interfaces for running preconfigured mHM.
6 : !> \details Interfaces to control the mHM run from outside (prepare domain, do timestep, ...).
7 : !> \authors Sebastian Mueller, Matthias Kelbling
8 : !> \version 0.1
9 : !> \date Jan 2022
10 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
11 : !! mHM is released under the LGPLv3+ license \license_note
12 : !> \ingroup f_mhm
13 : module mo_mhm_interface_run
14 :
15 : ! forces
16 : use mo_kind, only: i4, dp
17 : use mo_message, only: message, error_message
18 : use mo_string_utils, only : num2str
19 : ! mhm
20 : use mo_common_run_variables, only : run_cfg
21 : use mo_optimization_types, only : optidata_sim
22 : use mo_common_datetime_type, only : datetimeinfo
23 : use mo_common_mHM_mRM_variables, only : &
24 : resolutionRouting, &
25 : LCyearId, &
26 : mhmFileRestartIn, &
27 : mrmFileRestartIn, &
28 : nTstepDay, &
29 : optimize, &
30 : read_restart, &
31 : simPer, &
32 : timeStep, &
33 : warmingDays, &
34 : c2TSTu, &
35 : restart_reset_fluxes_states
36 : use mo_common_variables, only : &
37 : global_parameters, &
38 : level0, &
39 : level1, &
40 : domainMeta, &
41 : processMatrix
42 : use mo_global_variables, only : &
43 : meteo_handler, &
44 : L1_Throughfall, &
45 : L1_aETCanopy, &
46 : L1_aETSealed, &
47 : L1_aETSoil, &
48 : L1_baseflow, &
49 : L1_fastRunoff, &
50 : L1_infilSoil, &
51 : L1_inter, &
52 : L1_melt, &
53 : L1_neutrons, &
54 : L1_percol, &
55 : L1_pet_calc, &
56 : L1_temp_calc, &
57 : L1_prec_calc, &
58 : L1_preEffect, &
59 : L1_rain, &
60 : L1_runoffSeal, &
61 : L1_satSTW, &
62 : L1_sealSTW, &
63 : L1_slowRunoff, &
64 : L1_snow, &
65 : L1_snowPack, &
66 : L1_soilMoist, &
67 : L1_total_runoff, &
68 : L1_unsatSTW, &
69 : L1_twsaObs, &
70 : L1_etObs, &
71 : L1_smObs, &
72 : L1_neutronsObs, &
73 : evap_coeff, &
74 : nSoilHorizons_sm_input, &
75 : neutron_integral_AFast, &
76 : outputFlxState, &
77 : timeStep_model_outputs, &
78 : BFI_qBF_sum, &
79 : BFI_qT_sum
80 : use mo_init_states, only : variables_default_init, fluxes_states_default_init
81 : use mo_julian, only : caldat, julday
82 : use mo_string_utils, only : num2str
83 : use mo_mhm, only : mhm
84 : use mo_restart, only : read_restart_states
85 : use mo_write_fluxes_states, only : mHM_updateDataset, mHM_OutputDataset
86 : use mo_mrm_write_fluxes_states, only : mRM_updateDataset, mRM_OutputDataset, GW_OutputDataset, GW_updateDataset
87 : use mo_constants, only : HourSecs
88 : use mo_common_variables, only : resolutionHydrology
89 : use mo_mrm_global_variables, only : &
90 : InflowGauge, &
91 : L11_C1, &
92 : L11_C2, &
93 : L11_L1_Id, &
94 : L11_TSrout, &
95 : L11_fromN, &
96 : L11_length, &
97 : sink_cells, &
98 : L11_nLinkFracFPimp, &
99 : L11_nOutlets, &
100 : L11_netPerm, &
101 : L11_qMod, &
102 : L11_qOUT, &
103 : L11_qTIN, &
104 : L11_qTR, &
105 : L11_slope, &
106 : L11_toN, &
107 : L1_L11_Id, &
108 : domain_mrm, &
109 : level11, &
110 : mRM_runoff, &
111 : outputFlxState_mrm, &
112 : timeStep_model_outputs_mrm, &
113 : gw_coupling, &
114 : L0_river_head_mon_sum, &
115 : riv_temp_pcs
116 : use mo_mpr_global_variables, only : &
117 : L1_HarSamCoeff, &
118 : L1_PrieTayAlpha, &
119 : L1_aeroResist, &
120 : L1_alpha, &
121 : L1_degDay, &
122 : L1_degDayInc, &
123 : L1_degDayMax, &
124 : L1_degDayNoPre, &
125 : L1_fAsp, &
126 : L1_fRoots, &
127 : L1_fSealed, &
128 : L1_jarvis_thresh_c1, &
129 : L1_kBaseFlow, &
130 : L1_kPerco, &
131 : L1_kSlowFlow, &
132 : L1_karstLoss, &
133 : L1_kfastFlow, &
134 : L1_maxInter, &
135 : L1_petLAIcorFactor, &
136 : L1_sealedThresh, &
137 : L1_soilMoistExp, &
138 : L1_soilMoistFC, &
139 : L1_soilMoistSat, &
140 : L1_surfResist, &
141 : L1_tempThresh, &
142 : L1_unsatThresh, &
143 : L1_wiltingPoint, &
144 : L1_No_Count, &
145 : L1_bulkDens, &
146 : L1_latticeWater, &
147 : L1_COSMICL3, &
148 : HorizonDepth_mHM, &
149 : nSoilHorizons_mHM
150 : use mo_mrm_init, only : variables_default_init_routing, fluxes_states_default_init_routing
151 : use mo_mrm_mpr, only : mrm_update_param
152 : use mo_mrm_restart, only : mrm_read_restart_states
153 : use mo_mrm_routing, only : mrm_routing
154 : use mo_utils, only : ge
155 : use mo_mrm_river_head, only: calc_river_head
156 : use mo_mpr_eval, only : mpr_eval
157 :
158 : contains
159 :
160 : !> \brief prepare single run of mHM
161 122 : subroutine mhm_interface_run_prepare(parameterset, opti_domain_indices, runoff_present, BFI_present)
162 : implicit none
163 : !> a set of global parameter (gamma) to run mHM, DIMENSION [no. of global_Parameters]
164 : real(dp), dimension(:), optional, intent(in) :: parameterset
165 : !> selected domains for optimization
166 : integer(i4), dimension(:), optional, intent(in) :: opti_domain_indices
167 : !> whether runoff is present
168 : logical, optional, intent(in) :: runoff_present
169 : !> whether BFI is present
170 : logical, optional, intent(in) :: BFI_present
171 :
172 : integer(i4) :: i, iDomain, domainID
173 :
174 : ! run_cfg%output_runoff = .false. by default
175 61 : if (present(runoff_present)) run_cfg%output_runoff = runoff_present
176 : ! run_cfg%output_BFI = .false. by default
177 61 : if (present(BFI_present)) run_cfg%output_BFI = BFI_present
178 :
179 : ! store current parameter set
180 183 : allocate(run_cfg%parameterset(size(global_parameters, dim=1)))
181 61 : if (.not. present(parameterset) .and. optimize) then
182 0 : call error_message("mhm_interface_run_prepare: Can't optimize without parameter!")
183 61 : else if (.not. present(parameterset)) then
184 0 : run_cfg%parameterset = global_parameters(:, 3)
185 : else
186 3283 : run_cfg%parameterset = parameterset
187 : end if
188 :
189 : ! store domain indices for domain loop
190 61 : if (optimize .and. present(opti_domain_indices)) then
191 28 : run_cfg%nDomains = size(opti_domain_indices)
192 84 : allocate(run_cfg%domain_indices(run_cfg%nDomains))
193 98 : run_cfg%domain_indices = opti_domain_indices
194 : else
195 33 : run_cfg%nDomains = domainMeta%nDomains
196 99 : allocate(run_cfg%domain_indices(run_cfg%nDomains))
197 186 : run_cfg%domain_indices = [(i, i=1, run_cfg%nDomains)]
198 : end if
199 :
200 : !----------------------------------------------------------
201 : ! Check optionals and initialize
202 : !----------------------------------------------------------
203 61 : if (run_cfg%output_runoff) then
204 50 : do i = 1, run_cfg%nDomains
205 25 : iDomain = run_cfg%get_domain_index(i)
206 25 : domainID = domainMeta%indices(iDomain)
207 50 : if (.not. domainMeta%doRouting(iDomain)) then
208 : call error_message("***ERROR: runoff for domain", trim(num2str(domainID)),&
209 0 : "can not be produced, since routing process is off in Process Matrix")
210 : end if
211 : end do
212 : end if
213 :
214 : ! prepare BFI calculation
215 61 : if (run_cfg%output_BFI) then
216 0 : allocate(BFI_qBF_sum(run_cfg%nDomains))
217 0 : allocate(BFI_qT_sum(run_cfg%nDomains))
218 0 : BFI_qBF_sum = 0.0_dp
219 0 : BFI_qT_sum = 0.0_dp
220 : end if
221 :
222 61 : if (read_restart) then
223 2 : do i = 1, run_cfg%nDomains
224 1 : iDomain = run_cfg%get_domain_index(i)
225 1 : domainID = domainMeta%indices(iDomain)
226 : ! this reads the eff. parameters and optionally the states and fluxes
227 2 : call read_restart_states(iDomain, domainID, mhmFileRestartIn(iDomain))
228 : end do
229 1 : if (restart_reset_fluxes_states) then
230 0 : call message(' Resetting mHM states and fluxes from restart files ...')
231 0 : call fluxes_states_default_init()
232 : end if
233 : else
234 60 : call variables_default_init()
235 60 : call mpr_eval(run_cfg%parameterset)
236 60 : if (processMatrix(8, 1) > 0) call variables_default_init_routing()
237 : end if
238 :
239 305 : allocate(run_cfg%L1_fNotSealed(size(L1_fSealed, 1), size(L1_fSealed, 2), size(L1_fSealed, 3)))
240 15650 : run_cfg%L1_fNotSealed = 1.0_dp - L1_fSealed
241 :
242 61 : end subroutine mhm_interface_run_prepare
243 :
244 : !> \brief get number of domains for looping
245 61 : subroutine mhm_interface_run_get_ndomains(ndomains)
246 : implicit none
247 : integer(i4), intent(inout) :: ndomains !< number of domains
248 61 : ndomains = run_cfg%nDomains
249 61 : end subroutine mhm_interface_run_get_ndomains
250 :
251 : !> \brief prepare single domain to run mHM on
252 164 : subroutine mhm_interface_run_prepare_domain(domain, etOptiSim, twsOptiSim, neutronsOptiSim, smOptiSim)
253 : implicit none
254 : !> domain loop counter
255 : integer(i4), intent(in), optional :: domain
256 : !> returns soil moisture time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
257 : type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim
258 : !> dim1=ncells, dim2=time
259 : type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim
260 : !> returns evapotranspiration time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
261 : type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim
262 : !> returns tws time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
263 : type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim
264 :
265 : integer(i4) :: iDomain, domainID
266 :
267 : ! get domain index
268 82 : run_cfg%selected_domain = 1_i4
269 82 : if (present(domain)) run_cfg%selected_domain = domain
270 82 : iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
271 82 : domainID = domainMeta%indices(iDomain)
272 :
273 : ! evapotranspiration optimization
274 82 : if (present(etOptiSim)) call etOptiSim(iDomain)%init(L1_etObs(iDomain))
275 : ! total water storage optimization
276 82 : if (present(twsOptiSim)) call twsOptiSim(iDomain)%init(L1_twsaObs(iDomain))
277 : ! neutrons optimization
278 82 : if (present(neutronsOptiSim)) call neutronsOptiSim(iDomain)%init(L1_neutronsObs(iDomain))
279 : ! sm optimization
280 82 : if (present(smOptiSim)) call smOptiSim(iDomain)%init(L1_smObs(iDomain))
281 :
282 : ! get Domain information
283 82 : run_cfg%nCells = level1(iDomain)%nCells
284 82 : run_cfg%mask1 => level1(iDomain)%mask
285 82 : run_cfg%s1 = level1(iDomain)%iStart
286 82 : run_cfg%e1 = level1(iDomain)%iEnd
287 :
288 82 : if (domainMeta%doRouting(iDomain)) then
289 : ! ----------------------------------------
290 : ! initialize factor between routing resolution and hydrologic model resolution
291 : ! ----------------------------------------
292 46 : run_cfg%tsRoutFactor = 1_i4
293 138 : allocate(run_cfg%InflowDischarge(size(InflowGauge%Q, dim = 2)))
294 95 : run_cfg%InflowDischarge = 0._dp
295 :
296 : ! read states from restart
297 46 : if (read_restart) then
298 1 : call mrm_read_restart_states(iDomain, domainID, mrmFileRestartIn(iDomain))
299 1 : if (restart_reset_fluxes_states) then
300 0 : call message(' Resetting mRM states and fluxes from restart files for domain ', num2str(iDomain), ' ...')
301 0 : call fluxes_states_default_init_routing(iDomain)
302 : end if
303 : end if
304 : ! get Domain information at L11 if routing is activated
305 46 : run_cfg%s11 = level11(iDomain)%iStart
306 46 : run_cfg%e11 = level11(iDomain)%iEnd
307 46 : run_cfg%mask11 => level11(iDomain)%mask
308 :
309 : ! initialize routing parameters (has to be called for routing option 2)
310 46 : if ((processMatrix(8, 1) .eq. 2) .or. (processMatrix(8, 1) .eq. 3)) &
311 : call mrm_update_param( &
312 10 : iDomain, run_cfg%parameterset(processMatrix(8, 3) - processMatrix(8, 2) + 1 : processMatrix(8, 3)))
313 : ! initialize variable for runoff for routing
314 138 : allocate(run_cfg%RunToRout(run_cfg%e1 - run_cfg%s1 + 1))
315 2180 : run_cfg%RunToRout = 0._dp
316 :
317 46 : if ( riv_temp_pcs%active ) then
318 : ! set indices for current L11 domain
319 1 : riv_temp_pcs%s11 = run_cfg%s11
320 1 : riv_temp_pcs%e11 = run_cfg%e11
321 : ! allocate current L1 lateral components
322 1 : call riv_temp_pcs%alloc_lateral(run_cfg%nCells)
323 : end if
324 : end if
325 :
326 : ! init datetime variable
327 82 : call run_cfg%domainDateTime%init(iDomain)
328 :
329 : ! init time-loop
330 82 : run_cfg%time_step = 0_i4
331 :
332 143 : end subroutine mhm_interface_run_prepare_domain
333 :
334 : !> \brief check if current time loop is finished
335 2020272 : subroutine mhm_interface_run_finished(time_loop_finished)
336 : implicit none
337 : logical, intent(inout) :: time_loop_finished !< flag to indicate end of timeloop
338 2020272 : time_loop_finished = run_cfg%time_step == run_cfg%domainDateTime%nTimeSteps
339 82 : end subroutine mhm_interface_run_finished
340 :
341 : !> \brief do one time-step on current domain
342 2020272 : subroutine mhm_interface_run_do_time_step()
343 : implicit none
344 :
345 : integer(i4) :: iDomain, domainID, tt, jj, s1, e1
346 :
347 : ! increment time step count (first input is 0)
348 2020272 : run_cfg%time_step = run_cfg%time_step + 1_i4
349 : ! current time counter and domain indices
350 2020272 : tt = run_cfg%time_step
351 2020272 : s1 = run_cfg%s1
352 2020272 : e1 = run_cfg%e1
353 2020272 : iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
354 2020272 : domainID = domainMeta%indices(iDomain)
355 :
356 2020272 : call run_cfg%domainDateTime%update_LAI_timestep()
357 :
358 : ! update the meteo-handler
359 : call meteo_handler%update_timestep( &
360 : tt=tt, &
361 : time=run_cfg%domainDateTime%newTime - 0.5_dp, &
362 : iDomain=iDomain, &
363 : level1=level1, &
364 2020272 : simPer=simPer)
365 :
366 : ! get corrected pet
367 0 : call meteo_handler%get_corrected_pet(pet_calc=L1_pet_calc(s1 : e1), &
368 : ! pet calculation dependencies
369 0 : petLAIcorFactorL1 = L1_petLAIcorFactor(s1 : e1, run_cfg%domainDateTime%iLAI, run_cfg%domainDateTime%yId), &
370 0 : fAsp = L1_fAsp(s1 : e1, 1, 1), &
371 0 : HarSamCoeff = L1_HarSamCoeff(s1 : e1, 1, 1), &
372 0 : latitude = pack(level1(iDomain)%y, level1(iDomain)%mask), &
373 0 : PrieTayAlpha = L1_PrieTayAlpha(s1 : e1, run_cfg%domainDateTime%iLAI, 1), &
374 0 : aeroResist = L1_aeroResist(s1 : e1, run_cfg%domainDateTime%iLAI, run_cfg%domainDateTime%yId), &
375 2020272 : surfResist = L1_surfResist(s1 : e1, run_cfg%domainDateTime%iLAI, 1))
376 :
377 : ! get temperature and precipitation
378 2020272 : call meteo_handler%get_temp(temp_calc=L1_temp_calc(s1 : e1))
379 2020272 : call meteo_handler%get_prec(prec_calc=L1_prec_calc(s1 : e1))
380 :
381 : ! -------------------------------------------------------------------------
382 : ! ARGUMENT LIST KEY FOR mHM
383 : ! -------------------------------------------------------------------------
384 : ! C CONFIGURATION
385 : ! F FORCING DATA L2
386 : ! Q INFLOW FROM UPSTREAM AREAS
387 : ! L0 MORPHOLOGIC DATA L0
388 : ! L1 MORPHOLOGIC DATA L1
389 : ! L11 MORPHOLOGIC DATA L11
390 : ! P GLOBAL PARAMETERS
391 : ! E EFFECTIVE PARAMETER FIELDS (L1, L11 levels)
392 : ! S STATE VARIABLES L1
393 : ! X FLUXES (L1, L11 levels)
394 : ! --------------------------------------------------------------------------
395 : call mhm( &
396 : read_states = read_restart, & ! IN C
397 : tt = tt, &
398 : time = run_cfg%domainDateTime%newTime - 0.5_dp, &
399 : processMatrix = processMatrix, &
400 : horizon_depth = HorizonDepth_mHM, & ! IN C
401 : nCells1 = run_cfg%nCells, &
402 : nHorizons_mHM = nSoilHorizons_mHM, &
403 : c2TSTu = c2TSTu, & ! IN C
404 : neutron_integral_AFast = neutron_integral_AFast, & ! IN C
405 : evap_coeff = evap_coeff, &
406 0 : fSealed1 = L1_fSealed(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT L1
407 0 : interc = L1_inter(s1 : e1), &
408 0 : snowpack = L1_snowPack(s1 : e1), &
409 0 : sealedStorage = L1_sealSTW(s1 : e1), & ! INOUT S
410 0 : soilMoisture = L1_soilMoist(s1 : e1, :), &
411 0 : unsatStorage = L1_unsatSTW(s1 : e1), &
412 0 : satStorage = L1_satSTW(s1 : e1), & ! INOUT S
413 0 : neutrons = L1_neutrons(s1 : e1), & ! INOUT S
414 0 : pet_calc = L1_pet_calc(s1 : e1), & ! INOUT X
415 0 : temp_calc = L1_temp_calc(s1 : e1), & ! INOUT X
416 0 : prec_calc = L1_prec_calc(s1 : e1), & ! INOUT X
417 0 : aet_soil = L1_aETSoil(s1 : e1, :), &
418 0 : aet_canopy = L1_aETCanopy(s1 : e1), &
419 0 : aet_sealed = L1_aETSealed(s1 : e1), & ! INOUT X
420 0 : baseflow = L1_baseflow(s1 : e1), &
421 0 : infiltration = L1_infilSoil(s1 : e1, :), &
422 0 : fast_interflow = L1_fastRunoff(s1 : e1), & ! INOUT X
423 0 : melt = L1_melt(s1 : e1), &
424 0 : perc = L1_percol(s1 : e1), &
425 0 : prec_effect = L1_preEffect(s1 : e1), &
426 0 : rain = L1_rain(s1 : e1), & ! INOUT X
427 0 : runoff_sealed = L1_runoffSeal(s1 : e1), &
428 0 : slow_interflow = L1_slowRunoff(s1 : e1), &
429 0 : snow = L1_snow(s1 : e1), & ! INOUT X
430 0 : throughfall = L1_Throughfall(s1 : e1), &
431 0 : total_runoff = L1_total_runoff(s1 : e1), & ! INOUT X
432 : ! MPR
433 0 : alpha = L1_alpha(s1 : e1, 1, run_cfg%domainDateTime%yId), &
434 0 : deg_day_incr = L1_degDayInc(s1 : e1, 1, run_cfg%domainDateTime%yId), &
435 0 : deg_day_max = L1_degDayMax(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
436 0 : deg_day_noprec = L1_degDayNoPre(s1 : e1, 1, run_cfg%domainDateTime%yId), &
437 0 : deg_day = L1_degDay(s1 : e1, 1, 1), & ! INOUT E1
438 0 : frac_roots = L1_fRoots(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
439 0 : interc_max = L1_maxInter(s1 : e1, run_cfg%domainDateTime%iLAI, 1), &
440 0 : karst_loss = L1_karstLoss(s1 : e1, 1, 1), & ! INOUT E1
441 0 : k0 = L1_kFastFlow(s1 : e1, 1, run_cfg%domainDateTime%yId), &
442 0 : k1 = L1_kSlowFlow(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
443 0 : k2 = L1_kBaseFlow(s1 : e1, 1, run_cfg%domainDateTime%yId), &
444 0 : kp = L1_kPerco(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
445 0 : soil_moist_FC = L1_soilMoistFC(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
446 0 : soil_moist_sat = L1_soilMoistSat(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
447 0 : soil_moist_exponen = L1_soilMoistExp(s1 : e1, :, run_cfg%domainDateTime%yId), &
448 0 : jarvis_thresh_c1 = L1_jarvis_thresh_c1(s1 : e1, 1, 1), & ! INOUT E1
449 0 : temp_thresh = L1_tempThresh(s1 : e1, 1, run_cfg%domainDateTime%yId), &
450 0 : unsat_thresh = L1_unsatThresh(s1 : e1, 1, 1), & ! INOUT E1
451 0 : water_thresh_sealed = L1_sealedThresh(s1 : e1, 1, 1), & ! INOUT E1
452 0 : wilting_point = L1_wiltingPoint(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
453 : ! >> neutron count
454 0 : No_count = L1_No_Count(s1:e1, 1, 1), & ! INOUT E1
455 0 : bulkDens = L1_bulkDens(s1:e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
456 0 : latticeWater = L1_latticeWater(s1:e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
457 0 : COSMICL3 = L1_COSMICL3(s1:e1, :, run_cfg%domainDateTime%yId) & ! INOUT E1
458 2020272 : )
459 :
460 : ! call mRM routing
461 2020272 : run_cfg%doRoute = .false.
462 2020272 : if (domainMeta%doRouting(iDomain)) then
463 : ! set discharge timestep
464 779928 : run_cfg%iDischargeTS = ceiling(real(tt, dp) / real(nTstepDay, dp))
465 : ! set input variables for routing
466 779928 : if (processMatrix(8, 1) .eq. 1) then
467 : ! >>>
468 : ! >>> original Muskingum routing, executed every time
469 : ! >>>
470 495096 : run_cfg%doRoute = .True.
471 495096 : run_cfg%tsRoutFactorIn = 1._dp
472 495096 : run_cfg%timestep_rout = timestep
473 25781256 : run_cfg%RunToRout = L1_total_runoff(s1 : e1) ! runoff [mm TS-1] mm per timestep
474 1524600 : run_cfg%InflowDischarge = InflowGauge%Q(run_cfg%iDischargeTS, :) ! inflow discharge in [m3 s-1]
475 : !
476 284832 : else if ((processMatrix(8, 1) .eq. 2) .or. &
477 : (processMatrix(8, 1) .eq. 3)) then
478 : ! >>>
479 : ! >>> adaptive timestep
480 : ! >>>
481 : run_cfg%doRoute = .False.
482 : ! calculate factor
483 284832 : run_cfg%tsRoutFactor = L11_tsRout(iDomain) / (timestep * HourSecs)
484 : ! print *, 'routing factor: ', tsRoutFactor
485 : ! prepare routing call
486 284832 : if (run_cfg%tsRoutFactor .lt. 1._dp) then
487 : ! ----------------------------------------------------------------
488 : ! routing timesteps are shorter than hydrologic time steps
489 : ! ----------------------------------------------------------------
490 : ! set all input variables
491 0 : run_cfg%tsRoutFactorIn = run_cfg%tsRoutFactor
492 0 : run_cfg%RunToRout = L1_total_runoff(s1 : e1) ! runoff [mm TS-1] mm per timestep
493 0 : run_cfg%InflowDischarge = InflowGauge%Q(run_cfg%iDischargeTS, :) ! inflow discharge in [m3 s-1]
494 0 : run_cfg%timestep_rout = timestep
495 0 : run_cfg%doRoute = .True.
496 : else
497 : ! ----------------------------------------------------------------
498 : ! routing timesteps are longer than hydrologic time steps
499 : ! ----------------------------------------------------------------
500 : ! set all input variables
501 284832 : run_cfg%tsRoutFactorIn = run_cfg%tsRoutFactor
502 : ! Runoff is accumulated in [mm]
503 10253952 : run_cfg%RunToRout = run_cfg%RunToRout + L1_total_runoff(s1 : e1)
504 854496 : run_cfg%InflowDischarge = run_cfg%InflowDischarge + InflowGauge%Q(run_cfg%iDischargeTS, :)
505 : ! reset tsRoutFactorIn if last period did not cover full period
506 284832 : if ((tt == run_cfg%domainDateTime%nTimeSteps) .and. (mod(tt, nint(run_cfg%tsRoutFactorIn)) /= 0_i4)) &
507 0 : run_cfg%tsRoutFactorIn = mod(tt, nint(run_cfg%tsRoutFactorIn))
508 284832 : if ((mod(tt, nint(run_cfg%tsRoutFactorIn)) .eq. 0_i4) .or. (tt .eq. run_cfg%domainDateTime%nTimeSteps)) then
509 : ! Inflow discharge is given as flow-rate and has to be converted to [m3]
510 278280 : run_cfg%InflowDischarge = run_cfg%InflowDischarge / run_cfg%tsRoutFactorIn
511 139140 : run_cfg%timestep_rout = timestep * nint(run_cfg%tsRoutFactorIn, i4)
512 139140 : run_cfg%doRoute = .True.
513 : end if
514 : end if
515 : end if
516 : ! prepare temperature routing
517 779928 : if ( riv_temp_pcs%active ) then
518 : ! init riv-temp from current air temp
519 13104 : if ( tt .eq. 1_i4 ) call riv_temp_pcs%init_riv_temp( &
520 0 : temp_air = L1_temp_calc(s1 : e1), &
521 0 : efecarea = level1(iDomain)%CellArea * 1.E-6_dp, &
522 0 : L1_L11_Id = L1_L11_Id(s1 : e1), &
523 0 : L11_areacell = level11(iDomain)%CellArea * 1.E-6_dp, &
524 0 : L11_L1_Id = L11_L1_Id(run_cfg%s11 : run_cfg%e11), &
525 0 : map_flag = ge(resolutionRouting(iDomain), resolutionHydrology(iDomain)) &
526 69 : )
527 : ! get riv-temp specific meteo arrays
528 13104 : call meteo_handler%get_ssrd(riv_temp_pcs%L1_ssrd_calc)
529 13104 : call meteo_handler%get_strd(riv_temp_pcs%L1_strd_calc)
530 13104 : call meteo_handler%get_tann(riv_temp_pcs%L1_tann_calc)
531 : ! accumulate source Energy at L1 level
532 : call riv_temp_pcs%acc_source_E( &
533 0 : fSealed_area_fraction = L1_fSealed(s1 : e1, 1, run_cfg%domainDateTime%yId), &
534 0 : fast_interflow = L1_fastRunoff(s1 : e1), &
535 0 : slow_interflow = L1_slowRunoff(s1 : e1), &
536 0 : baseflow = L1_baseflow(s1 : e1), &
537 0 : direct_runoff = L1_runoffSeal(s1 : e1), &
538 0 : temp_air = L1_temp_calc(s1 : e1) &
539 13104 : )
540 : ! if routing should be performed, scale source energy to L11 level
541 13104 : if ( run_cfg%doRoute ) call riv_temp_pcs%finalize_source_E( &
542 0 : efecarea = level1(iDomain)%CellArea * 1.E-6_dp, &
543 0 : L1_L11_Id = L1_L11_Id(s1 : e1), &
544 0 : L11_areacell = level11(iDomain)%CellArea * 1.E-6_dp, &
545 0 : L11_L1_Id = L11_L1_Id(run_cfg%s11 : run_cfg%e11), &
546 : timestep = run_cfg%timestep_rout, &
547 0 : map_flag = ge(resolutionRouting(iDomain), resolutionHydrology(iDomain)) &
548 452088 : )
549 : end if
550 : ! -------------------------------------------------------------------
551 : ! execute routing
552 : ! -------------------------------------------------------------------
553 779928 : if (run_cfg%doRoute) call mRM_routing(&
554 : ! general INPUT variables
555 : read_restart, &
556 : processMatrix(8, 1), & ! parse process Case to be used
557 0 : run_cfg%parameterset(processMatrix(8, 3) - processMatrix(8, 2) + 1 : processMatrix(8, 3)), & ! routing par.
558 : run_cfg%RunToRout, & ! runoff [mm TS-1] mm per timestep old: L1_total_runoff_in(run_cfg%s1:run_cfg%e1, tt), &
559 0 : level1(iDomain)%CellArea * 1.E-6_dp, &
560 0 : L1_L11_Id(s1 : e1), &
561 0 : level11(iDomain)%CellArea * 1.E-6_dp, &
562 0 : L11_L1_Id(run_cfg%s11 : run_cfg%e11), &
563 0 : L11_netPerm(run_cfg%s11 : run_cfg%e11), & ! routing order at L11
564 0 : L11_fromN(run_cfg%s11 : run_cfg%e11), & ! link source at L11
565 0 : L11_toN(run_cfg%s11 : run_cfg%e11), & ! link target at L11
566 0 : L11_nOutlets(iDomain), & ! number of outlets
567 : run_cfg%timestep_rout, & ! timestep of runoff to rout [h]
568 : run_cfg%tsRoutFactorIn, & ! simulate timestep in [h]
569 0 : level11(iDomain)%nCells, & ! number of Nodes
570 0 : domain_mrm(iDomain)%nInflowGauges, &
571 : domain_mrm(iDomain)%InflowGaugeIndexList(:), &
572 : domain_mrm(iDomain)%InflowGaugeHeadwater(:), &
573 : domain_mrm(iDomain)%InflowGaugeNodeList(:), &
574 : run_cfg%InflowDischarge, &
575 : domain_mrm(iDomain)%nGauges, &
576 : domain_mrm(iDomain)%gaugeIndexList(:), &
577 : domain_mrm(iDomain)%gaugeNodeList(:), &
578 0 : ge(resolutionRouting(iDomain), resolutionHydrology(iDomain)), &
579 : ! original routing specific input variables
580 0 : L11_length(run_cfg%s11 : run_cfg%e11 - 1), & ! link length
581 0 : L11_slope(run_cfg%s11 : run_cfg%e11 - 1), &
582 0 : L11_nLinkFracFPimp(run_cfg%s11 : run_cfg%e11, run_cfg%domainDateTime%yId), & ! fraction of impervious layer at L11 scale
583 0 : sink_cells(iDomain)%ids, &
584 : ! general INPUT/OUTPUT variables
585 0 : L11_C1(run_cfg%s11 : run_cfg%e11), & ! first muskingum parameter
586 0 : L11_C2(run_cfg%s11 : run_cfg%e11), & ! second muskigum parameter
587 0 : L11_qOUT(run_cfg%s11 : run_cfg%e11), & ! routed runoff flowing out of L11 cell
588 0 : L11_qTIN(run_cfg%s11 : run_cfg%e11, :), & ! inflow water into the reach at L11
589 0 : L11_qTR(run_cfg%s11 : run_cfg%e11, :), & !
590 0 : L11_qMod(run_cfg%s11 : run_cfg%e11), &
591 : mRM_runoff(tt, :) &
592 58038084 : )
593 : ! -------------------------------------------------------------------
594 : ! groundwater coupling
595 : ! -------------------------------------------------------------------
596 779928 : if (gw_coupling) then
597 : call calc_river_head(iDomain, L11_Qmod, L0_river_head_mon_sum)
598 : end if
599 : ! -------------------------------------------------------------------
600 : ! reset variables
601 : ! -------------------------------------------------------------------
602 779928 : if (processMatrix(8, 1) .eq. 1) then
603 : ! reset Input variables
604 1029504 : run_cfg%InflowDischarge = 0._dp
605 25286160 : run_cfg%RunToRout = 0._dp
606 284832 : else if ((processMatrix(8, 1) .eq. 2) .or. (processMatrix(8, 1) .eq. 3)) then
607 284832 : if ((.not. (run_cfg%tsRoutFactorIn .lt. 1._dp)) .and. run_cfg%doRoute) then
608 423972 : do jj = 1, nint(run_cfg%tsRoutFactorIn) ! BUG: this should start at 2
609 708804 : mRM_runoff(tt - jj + 1, :) = mRM_runoff(tt, :)
610 : end do
611 : ! reset Input variables
612 278280 : run_cfg%InflowDischarge = 0._dp
613 4869900 : run_cfg%RunToRout = 0._dp
614 : ! reset lateral fluxes and time-step counter if routing was done
615 139140 : if ( riv_temp_pcs%active ) call riv_temp_pcs%reset_timestep()
616 : end if
617 : ! if routing is done every time-step, reset river-temp time step
618 284832 : if (run_cfg%tsRoutFactor .lt. 1._dp .and. riv_temp_pcs%active ) call riv_temp_pcs%reset_timestep()
619 : end if
620 : end if
621 :
622 : ! output only for evaluation period
623 2020272 : run_cfg%domainDateTime%tIndex_out = (tt - warmingDays(iDomain) * nTstepDay) ! tt if write out of warming period
624 :
625 2020272 : call run_cfg%domainDateTime%increment()
626 :
627 : ! update the year-dependent domainDateTime%yId (land cover id)
628 2020272 : if (run_cfg%domainDateTime%is_new_year .and. tt < run_cfg%domainDateTime%nTimeSteps) then
629 165 : run_cfg%domainDateTime%yId = LCyearId(run_cfg%domainDateTime%year, iDomain)
630 : end if
631 :
632 : ! calculate BFI releated after warming days if wanted
633 2020272 : if ( run_cfg%output_BFI .and. (run_cfg%domainDateTime%tIndex_out > 0_i4) ) then
634 0 : BFI_qBF_sum(iDomain) = BFI_qBF_sum(iDomain) &
635 0 : + sum(L1_baseflow(s1 : e1) * level1(iDomain)%CellArea) / level1(iDomain)%nCells
636 0 : BFI_qT_sum(iDomain) = BFI_qT_sum(iDomain) &
637 0 : + sum(L1_total_runoff(s1 : e1) * level1(iDomain)%CellArea) / level1(iDomain)%nCells
638 : end if
639 :
640 4040544 : end subroutine mhm_interface_run_do_time_step
641 :
642 : !> \brief write output after current time-step
643 2020272 : subroutine mhm_interface_run_write_output()
644 : implicit none
645 : integer(i4) :: iDomain, tt, s0, e0
646 :
647 : ! get time step
648 2020272 : tt = run_cfg%time_step
649 :
650 : ! get domain index
651 2020272 : iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
652 :
653 2020272 : if ( (.not. optimize) .and. (run_cfg%domainDateTime%tIndex_out > 0_i4)) then
654 350760 : if (any(outputFlxState_mrm) .AND. (domainMeta%doRouting(iDomain))) then
655 :
656 35040 : if (run_cfg%domainDateTime%tIndex_out == 1) then
657 4 : run_cfg%nc_mrm = mRM_OutputDataset(iDomain, run_cfg%mask11)
658 : end if
659 :
660 : ! update Dataset (riv-temp as optional input)
661 35040 : if ( riv_temp_pcs%active ) then
662 : call mRM_updateDataset(run_cfg%nc_mrm, &
663 8760 : L11_qmod(run_cfg%s11 : run_cfg%e11), riv_temp_pcs%river_temp(riv_temp_pcs%s11 : riv_temp_pcs%e11))
664 : else
665 26280 : call mRM_updateDataset(run_cfg%nc_mrm, L11_qmod(run_cfg%s11 : run_cfg%e11))
666 : end if
667 :
668 : ! write data
669 35040 : if (run_cfg%domainDateTime%writeout(timeStep_model_outputs_mrm, tt)) then
670 732 : call run_cfg%nc_mrm%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
671 : end if
672 :
673 35040 : if(tt == run_cfg%domainDateTime%nTimeSteps) then
674 4 : call run_cfg%nc_mrm%close()
675 : end if
676 :
677 35040 : if ( gw_coupling ) then
678 : ! create
679 0 : if (run_cfg%domainDateTime%tIndex_out == 1) then
680 0 : run_cfg%nc_gw = GW_OutputDataset(iDomain, level0(iDomain)%mask)
681 : end if
682 : ! add data
683 0 : s0 = level0(iDomain)%iStart
684 0 : e0 = level0(iDomain)%iEnd
685 0 : call GW_updateDataset(run_cfg%nc_gw, L0_river_head_mon_sum(s0 : e0))
686 : ! write
687 0 : if (run_cfg%domainDateTime%writeout(-2, tt)) then ! -2 for monthly
688 0 : call run_cfg%nc_gw%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
689 : end if
690 : ! close
691 0 : if(tt == run_cfg%domainDateTime%nTimeSteps) then
692 0 : call run_cfg%nc_gw%close()
693 : end if
694 : end if
695 :
696 : end if
697 :
698 587280 : if (any(outputFlxState)) then
699 :
700 131520 : if (run_cfg%domainDateTime%tIndex_out == 1) then
701 38 : run_cfg%nc_mhm = mHM_OutputDataset(iDomain, run_cfg%mask1)
702 : end if
703 :
704 : call mHM_updateDataset(&
705 : run_cfg%nc_mhm, &
706 : run_cfg%s1, run_cfg%e1, &
707 : L1_fSealed(:, 1, run_cfg%domainDateTime%yId), &
708 : run_cfg%L1_fNotSealed(:, 1, run_cfg%domainDateTime%yId), &
709 : L1_inter, &
710 : L1_snowPack, &
711 : L1_soilMoist, &
712 : L1_soilMoistSat(:, :, run_cfg%domainDateTime%yId), &
713 : L1_sealSTW, &
714 : L1_unsatSTW, &
715 : L1_satSTW, &
716 : L1_neutrons, &
717 : L1_pet_calc, &
718 : L1_aETSoil, &
719 : L1_aETCanopy, &
720 : L1_aETSealed, &
721 : L1_total_runoff, &
722 : L1_runoffSeal, &
723 : L1_fastRunoff, &
724 : L1_slowRunoff, &
725 : L1_baseflow, &
726 : L1_percol, &
727 : L1_infilSoil, &
728 : L1_preEffect, &
729 : L1_melt &
730 131520 : )
731 :
732 : ! write data
733 131520 : if (run_cfg%domainDateTime%writeout(timeStep_model_outputs, tt)) then
734 93 : call run_cfg%nc_mhm%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
735 : end if
736 :
737 131520 : if(tt == run_cfg%domainDateTime%nTimeSteps) then
738 15 : call run_cfg%nc_mhm%close()
739 : end if
740 :
741 : end if
742 : end if ! <-- if (.not. optimize)
743 :
744 2020272 : end subroutine mhm_interface_run_write_output
745 :
746 : !> \brief add simulation data to optimization data types
747 4040544 : subroutine mhm_interface_run_update_optisim(etOptiSim, twsOptiSim, neutronsOptiSim, smOptiSim)
748 : implicit none
749 : !> returns soil moisture time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
750 : type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim
751 : !> dim1=ncells, dim2=time
752 : type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim
753 : !> returns evapotranspiration time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
754 : type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim
755 : !> returns tws time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
756 : type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim
757 :
758 : integer(i4) :: gg, s1, e1, lcId
759 :
760 : integer(i4) :: iDomain, tt
761 :
762 : ! get time step
763 2020272 : tt = run_cfg%time_step
764 :
765 : ! get domain index
766 2020272 : iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
767 :
768 : ! slice settings
769 2020272 : s1 = run_cfg%s1
770 2020272 : e1 = run_cfg%e1
771 2020272 : lcId = run_cfg%domainDateTime%yId
772 :
773 : !----------------------------------------------------------------------
774 : ! FOR SOIL MOISTURE
775 : ! NOTE:: modeled soil moisture is averaged according to input time step
776 : ! soil moisture (timeStep_sm_input)
777 : !----------------------------------------------------------------------
778 2020272 : if (present(smOptiSim)) then
779 : ! only for evaluation period - ignore warming days
780 78624 : if ((tt - warmingDays(iDomain) * nTstepDay) .GT. 0) then
781 : ! decide for daily, monthly or yearly aggregation
782 52560 : call smOptiSim(iDomain)%average_per_timestep(L1_smObs(iDomain)%timeStepInput, &
783 105120 : run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
784 : ! last timestep is already done - write_counter exceeds size(smOptiSim(iDomain)%dataSim, dim=2)
785 52560 : if (tt /= run_cfg%domainDateTime%nTimeSteps) then
786 : ! aggregate soil moisture to needed time step for optimization
787 52554 : call smOptiSim(iDomain)%average_add(sum(L1_soilMoist(:, 1 : nSoilHorizons_sm_input), dim = 2) / &
788 5465616 : sum(L1_soilMoistSat(:, 1 : nSoilHorizons_sm_input, lcId), dim = 2))
789 : end if
790 : end if
791 : end if
792 :
793 : !----------------------------------------------------------------------
794 : ! FOR NEUTRONS
795 : ! NOTE:: modeled neutrons are averaged daily
796 : !----------------------------------------------------------------------
797 2020272 : if (present(neutronsOptiSim)) then
798 : ! only for evaluation period - ignore warming days
799 78624 : if ((tt - warmingDays(iDomain) * nTstepDay) .GT. 0) then
800 : ! decide for daily, monthly or yearly aggregation
801 : ! daily
802 52560 : if (run_cfg%domainDateTime%is_new_day) then
803 2190 : call neutronsOptiSim(iDomain)%average()
804 : end if
805 :
806 : ! last timestep is already done - write_counter exceeds size(sm_opti, dim=2)
807 52560 : if (tt /= run_cfg%domainDateTime%nTimeSteps) then
808 : ! aggregate neutrons to needed time step for optimization
809 52554 : call neutronsOptiSim(iDomain)%average_add(L1_neutrons(s1 : e1))
810 : end if
811 : end if
812 : end if
813 :
814 : !----------------------------------------------------------------------
815 : ! FOR EVAPOTRANSPIRATION
816 : ! NOTE:: modeled evapotranspiration is averaged according to input time step
817 : ! evapotranspiration (timeStep_et_input)
818 : !----------------------------------------------------------------------
819 2020272 : if (present(etOptiSim)) then
820 : ! only for evaluation period - ignore warming days
821 736344 : if ((tt - warmingDays(iDomain) * nTstepDay) .GT. 0) then
822 : ! decide for daily, monthly or yearly aggregation
823 736344 : call etOptiSim(iDomain)%increment_counter(L1_etObs(iDomain)%timeStepInput, &
824 1472688 : run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
825 :
826 : ! last timestep is already done - write_counter exceeds size(etOptiSim(iDomain)%dataSim, dim=2)
827 736344 : if (tt /= run_cfg%domainDateTime%nTimeSteps) then
828 : ! aggregate evapotranspiration to needed time step for optimization
829 0 : call etOptiSim(iDomain)%add(sum(L1_aETSoil(s1 : e1, :), dim = 2) * &
830 0 : run_cfg%L1_fNotSealed(s1 : e1, 1, lcId) + &
831 0 : L1_aETCanopy(s1 : e1) + &
832 75841269 : L1_aETSealed(s1 : e1) * L1_fSealed(s1 : e1, 1, lcId))
833 : end if
834 : end if
835 : end if
836 :
837 : !----------------------------------------------------------------------
838 : ! FOR TWS
839 : ! NOTE:: modeled tws is averaged according to input time step
840 : ! (timeStepInput)
841 : !----------------------------------------------------------------------
842 2020272 : if (present(twsOptiSim)) then
843 : ! only for evaluation period - ignore warming days
844 841464 : if ((tt - warmingDays(iDomain) * nTstepDay) > 0) then
845 : ! decide for daily, monthly or yearly aggregation
846 841464 : call twsOptiSim(iDomain)%average_per_timestep(L1_twsaObs(iDomain)%timeStepInput, &
847 1682928 : run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
848 :
849 : ! last timestep is already done - write_counter exceeds size(twsOptiSim(iDomain)%dataSim, dim=2)
850 841464 : if (tt /= run_cfg%domainDateTime%nTimeSteps) then
851 : ! aggregate evapotranspiration to needed time step for optimization
852 0 : call twsOptiSim(iDomain)%average_add(L1_inter(s1 : e1) + L1_snowPack(s1 : e1) + L1_sealSTW(s1 : e1) + &
853 29450295 : L1_unsatSTW(s1 : e1) + L1_satSTW(s1 : e1))
854 2524311 : do gg = 1, nSoilHorizons_mHM
855 2524311 : call twsOptiSim(iDomain)%add(L1_soilMoist (s1 : e1, gg))
856 : end do
857 : end if
858 : end if
859 : end if
860 :
861 : ! TODO-RIV-TEMP: add OptiSim for river temperature
862 :
863 2020272 : end subroutine mhm_interface_run_update_optisim
864 :
865 : !> \brief finalize current domain after running
866 82 : subroutine mhm_interface_run_finalize_domain()
867 : implicit none
868 :
869 : integer(i4) :: iDomain
870 :
871 : ! get domain index
872 82 : iDomain = run_cfg%get_domain_index(run_cfg%selected_domain)
873 :
874 46 : if (allocated(run_cfg%InflowDischarge)) deallocate(run_cfg%InflowDischarge)
875 82 : if (domainMeta%doRouting(iDomain)) then
876 : ! clean runoff variable
877 46 : deallocate(run_cfg%RunToRout)
878 46 : if ( riv_temp_pcs%active ) call riv_temp_pcs%dealloc_lateral()
879 : end if
880 :
881 2020272 : end subroutine mhm_interface_run_finalize_domain
882 :
883 : !> \brief finalize run
884 61 : subroutine mhm_interface_run_finalize(runoff, BFI)
885 : implicit none
886 : !> returns runoff time series, DIMENSION [nTimeSteps, nGaugesTotal]
887 : real(dp), dimension(:, :), allocatable, optional, intent(out) :: runoff
888 : !> baseflow index, dim1=domainID
889 : real(dp), dimension(:), allocatable, optional, intent(out) :: BFI
890 :
891 508027 : if (present(runoff) .and. (processMatrix(8, 1) > 0)) runoff = mRM_runoff
892 61 : if (present(BFI)) then
893 0 : BFI = BFI_qBF_sum / BFI_qT_sum
894 0 : deallocate(BFI_qBF_sum)
895 0 : deallocate(BFI_qT_sum)
896 : end if
897 :
898 61 : if (allocated(run_cfg%parameterset)) deallocate(run_cfg%parameterset)
899 61 : if (allocated(run_cfg%L1_fNotSealed)) deallocate(run_cfg%L1_fNotSealed)
900 61 : if (allocated(run_cfg%domain_indices)) deallocate(run_cfg%domain_indices)
901 :
902 82 : end subroutine mhm_interface_run_finalize
903 :
904 : end module mo_mhm_interface_run
|
