Line data Source code
1 : !> \file mo_mrm_objective_function_runoff.f90
2 : !> \brief \copybrief mo_mrm_objective_function_runoff
3 : !> \details \copydetails mo_mrm_objective_function_runoff
4 :
5 : !> \brief Objective Functions for Optimization of mHM/mRM against runoff.
6 : !> \details This module provides a wrapper for several objective functions used to optimize mRM/mHM against
7 : !> runoff.
8 : !!
9 : !> If the objective contains besides runoff another variable like TWS move it to mHM/mo_objective_function.f90.
10 : !> If it is only regarding runoff implement it here.
11 : !!
12 : !! All the objective functions are supposed to be minimized!
13 : !! 1. SO: Q: 1.0 - NSE
14 : !! 2. SO: Q: 1.0 - lnNSE
15 : !! 3. SO: Q: 1.0 - 0.5*(NSE+lnNSE)
16 : !! 4. SO: Q: -1.0 * loglikelihood with trend removed from absolute errors and then lag(1)-autocorrelation removed
17 : !! 5. SO: Q: ((1-NSE)**6+(1-lnNSE)**6)**(1/6)
18 : !! 6. SO: Q: SSE
19 : !! 7. SO: Q: -1.0 * loglikelihood with trend removed from absolute errors
20 : !! 8. SO: Q: -1.0 * loglikelihood with trend removed from the relative errors and then lag(1)-autocorrelation removed
21 : !! 9. SO: Q: 1.0 - KGE (Kling-Gupta efficiency measure)
22 : !! 14. SO: Q: sum[((1.0-KGE_i)/ nGauges)**6]**(1/6) > combination of KGE of every gauging station based on a power-6 norm
23 : !! 16. MO: Q: 1st objective: 1.0 - NSE
24 : !! Q: 2nd objective: 1.0 - lnNSE
25 : !! 18. MO: Q: 1st objective: 1.0 - lnNSE(Q_highflow) (95% percentile)
26 : !! Q: 2nd objective: 1.0 - lnNSE(Q_lowflow) (5% of data range)
27 : !! 19. MO: Q: 1st objective: 1.0 - lnNSE(Q_highflow) (non-low flow)
28 : !! Q: 2nd objective: 1.0 - lnNSE(Q_lowflow) (5% of data range)eshold for Q
29 : !! 20. MO: Q: 1st objective: absolute difference in FDC's low-segment volume
30 : !! Q: 2nd objective: 1.0 - NSE of discharge of months DJF
31 : !! 31. SO: Q: 1.0 - wNSE - weighted NSE
32 : !! 32. SO: Q: SSE of boxcox-transformed streamflow
33 : !> \changelog
34 : !! - Stephan Thober Oct 2015
35 : !! - adapted for mRM
36 : !! - Juliane Mai Nov 2015
37 : !! - introducing multi
38 : !! - and single-objective
39 : !! - first multi-objective function (16), but not used yet
40 : !! - Juliane Mai Feb 2016
41 : !! - multi-objective function (18) using lnNSE(highflows) and lnNSE(lowflows)
42 : !! - multi-objective function (19) using lnNSE(highflows) and lnNSE(lowflows)
43 : !! - multi-objective function (20) using FDC and discharge of months DJF
44 : !! - Stephan Thober,Bjoern Guse May 2018
45 : !! - single objective function (21) using weighted NSE following (Hundecha and Bardossy, 2004)
46 : !! - Robert Schweppe Jun 2018
47 : !! - refactoring and reformatting
48 : !! - Stephan Thober Aug 2019
49 : !! - added OF 32: SSE of boxcox-transformed streamflow
50 : !> \authors Juliane Mai
51 : !> \date Dec 2012
52 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
53 : !! mHM is released under the LGPLv3+ license \license_note
54 : !> \ingroup f_mrm
55 : MODULE mo_mrm_objective_function_runoff
56 :
57 : ! This module provides objective functions for optimization of the UFZ CHS mesoscale hydrologic model mHM.
58 :
59 : ! Written Juliane Mai, Dec 2012
60 : ! Modified Stephan Thober, Oct 2015 - removed all none runoff objectives,
61 : ! these can be found mo_objective_functions_sm
62 :
63 : USE mo_kind, ONLY : i4, dp
64 : use mo_optimization_utils, only : eval_interface
65 : use mo_message, only: message, error_message
66 : use mo_string_utils, only : num2str
67 :
68 : IMPLICIT NONE
69 :
70 : PRIVATE
71 :
72 : PUBLIC :: single_objective_runoff ! single-objective function wrapper
73 : PUBLIC :: multi_objective_runoff ! multi-objective function wrapper
74 : PUBLIC :: extract_runoff ! extract runoff period specified in mhm.nml from available runoff time series
75 : #ifdef MPI
76 : PUBLIC :: single_objective_runoff_master, single_objective_runoff_subprocess ! objective function wrapper for soil moisture only
77 : #endif
78 :
79 : ! ------------------------------------------------------------------
80 :
81 : CONTAINS
82 :
83 : ! ------------------------------------------------------------------
84 :
85 : ! NAME
86 : ! single_objective_runoff
87 :
88 : ! PURPOSE
89 : !> \brief Wrapper for objective functions optimizing agains runoff.
90 :
91 : !> \details The functions selects the objective function case defined in a namelist,
92 : !> i.e. the global variable \e opti\_function.
93 : !> It return the objective function value for a specific parameter set.
94 :
95 : ! INTENT(IN)
96 : !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
97 : !> \param[in] "procedure(eval_interface) :: eval"
98 :
99 : ! INTENT(IN), OPTIONAL
100 : !> \param[in] "real(dp), optional :: arg1"
101 :
102 : ! INTENT(OUT), OPTIONAL
103 : !> \param[out] "real(dp), optional :: arg2"
104 : !> \param[out] "real(dp), optional :: arg3"
105 :
106 : ! RETURN
107 : !> \return real(dp) :: objective — objective function value
108 : !> (which will be e.g. minimized by an optimization routine like DDS)
109 :
110 : ! HISTORY
111 : !> \authors Juliane Mai
112 :
113 : !> \date Dec 2012
114 :
115 : ! Modifications:
116 : ! Stephan Thober Oct 2015 - only runoff objective functions
117 : ! Stephan Thober, Bjoern Guse May 2018 - added weighted objective function
118 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
119 :
120 :
121 6 : FUNCTION single_objective_runoff(parameterset, eval, arg1, arg2, arg3)
122 :
123 : use mo_common_mHM_mRM_variables, only : opti_function, opti_method
124 :
125 : implicit none
126 :
127 : REAL(dp), DIMENSION(:), INTENT(IN) :: parameterset
128 :
129 : procedure(eval_interface), INTENT(IN), POINTER :: eval
130 :
131 : real(dp), optional, intent(in) :: arg1
132 :
133 : real(dp), optional, intent(out) :: arg2
134 :
135 : real(dp), optional, intent(out) :: arg3
136 :
137 : REAL(dp) :: single_objective_runoff
138 :
139 :
140 : !write(*,*) 'parameterset: ',parameterset(:)
141 0 : select case (opti_function)
142 : case (1)
143 : ! 1.0-nse
144 0 : single_objective_runoff = objective_nse(parameterset, eval)
145 : case (2)
146 : ! 1.0-lnnse
147 0 : single_objective_runoff = objective_lnnse(parameterset, eval)
148 : case (3)
149 : ! 1.0-0.5*(nse+lnnse)
150 6 : single_objective_runoff = objective_equal_nse_lnnse(parameterset, eval)
151 : case (4)
152 0 : if (opti_method .eq. 0_i4) then
153 : ! MCMC
154 0 : single_objective_runoff = loglikelihood_stddev(parameterset, eval, arg1, arg2, arg3)
155 : else
156 : ! -loglikelihood with trend removed from absolute errors and then lag(1)-autocorrelation removed
157 0 : single_objective_runoff = - loglikelihood_stddev(parameterset, eval, 1.0_dp)
158 : end if
159 : case (5)
160 : ! ((1-NSE)**6+(1-lnNSE)**6)**(1/6)
161 0 : single_objective_runoff = objective_power6_nse_lnnse(parameterset, eval)
162 : case (6)
163 : ! SSE
164 0 : single_objective_runoff = objective_sse(parameterset, eval)
165 : case (7)
166 : ! -loglikelihood with trend removed from absolute errors
167 0 : single_objective_runoff = -loglikelihood_trend_no_autocorr(parameterset, eval, 1.0_dp)
168 : case (8)
169 0 : if (opti_method .eq. 0_i4) then
170 : ! MCMC
171 0 : single_objective_runoff = loglikelihood_evin2013_2(parameterset, eval, regularize = .true.)
172 : else
173 : ! -loglikelihood of approach 2 of Evin et al. (2013),
174 : ! i.e. linear error model with lag(1)-autocorrelation on relative errors
175 0 : single_objective_runoff = -loglikelihood_evin2013_2(parameterset, eval)
176 : end if
177 :
178 : case (9)
179 : ! KGE
180 0 : single_objective_runoff = objective_kge(parameterset, eval)
181 : case (14)
182 : ! combination of KGE of every gauging station based on a power-6 norm \n
183 : ! sum[((1.0-KGE_i)/ nGauges)**6]**(1/6)
184 0 : single_objective_runoff = objective_multiple_gauges_kge_power6(parameterset, eval)
185 : case (31)
186 : ! weighted NSE with observed streamflow
187 0 : single_objective_runoff = objective_weighted_nse(parameterset, eval)
188 : case (32)
189 : ! sum of squared errors (SSE) of boxcox_transformed streamflow
190 0 : single_objective_runoff = objective_sse_boxcox(parameterset, eval)
191 : case default
192 6 : call error_message("Error objective: This opti_function is either not implemented yet or is not a single-objective one.")
193 : end select
194 :
195 6 : END FUNCTION single_objective_runoff
196 :
197 :
198 : ! ------------------------------------------------------------------
199 :
200 : ! NAME
201 : ! single_objective_runoff_master
202 :
203 : ! PURPOSE
204 : !> \brief Wrapper for objective functions optimizing agains runoff.
205 :
206 : !> \details The functions selects the objective function case defined in a namelist,
207 : !> i.e. the global variable \e opti\_function.
208 : !> It return the objective function value for a specific parameter set.
209 :
210 : ! INTENT(IN)
211 : !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
212 : !> \param[in] "procedure(eval_interface) :: eval"
213 :
214 : ! INTENT(IN), OPTIONAL
215 : !> \param[in] "real(dp), optional :: arg1"
216 :
217 : ! INTENT(OUT), OPTIONAL
218 : !> \param[out] "real(dp), optional :: arg2"
219 : !> \param[out] "real(dp), optional :: arg3"
220 :
221 : ! RETURN
222 : !> \return real(dp) :: objective — objective function value
223 : !> (which will be e.g. minimized by an optimization routine like DDS)
224 :
225 : ! HISTORY
226 : !> \authors Juliane Mai
227 :
228 : !> \date Dec 2012
229 :
230 : ! Modifications:
231 : ! Stephan Thober Oct 2015 - only runoff objective functions
232 : ! Stephan Thober, Bjoern Guse May 2018 - added weighted objective function
233 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
234 :
235 : #ifdef MPI
236 : FUNCTION single_objective_runoff_master(parameterset, eval, arg1, arg2, arg3)
237 :
238 : use mo_common_mHM_mRM_variables, only : opti_function, opti_method
239 : use mo_common_mpi_tools, only : distribute_parameterset
240 : use mo_mrm_global_variables, only: nGaugesTotal
241 : use mo_common_variables, only : domainMeta
242 : use mpi_f08
243 :
244 : implicit none
245 :
246 : REAL(dp), DIMENSION(:), INTENT(IN) :: parameterset
247 :
248 : procedure(eval_interface), INTENT(IN), POINTER :: eval
249 :
250 : real(dp), optional, intent(in) :: arg1
251 :
252 : real(dp), optional, intent(out) :: arg2
253 :
254 : real(dp), optional, intent(out) :: arg3
255 :
256 : REAL(dp) :: single_objective_runoff_master
257 :
258 : REAL(dp) :: partial_objective
259 :
260 : ! for sixth root
261 : real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
262 :
263 : integer(i4) :: iproc, nproc
264 :
265 : integer(i4) :: ierror
266 :
267 : type(MPI_Status) :: status
268 :
269 : !write(*,*) 'parameterset: ',parameterset(:)
270 : call distribute_parameterset(parameterset)
271 : select case (opti_function)
272 : case(1 : 3, 5, 6, 9, 31)
273 : call MPI_Comm_size(domainMeta%comMaster, nproc, ierror)
274 : single_objective_runoff_master = 0.0_dp
275 : do iproc = 1, nproc - 1
276 : call MPI_Recv(partial_objective, 1, MPI_DOUBLE_PRECISION, iproc, 0, domainMeta%comMaster, status, ierror)
277 : single_objective_runoff_master = single_objective_runoff_master + partial_objective
278 : end do
279 : single_objective_runoff_master = 1.0_dp - single_objective_runoff_master / real(nGaugesTotal, dp)
280 : case(14)
281 : call MPI_Comm_size(domainMeta%comMaster, nproc, ierror)
282 : single_objective_runoff_master = 0.0_dp
283 : do iproc = 1, nproc - 1
284 : call MPI_Recv(partial_objective, 1, MPI_DOUBLE_PRECISION, iproc, 0, domainMeta%comMaster, status, ierror)
285 : single_objective_runoff_master = single_objective_runoff_master + partial_objective
286 : end do
287 : single_objective_runoff_master = single_objective_runoff_master**onesixth
288 : case(4, 7, 8)
289 : call message("case 4, 7, 8 are not implemented in parallel yet")
290 : case default
291 : call error_message("Error single_objective_runoff_master:", raise=.false.)
292 : call error_message("This opti_function is either not implemented yet or is not a single-objective one.")
293 : end select
294 :
295 : select case (opti_function)
296 : case (1)
297 : ! 1.0-nse
298 : call message('objective_nse (i.e., 1 - NSE) = ', num2str(single_objective_runoff_master))
299 : case (2)
300 : ! 1.0-lnnse
301 : call message('objective_lnnse = ', num2str(single_objective_runoff_master))
302 : case (3)
303 : ! 1.0-0.5*(nse+lnnse)
304 : call message('objective_equal_nse_lnnse = ', num2str(single_objective_runoff_master))
305 : case (4)
306 : call error_message("case 4, loglikelihood_stddev not implemented in parallel yet")
307 : case (5)
308 : ! ((1-NSE)**6+(1-lnNSE)**6)**(1/6)
309 : call message('objective_power6_nse_lnnse = ', num2str(single_objective_runoff_master))
310 : case (6)
311 : ! SSE
312 : call message('objective_sse = ', num2str(single_objective_runoff_master))
313 : case (7)
314 : ! -loglikelihood with trend removed from absolute errors
315 : call error_message("case 7, single_objective_runoff_master not implemented in parallel yet")
316 : case (8)
317 : call error_message("case 8, loglikelihood_evin2013_2 not implemented in parallel yet")
318 : case (9)
319 : ! KGE
320 : call message('objective_kge (i.e., 1 - KGE) = ', num2str(single_objective_runoff_master))
321 : case (14)
322 : ! combination of KGE of every gauging station based on a power-6 norm \n
323 : ! sum[((1.0-KGE_i)/ nGauges)**6]**(1/6)
324 : call message('objective_multiple_gauges_kge_power6 = ', num2str(single_objective_runoff_master))
325 : case (31)
326 : ! weighted NSE with observed streamflow
327 : call message('objective_weighted_nse (i.e., 1 - wNSE) = ', num2str(single_objective_runoff_master))
328 : case (32)
329 : ! SSE of boxcox-transformed streamflow
330 : call message('sse_boxcox_streamflow = ', num2str(single_objective_runoff_master))
331 : case default
332 : call error_message("Error single_objective_runoff_master:", raise=.false.)
333 : call error_message("This opti_function is either not implemented yet or is not a single-objective one.", raise=.false.)
334 : call error_message("This part of the code should never be executed.")
335 : end select
336 :
337 : END FUNCTION single_objective_runoff_master
338 :
339 :
340 : ! ------------------------------------------------------------------
341 :
342 : ! NAME
343 : ! single_objective_runoff_subprocess
344 :
345 : ! PURPOSE
346 : !> \brief Wrapper for objective functions optimizing agains runoff.
347 :
348 : !> \details The functions selects the objective function case defined in a namelist,
349 : !> i.e. the global variable \e opti\_function.
350 : !> It return the objective function value for a specific parameter set.
351 :
352 : ! INTENT(IN)
353 : !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
354 : !> \param[in] "procedure(eval_interface) :: eval"
355 :
356 : ! INTENT(IN), OPTIONAL
357 : !> \param[in] "real(dp), optional :: arg1"
358 :
359 : ! INTENT(OUT), OPTIONAL
360 : !> \param[out] "real(dp), optional :: arg2"
361 : !> \param[out] "real(dp), optional :: arg3"
362 :
363 : ! RETURN
364 : !> \return real(dp) :: objective — objective function value
365 : !> (which will be e.g. minimized by an optimization routine like DDS)
366 :
367 : ! HISTORY
368 : !> \authors Juliane Mai
369 :
370 : !> \date Dec 2012
371 :
372 : ! Modifications:
373 : ! Stephan Thober Oct 2015 - only runoff objective functions
374 : ! Stephan Thober, Bjoern Guse May 2018 - added weighted objective function
375 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
376 :
377 :
378 : subroutine single_objective_runoff_subprocess(eval, arg1, arg2, arg3)
379 :
380 : use mo_common_mHM_mRM_variables, only : opti_function, opti_method
381 : use mo_common_mpi_tools, only : get_parameterset
382 : use mo_common_variables, only : domainMeta
383 : use mpi_f08
384 :
385 : implicit none
386 :
387 : procedure(eval_interface), INTENT(IN), POINTER :: eval
388 :
389 : real(dp), optional, intent(in) :: arg1
390 :
391 : real(dp), optional, intent(out) :: arg2
392 :
393 : real(dp), optional, intent(out) :: arg3
394 :
395 : REAL(dp) :: partial_single_objective_runoff
396 :
397 : REAL(dp), DIMENSION(:), allocatable :: parameterset
398 :
399 : integer(i4) :: ierror
400 :
401 : type(MPI_Status) :: status
402 :
403 : logical :: do_obj_loop
404 :
405 :
406 : do ! a do loop without condition runs until exit
407 : call MPI_Recv(do_obj_loop, 1, MPI_LOGICAL, 0, 0, domainMeta%comMaster, status, ierror)
408 :
409 : if (.not. do_obj_loop) exit
410 : !write(*,*) 'parameterset: ',parameterset(:)
411 : call get_parameterset(parameterset)
412 : select case (opti_function)
413 : case (1)
414 : ! 1.0-nse
415 : partial_single_objective_runoff = objective_nse(parameterset, eval)
416 : case (2)
417 : ! 1.0-lnnse
418 : partial_single_objective_runoff = objective_lnnse(parameterset, eval)
419 : case (3)
420 : ! 1.0-0.5*(nse+lnnse)
421 : partial_single_objective_runoff = objective_equal_nse_lnnse(parameterset, eval)
422 : case (4)
423 : if (opti_method .eq. 0_i4) then
424 : ! MCMC
425 : ! partial_single_objective_runoff = loglikelihood_stddev(parameterset, eval, arg1, arg2, arg3)
426 : call error_message("Error single_objective_runoff_subprocess: case 4 with optimethod 0 not supported")
427 : else
428 : ! -loglikelihood with trend removed from absolute errors and then lag(1)-autocorrelation removed
429 : partial_single_objective_runoff = - loglikelihood_stddev(parameterset, eval, 1.0_dp)
430 : end if
431 : case (5)
432 : ! ((1-NSE)**6+(1-lnNSE)**6)**(1/6)
433 : partial_single_objective_runoff = objective_power6_nse_lnnse(parameterset, eval)
434 : case (6)
435 : ! SSE
436 : partial_single_objective_runoff = objective_sse(parameterset, eval)
437 : case (7)
438 : ! -loglikelihood with trend removed from absolute errors
439 : ! partial_single_objective_runoff = -loglikelihood_trend_no_autocorr(parameterset, eval, 1.0_dp)
440 : call error_message("Error single_objective_runoff_subprocess: case 7 not supported")
441 : case (8)
442 : if (opti_method .eq. 0_i4) then
443 : ! MCMC
444 : ! partial_single_objective_runoff = loglikelihood_evin2013_2(parameterset, eval, regularize = .true.)
445 : call error_message("Error single_objective_runoff_subprocess: case 8 with optimethod 0 not supported")
446 : else
447 : ! -loglikelihood of approach 2 of Evin et al. (2013),
448 : ! i.e. linear error model with lag(1)-autocorrelation on relative errors
449 : partial_single_objective_runoff = -loglikelihood_evin2013_2(parameterset, eval)
450 : end if
451 :
452 : case (9)
453 : ! KGE
454 : partial_single_objective_runoff = objective_kge(parameterset, eval)
455 : case (14)
456 : ! combination of KGE of every gauging station based on a power-6 norm \n
457 : ! sum[((1.0-KGE_i)/ nGauges)**6]**(1/6)
458 : partial_single_objective_runoff = objective_multiple_gauges_kge_power6(parameterset, eval)
459 : case (31)
460 : ! weighted NSE with observed streamflow
461 : partial_single_objective_runoff = objective_weighted_nse(parameterset, eval)
462 : case (32)
463 : ! SSE of transformed streamflow
464 : partial_single_objective_runoff = objective_sse_boxcox(parameterset, eval)
465 : case default
466 : call error_message("Error single_objective_runoff_subprocess:", raise=.false.)
467 : call error_message("This opti_function is either not implemented yet or is not a single-objective one.")
468 : end select
469 :
470 : select case (opti_function)
471 : case (1 : 3, 5, 6, 9, 14, 31)
472 : call MPI_Send(partial_single_objective_runoff,1, MPI_DOUBLE_PRECISION,0,0,domainMeta%comMaster,ierror)
473 : case default
474 : call error_message("Error objective_subprocess: this part should not be executed -> error in the code.")
475 : end select
476 : deallocate(parameterset)
477 : end do
478 :
479 : END subroutine single_objective_runoff_subprocess
480 :
481 : #endif
482 :
483 : ! ------------------------------------------------------------------
484 :
485 : ! NAME
486 : ! multi_objective_runoff
487 :
488 : ! PURPOSE
489 : !> \brief Wrapper for multi-objective functions where at least one is regarding runoff.
490 :
491 : !> \details The functions selects the objective function case defined in a namelist,
492 : !> i.e. the global variable \e opti\_function.
493 : !> It return the multiple objective function values for a specific parameter set.
494 :
495 : ! INTENT(IN)
496 : !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
497 : !> \param[in] "procedure(eval_interface) :: eval"
498 :
499 : ! INTENT(OUT)
500 : !> \param[out] "REAL(dp), DIMENSION(:) :: multi_objectives"
501 :
502 : ! HISTORY
503 : !> \authors Juliane Mai
504 :
505 : !> \date Oct 2015
506 :
507 : ! Modifications:
508 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
509 :
510 0 : SUBROUTINE multi_objective_runoff(parameterset, eval, multi_objectives)
511 :
512 6 : use mo_common_mHM_mRM_variables, only : opti_function
513 :
514 : implicit none
515 :
516 : REAL(dp), DIMENSION(:), INTENT(IN) :: parameterset
517 :
518 : procedure(eval_interface), INTENT(IN), pointer :: eval
519 :
520 : REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: multi_objectives
521 :
522 :
523 0 : select case (opti_function)
524 : case (16)
525 : ! 1st objective: 1.0-nse
526 : ! 2nd objective: 1.0-lnnse
527 0 : multi_objectives = multi_objective_nse_lnnse(parameterset, eval)
528 : case (18)
529 : ! 1st objective: 1.0 - lnNSE(Q_highflow) (95% percentile)
530 : ! 2nd objective: 1.0 - lnNSE(Q_lowflow) (5% of data range)
531 0 : multi_objectives = multi_objective_lnnse_highflow_lnnse_lowflow(parameterset, eval)
532 : case (19)
533 : ! 1st objective: 1.0 - lnNSE(Q_highflow) (non-low flow)
534 : ! 2nd objective: 1.0 - lnNSE(Q_lowflow) (5% of data range)
535 0 : multi_objectives = multi_objective_lnnse_highflow_lnnse_lowflow_2(parameterset, eval)
536 : case (20)
537 : ! 1st objective: 1.0 - difference in FDC's low-segment volume
538 : ! 2nd objective: 1.0 - NSE of discharge of months DJF
539 0 : multi_objectives = multi_objective_ae_fdc_lsv_nse_djf(parameterset, eval)
540 : case default
541 0 : call error_message("Error objective: Either this opti_function is not implemented yet or it is not a multi-objective one.")
542 : end select
543 :
544 0 : END SUBROUTINE multi_objective_runoff
545 :
546 : ! ------------------------------------------------------------------
547 :
548 : ! NAME
549 : ! loglikelihood_stddev
550 :
551 : ! PURPOSE
552 : !> \brief Logarithmic likelihood function with removed linear trend and Lag(1)-autocorrelation.
553 :
554 : !> \details The logarithmis likelihood function is used when mHM runs in MCMC mode.
555 : !> It can also be used for optimization when selecting the likelihood in the
556 : !> namelist as \e opti\_function.
557 :
558 : ! INTENT(IN)
559 : !> \param[in] "real(dp), dimension(:) :: parameterset"
560 : !> \param[in] "procedure(eval_interface) :: eval"
561 : !> \param[in] "real(dp) :: stddev" standard deviation of data
562 :
563 : ! INTENT(OUT), OPTIONAL
564 : !> \param[out] "real(dp), optional :: stddev_new" standard deviation of errors with removed trend and
565 : !> correlationbetween model run using parameter set and observation
566 : !> \param[out] "real(dp), optional :: likeli_new" logarithmic likelihood determined with stddev_new instead of
567 : !> stddev
568 :
569 : ! RETURN
570 : !> \return real(dp) :: loglikelihood_stddev — logarithmic likelihood using given stddev
571 : !> but remove optimal trend and lag(1)-autocorrelation in errors
572 : !> (absolute between running model with parameterset and observation)
573 :
574 : ! HISTORY
575 : !> \authors Juliane Mai
576 :
577 : !> \date Dec 2012
578 :
579 : ! Modifications:
580 : ! Stephan Thober Jan 2015 - introduced extract_runoff
581 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2) to not interfere with mRM
582 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
583 :
584 0 : FUNCTION loglikelihood_stddev(parameterset, eval, stddev, stddev_new, likeli_new)
585 :
586 0 : use mo_append, only : append
587 : use mo_linfit, only : linfit
588 : use mo_moment, only : correlation, mean
589 :
590 : implicit none
591 :
592 : real(dp), dimension(:), intent(in) :: parameterset
593 :
594 : procedure(eval_interface), INTENT(IN), pointer :: eval
595 :
596 : ! standard deviation of data
597 : real(dp), intent(in) :: stddev
598 :
599 : ! standard deviation of errors with removed trend and correlationbetween model run using parameter set and
600 : ! observation
601 : real(dp), intent(out), optional :: stddev_new
602 :
603 : ! logarithmic likelihood determined with stddev_new instead of stddev
604 : real(dp), intent(out), optional :: likeli_new
605 :
606 : real(dp) :: loglikelihood_stddev
607 :
608 : ! modelled runoff for a given parameter set
609 : ! dim1=nTimeSteps, dim2=nGauges
610 0 : real(dp), dimension(:, :), allocatable :: runoff
611 :
612 : ! gauges counter
613 : integer(i4) :: gg
614 :
615 : ! aggregated simulated runoff
616 0 : real(dp), dimension(:), allocatable :: runoff_agg
617 :
618 : ! measured runoff
619 0 : real(dp), dimension(:), allocatable :: runoff_obs
620 :
621 : ! mask for aggregated measured runoff
622 0 : logical, dimension(:), allocatable :: runoff_obs_mask
623 :
624 : integer(i4) :: nmeas
625 :
626 0 : real(dp), dimension(:), allocatable :: errors
627 :
628 0 : real(dp), dimension(:), allocatable :: obs, calc, out
629 :
630 0 : real(dp) :: a, b, c
631 :
632 0 : real(dp) :: stddev_tmp
633 :
634 :
635 0 : call eval(parameterset, runoff = runoff)
636 :
637 : ! extract runoff and append it to obs and calc
638 0 : do gg = 1, size(runoff, dim = 2) ! second dimension equals nGaugesTotal
639 : ! extract runoff
640 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
641 : ! append it to variables
642 0 : call append(obs, runoff_obs)
643 0 : call append(calc, runoff_agg)
644 :
645 : end do
646 : ! ----------------------------------------
647 :
648 0 : nmeas = size(obs, dim = 1)
649 :
650 0 : allocate(out(nmeas), errors(nmeas))
651 0 : errors(:) = abs(calc(:) - obs(:))
652 :
653 : ! remove linear trend of errors - must be model NOT obs
654 : ! out = linfit(obs, errors, a=a, b=b, model2=.False.)
655 : ! errors(:) = errors(:) - (a + obs(:)*b)
656 0 : out = linfit(calc, errors, a = a, b = b, model2 = .False.)
657 0 : errors(:) = errors(:) - (a + calc(:) * b)
658 :
659 : ! remove lag(1)-autocorrelation of errors
660 0 : c = correlation(errors(2 : nmeas), errors(1 : nmeas - 1))
661 0 : errors(1 : nmeas - 1) = errors(2 : nmeas) - c * errors(1 : nmeas - 1)
662 0 : errors(nmeas) = 0.0_dp
663 :
664 : ! you have to take stddev=const because otherwise loglikelihood is always N
665 : ! in MCMC stddev gets updated only when a better likelihood is found.
666 0 : loglikelihood_stddev = sum(errors(:) * errors(:) / stddev**2)
667 0 : loglikelihood_stddev = -0.5_dp * loglikelihood_stddev
668 :
669 0 : call message('-loglikelihood_stddev = ', num2str(-loglikelihood_stddev))
670 :
671 0 : stddev_tmp = sqrt(sum((errors(:) - mean(errors)) * (errors(:) - mean(errors))) / real(nmeas - 1, dp))
672 0 : if (present(stddev_new)) then
673 0 : stddev_new = stddev_tmp
674 : end if
675 0 : if (present(likeli_new)) then
676 0 : likeli_new = sum(errors(:) * errors(:) / stddev_tmp**2)
677 0 : likeli_new = -0.5_dp * likeli_new
678 : end if
679 :
680 0 : deallocate(runoff, runoff_agg, runoff_obs_mask, runoff_obs)
681 0 : deallocate(obs, calc, out, errors)
682 :
683 0 : END FUNCTION loglikelihood_stddev
684 :
685 : ! ------------------------------------------------------------------
686 :
687 : ! NAME
688 : ! loglikelihood_evin2013_2
689 :
690 : ! PURPOSE
691 : !> \brief Logarithmised likelihood with linear error model and lag(1)-autocorrelation
692 : !> of the relative errors.
693 :
694 : !> \details This loglikelihood uses a linear error model and a lag(1)-autocorrelation
695 : !> on the relative errors. This is approach 2 of the paper Evin et al. (WRR, 2013).
696 :
697 : !> This is opti_function = 8.
698 :
699 : !> mHM then adds two extra (local) parameters for the error model in mhm_driver,
700 : !> which get optimised together with the other, global parameters.
701 : !> ADDITIONAL INFORMATION
702 : !> Evin et al., WRR 49, 4518-4524, 2013
703 :
704 : ! INTENT(IN)
705 : !> \param[in] "real(dp), dimension(:) :: parameterset"
706 : !> \param[in] "procedure(eval_interface) :: eval"
707 :
708 : ! INTENT(IN), OPTIONAL
709 : !> \param[in] "logical, optional :: regularize"
710 :
711 : ! RETURN
712 : !> \return real(dp) :: loglikelihood_evin2013_2 — logarithmic likelihood using given stddev
713 : !> but remove optimal trend and lag(1)-autocorrelation in errors
714 : !> (absolute between running model with parameterset and observation)
715 :
716 : ! HISTORY
717 : !> \authors Juliane Mai and Matthias Cuntz
718 :
719 : !> \date Mar 2014
720 :
721 : ! Modifications:
722 : ! Stephan Thober Jan 2015 - introduced extract_runoff
723 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2) to not interfere with mRM
724 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
725 :
726 0 : FUNCTION loglikelihood_evin2013_2(parameterset, eval, regularize)
727 :
728 0 : use mo_append, only : append
729 : use mo_common_variables, only : global_parameters
730 : use mo_constants, only : pi_dp
731 : use mo_moment, only : correlation
732 : use mo_utils, only : eq
733 :
734 : implicit none
735 :
736 : real(dp), dimension(:), intent(in) :: parameterset
737 :
738 : procedure(eval_interface), INTENT(IN), pointer :: eval
739 :
740 : logical, optional, intent(in) :: regularize
741 :
742 : real(dp) :: loglikelihood_evin2013_2
743 :
744 : ! modelled runoff for a given parameter set
745 : ! dim1=nTimeSteps, dim2=nGauges
746 0 : real(dp), dimension(:, :), allocatable :: runoff
747 :
748 : ! penalty term due to a parmeter set out of bound
749 0 : real(dp) :: penalty
750 :
751 : ! gauges counter
752 : integer(i4) :: gg
753 :
754 : ! aggregated simulated runoff
755 0 : real(dp), dimension(:), allocatable :: runoff_agg
756 :
757 : ! measured runoff
758 0 : real(dp), dimension(:), allocatable :: runoff_obs
759 :
760 : ! mask for measured runoff
761 0 : logical, dimension(:), allocatable :: runoff_obs_mask
762 :
763 : integer(i4) :: nmeas
764 :
765 0 : real(dp), dimension(:), allocatable :: errors, sigma, eta, y
766 :
767 0 : real(dp), dimension(:), allocatable :: obs, calc, out
768 :
769 0 : real(dp) :: a, b, c, vary, vary1, ln2pi, tmp
770 :
771 : integer(i4) :: npara
772 :
773 : logical :: iregularize
774 :
775 :
776 0 : iregularize = .false.
777 0 : if (present(regularize)) iregularize = regularize
778 :
779 0 : npara = size(parameterset)
780 0 : call eval(parameterset(1 : npara - 2), runoff = runoff)
781 :
782 : ! extract runoff and append it to obs and calc
783 0 : do gg = 1, size(runoff, dim = 2) ! second dimension equals nGaugesTotal
784 : ! extract runoff
785 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
786 : ! append it to variables
787 0 : call append(obs, runoff_obs)
788 0 : call append(calc, runoff_agg)
789 :
790 : end do
791 :
792 : ! ----------------------------------------
793 :
794 0 : nmeas = size(obs, dim = 1)
795 :
796 0 : allocate(out(nmeas), errors(nmeas), sigma(nmeas), eta(nmeas), y(nmeas))
797 : ! residual errors
798 0 : errors(:) = calc(:) - obs(:)
799 :
800 : ! linear error model
801 0 : a = parameterset(npara - 1)
802 0 : b = parameterset(npara)
803 0 : sigma(:) = a + b * calc(:)
804 : ! standardized residual errors (SRE)
805 0 : eta(:) = errors(:) / sigma(:)
806 :
807 : ! remove lag(1)-autocorrelation of SRE
808 0 : c = correlation(eta(2 : nmeas), eta(1 : nmeas - 1))
809 0 : y(1) = 0.0_dp ! only for completeness
810 0 : y(2 : nmeas) = eta(2 : nmeas) - c * eta(1 : nmeas - 1)
811 :
812 : ! likelihood of residual errors
813 0 : ln2pi = log(sqrt(2.0_dp * pi_dp))
814 0 : vary = 1.0_dp - c * c
815 0 : vary1 = 1.0_dp / vary
816 : loglikelihood_evin2013_2 = -ln2pi - 0.5_dp * eta(1) * eta(1) - log(sigma(1)) & ! li(eta(1))/sigma(1)
817 : - real(nmeas - 1, dp) * log(sqrt(2.0_dp * pi_dp * vary)) &
818 0 : - sum(0.5_dp * y(2 : nmeas) * y(2 : nmeas) * vary1) - sum(log(sigma(2 : nmeas)))
819 :
820 0 : if (iregularize) then
821 : ! Regularistion term as deviation from initial parameter value
822 : penalty = parameter_regularization(&
823 : parameterset(1 : npara - 2), & ! current parameter set
824 0 : global_parameters(1 : npara - 2, 3), & ! prior/initial parameter set
825 0 : global_parameters(1 : npara - 2, 1 : 2), & ! bounds
826 0 : eq(global_parameters(1 : npara - 2, 4), 1.0_dp)) ! used/unused
827 :
828 0 : tmp = loglikelihood_evin2013_2 + penalty
829 : call message( &
830 : '-loglikelihood_evin2013_2, + penalty, chi^2: ', &
831 : num2str(-loglikelihood_evin2013_2), &
832 : num2str(-tmp), &
833 0 : num2str(-tmp / real(nmeas, dp)))
834 0 : loglikelihood_evin2013_2 = tmp
835 : else
836 : call message( &
837 : '-loglikelihood_evin2013_2, chi^2: ', &
838 : num2str(-loglikelihood_evin2013_2), &
839 0 : num2str(-loglikelihood_evin2013_2 / real(nmeas, dp)))
840 : end if
841 :
842 0 : deallocate(runoff, runoff_agg, runoff_obs_mask, runoff_obs)
843 0 : deallocate(obs, calc, out, errors, sigma, eta, y)
844 :
845 0 : END FUNCTION loglikelihood_evin2013_2
846 :
847 : ! NAME
848 : ! parameter_regularization
849 :
850 : ! PURPOSE
851 : !> \brief TODO: add description
852 :
853 : !> \details TODO: add description
854 :
855 : ! INTENT(IN)
856 : !> \param[in] "real(dp), dimension(:) :: paraset"
857 : !> \param[in] "real(dp), dimension(size(paraset)) :: prior"
858 : !> \param[in] "real(dp), dimension(size(paraset), 2) :: bounds" (min, max)
859 : !> \param[in] "logical, dimension(size(paraset)) :: mask"
860 :
861 : ! HISTORY
862 : !> \authors Robert Schweppe
863 :
864 : !> \date Jun 2018
865 :
866 : ! Modifications:
867 :
868 0 : FUNCTION parameter_regularization(paraset, prior, bounds, mask)
869 :
870 0 : use mo_constants, only : pi_dp
871 :
872 : implicit none
873 :
874 : real(dp), dimension(:), intent(in) :: paraset
875 :
876 : real(dp), dimension(size(paraset)), intent(in) :: prior
877 :
878 : ! (min, max)
879 : real(dp), dimension(size(paraset), 2), intent(in) :: bounds
880 :
881 : logical, dimension(size(paraset)), intent(in) :: mask
882 :
883 : real(dp) :: parameter_regularization
884 :
885 : integer(i4) :: ipara
886 :
887 : integer(i4) :: npara
888 :
889 : real(dp), parameter :: onetwelveth = 1._dp / 12._dp
890 :
891 0 : real(dp), dimension(size(paraset)) :: sigma
892 :
893 :
894 0 : npara = size(paraset, 1)
895 :
896 0 : sigma = sqrt(onetwelveth * (bounds(:, 2) - bounds(:, 1))**2) ! standard deviation of uniform distribution
897 0 : parameter_regularization = -sum(log(sqrt(2.0_dp * pi_dp) * sigma), mask = mask)
898 :
899 0 : do ipara = 1, npara
900 0 : if (mask(ipara)) then
901 : ! if ((paraset(ipara) .lt. bounds(ipara,1)) .or. (paraset(ipara) .gt. bounds(ipara,2))) then
902 : ! ! outside bounds
903 : parameter_regularization = parameter_regularization - &
904 0 : 0.5_dp * ((paraset(ipara) - prior(ipara)) / sigma(ipara))**2
905 : ! else
906 : ! ! in bound
907 : ! parameter_regularization = 0.0_dp
908 : ! end if
909 : end if
910 : end do
911 :
912 0 : END FUNCTION parameter_regularization
913 :
914 : ! ------------------------------------------------------------------
915 :
916 : ! NAME
917 : ! loglikelihood_trend_no_autocorr
918 :
919 : ! PURPOSE
920 : !> \brief Logarithmic likelihood function with linear trend removed.
921 :
922 : !> \details The logarithmis likelihood function is used when mHM runs in MCMC mode.
923 : !> It can also be used for optimization when selecting the likelihood in the
924 : !> namelist as \e opti\_function.
925 :
926 : ! INTENT(IN)
927 : !> \param[in] "real(dp), dimension(:) :: parameterset"
928 : !> \param[in] "procedure(eval_interface) :: eval"
929 : !> \param[in] "real(dp) :: stddev_old" standard deviation of data
930 :
931 : ! INTENT(OUT), OPTIONAL
932 : !> \param[out] "real(dp), optional :: stddev_new" standard deviation of errors with removed trendbetween model
933 : !> run using parameter set and observation
934 : !> \param[out] "real(dp), optional :: likeli_new" logarithmic likelihood determined with stddev_new instead of
935 : !> stddev
936 :
937 : ! RETURN
938 : !> \return real(dp) :: loglikelihood_trend_no_autocorr — logarithmic likelihood using given stddev
939 : !> but remove optimal trend in errors
940 : !> (absolute between running model with parameterset and observation)
941 :
942 : ! HISTORY
943 : !> \authors Juliane Mai and Matthias Cuntz
944 :
945 : !> \date Mar 2014
946 :
947 : ! Modifications:
948 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2)
949 : ! to not interfere with mRM
950 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
951 :
952 0 : FUNCTION loglikelihood_trend_no_autocorr(parameterset, eval, stddev_old, stddev_new, likeli_new)
953 :
954 0 : use mo_append, only : append
955 : use mo_linfit, only : linfit
956 : use mo_moment, only : stddev
957 :
958 : implicit none
959 :
960 : real(dp), dimension(:), intent(in) :: parameterset
961 :
962 : procedure(eval_interface), INTENT(IN), pointer :: eval
963 :
964 : ! standard deviation of data
965 : real(dp), intent(in) :: stddev_old
966 :
967 : ! standard deviation of errors with removed trendbetween model run using parameter set and observation
968 : real(dp), intent(out), optional :: stddev_new
969 :
970 : ! logarithmic likelihood determined with stddev_new instead of stddev
971 : real(dp), intent(out), optional :: likeli_new
972 :
973 : real(dp) :: loglikelihood_trend_no_autocorr
974 :
975 : ! modelled runoff for a given parameter set
976 : ! dim1=nTimeSteps, dim2=nGauges
977 0 : real(dp), dimension(:, :), allocatable :: runoff
978 :
979 : ! gauges counter
980 : integer(i4) :: gg
981 :
982 : ! aggregated simulated runoff
983 0 : real(dp), dimension(:), allocatable :: runoff_agg
984 :
985 : ! measured runoff
986 0 : real(dp), dimension(:), allocatable :: runoff_obs
987 :
988 : ! mask for aggregated measured runoff
989 0 : logical, dimension(:), allocatable :: runoff_obs_mask
990 :
991 : integer(i4) :: nmeas
992 :
993 0 : real(dp), dimension(:), allocatable :: errors
994 :
995 0 : real(dp), dimension(:), allocatable :: obs, calc, out
996 :
997 0 : real(dp) :: a, b
998 :
999 0 : real(dp) :: stddev_tmp
1000 :
1001 :
1002 0 : call eval(parameterset, runoff = runoff)
1003 :
1004 : ! extract runoff and append it to obs and calc
1005 0 : do gg = 1, size(runoff, dim = 2) ! second dimension equals nGaugesTotal
1006 : ! extract runoff
1007 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1008 : ! append it to variables
1009 0 : call append(obs, runoff_obs)
1010 0 : call append(calc, runoff_agg)
1011 :
1012 : end do
1013 :
1014 : ! ----------------------------------------
1015 0 : nmeas = size(obs, dim = 1)
1016 :
1017 : ! allocate output variables
1018 0 : allocate(out(nmeas), errors(nmeas))
1019 0 : errors(:) = abs(calc(:) - obs(:))
1020 :
1021 : ! remove linear trend of errors - must be model NOT obs
1022 0 : out = linfit(calc, errors, a = a, b = b, model2 = .False.)
1023 0 : errors(:) = errors(:) - (a + calc(:) * b)
1024 :
1025 : ! you have to take stddev_old=const because otherwise loglikelihood_trend_no_autocorr is always N
1026 : ! in MCMC stddev_old gets updated only when a better likelihood is found.
1027 0 : loglikelihood_trend_no_autocorr = sum(errors(:) * errors(:) / stddev_old**2)
1028 0 : loglikelihood_trend_no_autocorr = -0.5_dp * loglikelihood_trend_no_autocorr
1029 :
1030 0 : call message('-loglikelihood_trend_no_autocorr = ', num2str(-loglikelihood_trend_no_autocorr))
1031 :
1032 0 : stddev_tmp = 1.0_dp ! initialization
1033 0 : if (present(stddev_new) .or. present(likeli_new)) then
1034 0 : stddev_tmp = stddev(errors(:))
1035 : end if
1036 0 : if (present(stddev_new)) then
1037 0 : stddev_new = stddev_tmp
1038 : end if
1039 0 : if (present(likeli_new)) then
1040 0 : likeli_new = sum(errors(:) * errors(:) / stddev_tmp**2)
1041 0 : likeli_new = -0.5_dp * likeli_new
1042 : end if
1043 :
1044 0 : deallocate(runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1045 0 : deallocate(obs, calc, out, errors)
1046 :
1047 0 : END FUNCTION loglikelihood_trend_no_autocorr
1048 :
1049 : ! ------------------------------------------------------------------
1050 :
1051 : ! NAME
1052 : ! objective_lnnse
1053 :
1054 : ! PURPOSE
1055 : !> \brief Objective function of logarithmic NSE.
1056 :
1057 : !> \details The objective function only depends on a parameter vector.
1058 : !> The model will be called with that parameter vector and
1059 : !> the model output is subsequently compared to observed data.
1060 : !> Therefore, the logarithmic Nash-Sutcliffe model efficiency coefficient \f$ lnNSE \f$
1061 : !> \f[ lnNSE = 1 - \frac{\sum_{i=1}^N (\ln Q_{obs}(i) - \ln Q_{model}(i))^2}
1062 : !> {\sum_{i=1}^N (\ln Q_{obs}(i) - \bar{\ln Q_{obs}})^2} \f]
1063 : !> is calculated.
1064 : !> \f[ obj\_value = lnNSE \f]
1065 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1066 :
1067 : ! INTENT(IN)
1068 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1069 : !> \param[in] "procedure(eval_interface) :: eval"
1070 :
1071 : ! RETURN
1072 : !> \return real(dp) :: objective_lnnse — objective function value
1073 : !> (which will be e.g. minimized by an optimization routine like DDS)
1074 :
1075 : ! HISTORY
1076 : !> \authors Juliane Mai
1077 :
1078 : !> \date May 2013
1079 :
1080 : ! Modifications:
1081 : ! Stephan Thober Jan 2015 - introduced extract_runoff
1082 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2) to not interfere with mRM
1083 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1084 :
1085 0 : FUNCTION objective_lnnse(parameterset, eval)
1086 :
1087 0 : use mo_errormeasures, only : lnnse
1088 :
1089 : implicit none
1090 :
1091 : real(dp), dimension(:), intent(in) :: parameterset
1092 :
1093 : procedure(eval_interface), INTENT(IN), pointer :: eval
1094 :
1095 : real(dp) :: objective_lnnse
1096 :
1097 : ! modelled runoff for a given parameter set
1098 : ! dim1=nTimeSteps, dim2=nGauges
1099 0 : real(dp), allocatable, dimension(:, :) :: runoff
1100 :
1101 : ! gauges counter
1102 : integer(i4) :: gg
1103 :
1104 : integer(i4) :: nGaugesTotal
1105 :
1106 : ! aggregated simulated runoff
1107 0 : real(dp), dimension(:), allocatable :: runoff_agg
1108 :
1109 : ! measured runoff
1110 0 : real(dp), dimension(:), allocatable :: runoff_obs
1111 :
1112 : ! mask for measured runoff
1113 0 : logical, dimension(:), allocatable :: runoff_obs_mask
1114 :
1115 :
1116 0 : call eval(parameterset, runoff = runoff)
1117 0 : nGaugesTotal = size(runoff, dim = 2)
1118 :
1119 0 : objective_lnnse = 0.0_dp
1120 0 : do gg = 1, nGaugesTotal
1121 : ! extract runoff
1122 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1123 : ! lnNSE
1124 : objective_lnnse = objective_lnnse + &
1125 0 : lnnse(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1126 : end do
1127 : #ifndef MPI
1128 : ! objective function value which will be minimized
1129 0 : objective_lnnse = 1.0_dp - objective_lnnse / real(nGaugesTotal, dp)
1130 :
1131 0 : call message('objective_lnnse = ', num2str(objective_lnnse))
1132 : ! pause
1133 : #endif
1134 :
1135 0 : deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
1136 :
1137 0 : END FUNCTION objective_lnnse
1138 :
1139 : ! ------------------------------------------------------------------
1140 :
1141 : ! NAME
1142 : ! objective_sse
1143 :
1144 : ! PURPOSE
1145 : !> \brief Objective function of SSE.
1146 :
1147 : !> \details The objective function only depends on a parameter vector.
1148 : !> The model will be called with that parameter vector and
1149 : !> the model output is subsequently compared to observed data.
1150 : !> Therefore, the sum squared errors
1151 : !> \f[ SSE = \sum_{i=1}^N (Q_{obs}(i) - Q_{model}(i))^2 \f]
1152 : !> is calculated and the objective function is
1153 : !> \f[ obj\_value = SSE \f]
1154 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1155 :
1156 : ! INTENT(IN)
1157 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1158 : !> \param[in] "procedure(eval_interface) :: eval"
1159 :
1160 : ! RETURN
1161 : !> \return real(dp) :: objective_sse — objective function value
1162 : !> (which will be e.g. minimized by an optimization routine like DDS)
1163 :
1164 : ! HISTORY
1165 : !> \authors Juliane Mai and Matthias Cuntz
1166 :
1167 : !> \date March 2014
1168 :
1169 : ! Modifications:
1170 : ! Stephan Thober Jan 2015 - introduced extract_runoff
1171 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2) to not interfere with mRM
1172 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1173 :
1174 0 : FUNCTION objective_sse(parameterset, eval)
1175 :
1176 0 : use mo_errormeasures, only : sse
1177 :
1178 : implicit none
1179 :
1180 : real(dp), dimension(:), intent(in) :: parameterset
1181 :
1182 : procedure(eval_interface), INTENT(IN), pointer :: eval
1183 :
1184 : real(dp) :: objective_sse
1185 :
1186 : ! modelled runoff for a given parameter set
1187 : ! dim1=nTimeSteps, dim2=nGauges
1188 0 : real(dp), allocatable, dimension(:, :) :: runoff
1189 :
1190 : ! gauges counter
1191 : integer(i4) :: gg
1192 :
1193 : integer(i4) :: nGaugesTotal
1194 :
1195 : ! aggregated simulated runoff
1196 0 : real(dp), dimension(:), allocatable :: runoff_agg
1197 :
1198 : ! measured runoff
1199 0 : real(dp), dimension(:), allocatable :: runoff_obs
1200 :
1201 : ! mask for measured runoff
1202 0 : logical, dimension(:), allocatable :: runoff_obs_mask
1203 :
1204 :
1205 0 : call eval(parameterset, runoff = runoff)
1206 0 : nGaugesTotal = size(runoff, dim = 2)
1207 :
1208 0 : objective_sse = 0.0_dp
1209 0 : do gg = 1, nGaugesTotal
1210 : !
1211 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1212 : !
1213 : objective_sse = objective_sse + &
1214 0 : sse(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1215 : end do
1216 : #ifndef MPI
1217 : ! objective_sse = objective_sse + sse(gauge%Q, runoff_model_agg) !, runoff_model_agg_mask)
1218 0 : objective_sse = objective_sse / real(nGaugesTotal, dp)
1219 :
1220 0 : call message('objective_sse = ', num2str(objective_sse))
1221 : ! pause
1222 : #endif
1223 :
1224 0 : deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
1225 :
1226 0 : END FUNCTION objective_sse
1227 :
1228 : ! ------------------------------------------------------------------
1229 :
1230 : ! NAME
1231 : ! objective_nse
1232 :
1233 : ! PURPOSE
1234 : !> \brief Objective function of NSE.
1235 :
1236 : !> \details The objective function only depends on a parameter vector.
1237 : !> The model will be called with that parameter vector and
1238 : !> the model output is subsequently compared to observed data.
1239 : !> Therefore, the Nash-Sutcliffe model efficiency coefficient \f$ NSE \f$
1240 : !> \f[ NSE = 1 - \frac{\sum_{i=1}^N (Q_{obs}(i) - Q_{model}(i))^2}
1241 : !> {\sum_{i=1}^N (Q_{obs}(i) - \bar{Q_{obs}})^2} \f]
1242 : !> is calculated and the objective function is
1243 : !> \f[ obj\_value = 1-NSE \f]
1244 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1245 :
1246 : ! INTENT(IN)
1247 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1248 : !> \param[in] "procedure(eval_interface) :: eval"
1249 :
1250 : ! RETURN
1251 : !> \return real(dp) :: objective_nse — objective function value
1252 : !> (which will be e.g. minimized by an optimization routine like DDS)
1253 :
1254 : ! HISTORY
1255 : !> \authors Juliane Mai
1256 :
1257 : !> \date May 2013
1258 :
1259 : ! Modifications:
1260 : ! Stephan Thober Jan 2015 - introduced extract runoff
1261 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2) to not interfere with mRM
1262 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1263 :
1264 0 : FUNCTION objective_nse(parameterset, eval)
1265 :
1266 0 : use mo_errormeasures, only : nse
1267 :
1268 : implicit none
1269 :
1270 : real(dp), dimension(:), intent(in) :: parameterset
1271 :
1272 : procedure(eval_interface), INTENT(IN), pointer :: eval
1273 :
1274 : real(dp) :: objective_nse
1275 :
1276 : ! modelled runoff for a given parameter set
1277 : ! dim1=nTimeSteps, dim2=nGauges
1278 0 : real(dp), allocatable, dimension(:, :) :: runoff
1279 :
1280 : ! gauges counter
1281 : integer(i4) :: gg
1282 :
1283 : integer(i4) :: nGaugesTotal
1284 :
1285 : ! aggregated simulated runoff
1286 0 : real(dp), dimension(:), allocatable :: runoff_agg
1287 :
1288 : ! measured runoff
1289 0 : real(dp), dimension(:), allocatable :: runoff_obs
1290 :
1291 : ! mask for aggregated measured runoff
1292 0 : logical, dimension(:), allocatable :: runoff_obs_mask
1293 :
1294 :
1295 0 : call eval(parameterset, runoff = runoff)
1296 0 : nGaugesTotal = size(runoff, dim = 2)
1297 :
1298 0 : objective_nse = 0.0_dp
1299 0 : do gg = 1, nGaugesTotal
1300 : !
1301 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1302 : !
1303 : objective_nse = objective_nse + &
1304 0 : nse(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1305 : end do
1306 : #ifndef MPI
1307 : ! objective_nse = objective_nse + nse(gauge%Q, runoff_model_agg) !, runoff_model_agg_mask)
1308 0 : objective_nse = 1.0_dp - objective_nse / real(nGaugesTotal, dp)
1309 :
1310 0 : call message('objective_nse (i.e., 1 - NSE) = ', num2str(objective_nse))
1311 : ! pause
1312 : #endif
1313 :
1314 0 : deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
1315 :
1316 0 : END FUNCTION objective_nse
1317 :
1318 : ! ------------------------------------------------------------------
1319 :
1320 : ! NAME
1321 : ! objective_equal_nse_lnnse
1322 :
1323 : ! PURPOSE
1324 : !> \brief Objective function equally weighting NSE and lnNSE.
1325 :
1326 : !> \details The objective function only depends on a parameter vector.
1327 : !> The model will be called with that parameter vector and
1328 : !> the model output is subsequently compared to observed data.
1329 : !> Therefore, the Nash-Sutcliffe model efficiency coefficient \f$ NSE \f$
1330 : !> \f[ NSE = 1 - \frac{\sum_{i=1}^N (Q_{obs}(i) - Q_{model}(i))^2}
1331 : !> {\sum_{i=1}^N (Q_{obs}(i) - \bar{Q_{obs}})^2} \f]
1332 : !> and the logarithmic Nash-Sutcliffe model efficiency coefficient \f$ lnNSE \f$
1333 : !> \f[ lnNSE = 1 - \frac{\sum_{i=1}^N (\ln Q_{obs}(i) - \ln Q_{model}(i))^2}
1334 : !> {\sum_{i=1}^N (\ln Q_{obs}(i) - \bar{\ln Q_{obs}})^2} \f]
1335 : !> are calculated and added up equally weighted:
1336 : !> \f[ obj\_value = \frac{1}{2} (NSE + lnNSE) \f]
1337 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1338 :
1339 : ! INTENT(IN)
1340 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1341 : !> \param[in] "procedure(eval_interface) :: eval"
1342 :
1343 : ! RETURN
1344 : !> \return real(dp) :: objective_equal_nse_lnnse — objective function value
1345 : !> (which will be e.g. minimized by an optimization routine like DDS)
1346 :
1347 : ! HISTORY
1348 : !> \authors Juliane Mai
1349 :
1350 : !> \date May 2013
1351 :
1352 : ! Modifications:
1353 : ! Stephan Thober Jan 2015 - introduced extract_runoff
1354 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2) to not interfere with mRM
1355 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1356 :
1357 6 : FUNCTION objective_equal_nse_lnnse(parameterset, eval)
1358 :
1359 0 : use mo_errormeasures, only : lnnse, nse
1360 :
1361 : implicit none
1362 :
1363 : real(dp), dimension(:), intent(in) :: parameterset
1364 :
1365 : procedure(eval_interface), INTENT(IN), pointer :: eval
1366 :
1367 : real(dp) :: objective_equal_nse_lnnse
1368 :
1369 : ! modelled runoff for a given parameter set
1370 : ! dim2=nGauges
1371 6 : real(dp), allocatable, dimension(:, :) :: runoff
1372 :
1373 : ! gauges counter
1374 : integer(i4) :: gg
1375 :
1376 : integer(i4) :: nGaugesTotal
1377 :
1378 : ! measured runoff
1379 6 : real(dp), dimension(:), allocatable :: runoff_obs
1380 :
1381 : ! aggregated simulated runoff
1382 6 : real(dp), dimension(:), allocatable :: runoff_agg
1383 :
1384 : ! mask for aggregated measured runoff
1385 6 : logical, dimension(:), allocatable :: runoff_obs_mask
1386 :
1387 :
1388 6 : call eval(parameterset, runoff = runoff)
1389 6 : nGaugesTotal = size(runoff, dim = 2)
1390 :
1391 6 : objective_equal_nse_lnnse = 0.0_dp
1392 12 : do gg = 1, nGaugesTotal
1393 : ! extract runoff
1394 6 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1395 : !
1396 : ! NSE
1397 : objective_equal_nse_lnnse = objective_equal_nse_lnnse + &
1398 6 : 0.5_dp * nse(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1399 : ! lnNSE
1400 : objective_equal_nse_lnnse = objective_equal_nse_lnnse + &
1401 18 : 0.5_dp * lnnse(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1402 : end do
1403 : #ifndef MPI
1404 : ! objective function value which will be minimized
1405 6 : objective_equal_nse_lnnse = 1.0_dp - objective_equal_nse_lnnse / real(nGaugesTotal, dp)
1406 :
1407 6 : call message('objective_equal_nse_lnnse = ', num2str(objective_equal_nse_lnnse))
1408 : #endif
1409 :
1410 : ! clean up
1411 6 : deallocate(runoff_agg, runoff_obs)
1412 6 : deallocate(runoff_obs_mask)
1413 :
1414 6 : END FUNCTION objective_equal_nse_lnnse
1415 :
1416 :
1417 : ! ------------------------------------------------------------------
1418 :
1419 : ! NAME
1420 : ! multi_objective_nse_lnnse
1421 :
1422 : ! PURPOSE
1423 : !> \brief Multi-objective function with NSE and lnNSE.
1424 :
1425 : !> \details The objective function only depends on a parameter vector.
1426 : !> The model will be called with that parameter vector and
1427 : !> the model output is subsequently compared to observed data.
1428 : !> Therefore, the Nash-Sutcliffe model efficiency coefficient \f$ NSE \f$
1429 : !> \f[ NSE = 1 - \frac{\sum_{i=1}^N (Q_{obs}(i) - Q_{model}(i))^2}
1430 : !> {\sum_{i=1}^N (Q_{obs}(i) - \bar{Q_{obs}})^2} \f]
1431 : !> and the logarithmic Nash-Sutcliffe model efficiency coefficient \f$ lnNSE \f$
1432 : !> \f[ lnNSE = 1 - \frac{\sum_{i=1}^N (\ln Q_{obs}(i) - \ln Q_{model}(i))^2}
1433 : !> {\sum_{i=1}^N (\ln Q_{obs}(i) - \bar{\ln Q_{obs}})^2} \f]
1434 : !> are calculated and both returned.
1435 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1436 : !> To calibrate this objective you need a multi-objective optimizer like PA-DDS.
1437 :
1438 : ! INTENT(IN)
1439 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1440 : !> \param[in] "procedure(eval_interface) :: eval"
1441 :
1442 : ! RETURN
1443 : !> \return real(dp), dimension(2) :: multi_objective_nse_lnnse — objective function value
1444 : !> (which will be e.g. minimized by an optimization routine like PA-DDS)
1445 :
1446 : ! HISTORY
1447 : !> \authors Juliane Mai
1448 :
1449 : !> \date Oct 2015
1450 :
1451 : ! Modifications:
1452 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1453 :
1454 0 : FUNCTION multi_objective_nse_lnnse(parameterset, eval)
1455 :
1456 6 : use mo_errormeasures, only : lnnse, nse
1457 :
1458 : implicit none
1459 :
1460 : real(dp), dimension(:), intent(in) :: parameterset
1461 :
1462 : procedure(eval_interface), INTENT(IN), pointer :: eval
1463 :
1464 : real(dp), dimension(2) :: multi_objective_nse_lnnse
1465 :
1466 : ! modelled runoff for a given parameter set
1467 : ! dim2=nGauges
1468 0 : real(dp), allocatable, dimension(:, :) :: runoff
1469 :
1470 : ! gauges counter
1471 : integer(i4) :: gg
1472 :
1473 : integer(i4) :: nGaugesTotal
1474 :
1475 : ! measured runoff
1476 0 : real(dp), dimension(:), allocatable :: runoff_obs
1477 :
1478 : ! aggregated simulated runoff
1479 0 : real(dp), dimension(:), allocatable :: runoff_agg
1480 :
1481 : ! mask for aggregated measured runoff
1482 0 : logical, dimension(:), allocatable :: runoff_obs_mask
1483 :
1484 :
1485 : ! call mhm_eval(parameterset, runoff=runoff)
1486 0 : call eval(parameterset, runoff = runoff)
1487 0 : nGaugesTotal = size(runoff, dim = 2)
1488 :
1489 0 : multi_objective_nse_lnnse = 0.0_dp
1490 0 : do gg = 1, nGaugesTotal
1491 : ! extract runoff
1492 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1493 : !
1494 : ! NSE
1495 : multi_objective_nse_lnnse(1) = multi_objective_nse_lnnse(1) + &
1496 0 : nse(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1497 : ! lnNSE
1498 : multi_objective_nse_lnnse(2) = multi_objective_nse_lnnse(2) + &
1499 0 : lnnse(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1500 : end do
1501 : ! objective function value which will be minimized
1502 0 : multi_objective_nse_lnnse(:) = 1.0_dp - multi_objective_nse_lnnse(:) / real(nGaugesTotal, dp)
1503 :
1504 : ! write(*,*) 'multi_objective_nse_lnnse = ',multi_objective_nse_lnnse
1505 :
1506 : ! clean up
1507 0 : deallocate(runoff_agg, runoff_obs)
1508 0 : deallocate(runoff_obs_mask)
1509 :
1510 0 : END FUNCTION multi_objective_nse_lnnse
1511 :
1512 : ! ------------------------------------------------------------------
1513 :
1514 : ! NAME
1515 : ! multi_objective_lnnse_highflow_lnnse_lowflow
1516 :
1517 : ! PURPOSE
1518 : !> \brief Multi-objective function with NSE and lnNSE.
1519 :
1520 : !> \details The objective function only depends on a parameter vector.
1521 : !> The model will be called with that parameter vector and
1522 : !> the model output is subsequently compared to observed data.
1523 :
1524 : !> A timepoint \f$t\f$ of the observed data is marked as a lowflow timepoint \f$t_{low}\f$ if
1525 : !> \f[ Q_{obs}(t) < min(Q_{obs}) + 0.05 * ( max(Q_{obs}) - min(Q_{obs}) )\f]
1526 : !> and a timepoint \f$t\f$ of the observed data is marked as a highflow timepoint \f$t_{high}\f$ if
1527 : !> \f[ t_{high} if Q_{obs}(i) > percentile(Q_{obs},95.)\f]
1528 : !> This timepoint identification is only performed for the observed data.
1529 :
1530 : !> The first objective is the logarithmic Nash-Sutcliffe model efficiency coefficient \f$ lnNSE_{high} \f$
1531 : !> of discharge values at high-flow timepoints
1532 : !> \f[ lnNSE_{high} = 1 - \frac{\sum_{i=1}^{N_{high}} (\ln Q_{obs}(i) - \ln Q_{model}(i))^2}
1533 : !> {\sum_{i=1}^N (\ln Q_{obs}(i) - \bar{\ln Q_{obs}})^2} \f] .
1534 : !> The second objective is the logarithmic Nash-Sutcliffe model efficiency coefficient \f$ lnNSE_{low} \f$
1535 : !> of discharge values at low-flow timepoints
1536 : !> \f[ lnNSE_{low} = 1 - \frac{\sum_{i=1}^{N_{low}} (\ln Q_{obs}(i) - \ln Q_{model}(i))^2}
1537 : !> {\sum_{i=1}^N (\ln Q_{obs}(i) - \bar{\ln Q_{obs}})^2} \f] .
1538 : !> Both objectives are returned.
1539 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1540 : !> To calibrate this objective you need a multi-objective optimizer like PA-DDS.
1541 :
1542 : ! INTENT(IN)
1543 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1544 : !> \param[in] "procedure(eval_interface) :: eval"
1545 :
1546 : ! RETURN
1547 : !> \return real(dp), dimension(2) :: multi_objective_lnnse_highflow_lnnse_lowflow — objective function
1548 : !> value
1549 : !> (which will be e.g. minimized by an optimization routine like PA-DDS)
1550 :
1551 : ! HISTORY
1552 : !> \authors Juliane Mai
1553 :
1554 : !> \date Oct 2015
1555 :
1556 : ! Modifications:
1557 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1558 :
1559 0 : FUNCTION multi_objective_lnnse_highflow_lnnse_lowflow(parameterset, eval)
1560 :
1561 0 : use mo_errormeasures, only : lnnse
1562 : use mo_percentile, only : percentile
1563 :
1564 : implicit none
1565 :
1566 : real(dp), dimension(:), intent(in) :: parameterset
1567 :
1568 : procedure(eval_interface), INTENT(IN), pointer :: eval
1569 :
1570 : real(dp), dimension(2) :: multi_objective_lnnse_highflow_lnnse_lowflow
1571 :
1572 : ! modelled runoff for a given parameter set
1573 0 : real(dp), dimension(:, :), allocatable :: runoff
1574 :
1575 : ! gauges counter
1576 : integer(i4) :: gg
1577 :
1578 : integer(i4) :: nGaugesTotal
1579 :
1580 : ! measured runoff
1581 0 : real(dp), dimension(:), allocatable :: runoff_obs
1582 :
1583 : ! aggregated simulated runoff
1584 0 : real(dp), dimension(:), allocatable :: runoff_agg
1585 :
1586 : ! mask for aggregated measured runoff
1587 0 : logical, dimension(:), allocatable :: runoff_obs_mask
1588 :
1589 : ! upper discharge value to determine lowflow timepoints
1590 0 : real(dp) :: q_low
1591 :
1592 : ! lower discharge value to determine highflow timepoints
1593 0 : real(dp) :: q_high
1594 :
1595 : ! total number of discharge values
1596 : integer(i4) :: nrunoff
1597 :
1598 : ! timepoint counter
1599 : integer(i4) :: tt
1600 :
1601 : ! mask to get lowflow values
1602 0 : logical, dimension(:), allocatable :: lowflow_mask
1603 :
1604 : ! mask to get highflow values
1605 0 : logical, dimension(:), allocatable :: highflow_mask
1606 :
1607 :
1608 : ! call mhm_eval(parameterset, runoff=runoff)
1609 0 : call eval(parameterset, runoff = runoff)
1610 0 : nGaugesTotal = size(runoff, dim = 2)
1611 :
1612 0 : multi_objective_lnnse_highflow_lnnse_lowflow = 0.0_dp
1613 0 : do gg = 1, nGaugesTotal
1614 : !
1615 : ! extract runoff
1616 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1617 0 : nrunoff = size(runoff_obs, dim = 1)
1618 : !print*, 'nrunoff = ',nrunoff
1619 : !
1620 : ! mask for highflow timepoints
1621 0 : if (allocated(highflow_mask)) deallocate(highflow_mask)
1622 0 : allocate(highflow_mask(nrunoff))
1623 0 : highflow_mask = .false.
1624 0 : q_high = percentile(runoff_obs, 95._dp, mask = runoff_obs_mask)
1625 0 : forall(tt = 1 : nrunoff) highflow_mask(tt) = (runoff_obs(tt) > q_high .and. runoff_obs_mask(tt))
1626 : !print*, 'nhigh = ',count(highflow_mask)
1627 : !
1628 : ! mask for lowflow timepoints
1629 0 : if (allocated(lowflow_mask)) deallocate(lowflow_mask)
1630 0 : allocate(lowflow_mask(nrunoff))
1631 0 : lowflow_mask = .false.
1632 0 : q_low = minval(runoff_obs, mask = runoff_obs_mask) &
1633 0 : + 0.05_dp * (maxval(runoff_obs, mask = runoff_obs_mask) - minval(runoff_obs, mask = runoff_obs_mask))
1634 0 : forall(tt = 1 : nrunoff) lowflow_mask(tt) = (runoff_obs(tt) < q_low .and. runoff_obs_mask(tt))
1635 : !print*, 'nlow = ',count(lowflow_mask)
1636 : !
1637 : ! lnNSE highflows
1638 : multi_objective_lnnse_highflow_lnnse_lowflow(1) = multi_objective_lnnse_highflow_lnnse_lowflow(1) + &
1639 0 : lnnse(runoff_obs, runoff_agg, mask = highflow_mask)
1640 : !
1641 : ! lnNSE lowflows
1642 : multi_objective_lnnse_highflow_lnnse_lowflow(2) = multi_objective_lnnse_highflow_lnnse_lowflow(2) + &
1643 0 : lnnse(runoff_obs, runoff_agg, mask = lowflow_mask)
1644 : end do
1645 : ! objective function value which will be minimized
1646 : multi_objective_lnnse_highflow_lnnse_lowflow(:) = 1.0_dp &
1647 0 : - multi_objective_lnnse_highflow_lnnse_lowflow(:) / real(nGaugesTotal, dp)
1648 :
1649 : ! write(*,*) 'multi_objective_lnnse_highflow_lnnse_lowflow = ',multi_objective_lnnse_highflow_lnnse_lowflow
1650 :
1651 : ! clean up
1652 0 : deallocate(runoff_agg, runoff_obs)
1653 0 : deallocate(runoff_obs_mask)
1654 :
1655 0 : END FUNCTION multi_objective_lnnse_highflow_lnnse_lowflow
1656 :
1657 : ! ------------------------------------------------------------------
1658 :
1659 : ! NAME
1660 : ! multi_objective_lnnse_highflow_lnnse_lowflow_2
1661 :
1662 : ! PURPOSE
1663 : !> \brief Multi-objective function with NSE and lnNSE.
1664 :
1665 : !> \details The objective function only depends on a parameter vector.
1666 : !> The model will be called with that parameter vector and
1667 : !> the model output is subsequently compared to observed data.
1668 :
1669 : !> A timepoint \f$t\f$ of the observed data is marked as a lowflow timepoint \f$t_{low}\f$ if
1670 : !> \f[ Q_{obs}(t) < min(Q_{obs}) + 0.05 * ( max(Q_{obs}) - min(Q_{obs}) )\f]
1671 : !> and all other timepoints are marked as a highflow timepoints \f$t_{high}\f$.
1672 : !> This timepoint identification is only performed for the observed data.
1673 :
1674 : !> The first objective is the logarithmic Nash-Sutcliffe model efficiency coefficient \f$ lnNSE_{high} \f$
1675 : !> of discharge values at high-flow timepoints
1676 : !> \f[ lnNSE_{high} = 1 - \frac{\sum_{i=1}^{N_{high}} (\ln Q_{obs}(i) - \ln Q_{model}(i))^2}
1677 : !> {\sum_{i=1}^N (\ln Q_{obs}(i) - \bar{\ln Q_{obs}})^2} \f] .
1678 : !> The second objective is the logarithmic Nash-Sutcliffe model efficiency coefficient \f$ lnNSE_{low} \f$
1679 : !> of discharge values at low-flow timepoints
1680 : !> \f[ lnNSE_{low} = 1 - \frac{\sum_{i=1}^{N_{low}} (\ln Q_{obs}(i) - \ln Q_{model}(i))^2}
1681 : !> {\sum_{i=1}^N (\ln Q_{obs}(i) - \bar{\ln Q_{obs}})^2} \f] .
1682 : !> Both objectives are returned.
1683 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1684 : !> To calibrate this objective you need a multi-objective optimizer like PA-DDS.
1685 :
1686 : ! INTENT(IN)
1687 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1688 : !> \param[in] "procedure(eval_interface) :: eval"
1689 :
1690 : ! RETURN
1691 : !> \return real(dp), dimension(2) :: multi_objective_lnnse_highflow_lnnse_lowflow_2 — objective function
1692 : !> value
1693 : !> (which will be e.g. minimized by an optimization routine like PA-DDS)
1694 :
1695 : ! HISTORY
1696 : !> \authors Juliane Mai
1697 :
1698 : !> \date Oct 2015
1699 :
1700 : ! Modifications:
1701 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1702 :
1703 0 : FUNCTION multi_objective_lnnse_highflow_lnnse_lowflow_2(parameterset, eval)
1704 :
1705 0 : use mo_errormeasures, only : lnnse
1706 :
1707 : implicit none
1708 :
1709 : real(dp), dimension(:), intent(in) :: parameterset
1710 :
1711 : procedure(eval_interface), INTENT(IN), pointer :: eval
1712 :
1713 : real(dp), dimension(2) :: multi_objective_lnnse_highflow_lnnse_lowflow_2
1714 :
1715 : ! modelled runoff for a given parameter set
1716 0 : real(dp), dimension(:, :), allocatable :: runoff
1717 :
1718 : ! gauges counter
1719 : integer(i4) :: gg
1720 :
1721 : integer(i4) :: nGaugesTotal
1722 :
1723 : ! measured runoff
1724 0 : real(dp), dimension(:), allocatable :: runoff_obs
1725 :
1726 : ! aggregated simulated runoff
1727 0 : real(dp), dimension(:), allocatable :: runoff_agg
1728 :
1729 : ! mask for aggregated measured runoff
1730 0 : logical, dimension(:), allocatable :: runoff_obs_mask
1731 :
1732 : ! upper discharge value to determine lowflow timepoints
1733 0 : real(dp) :: q_low
1734 :
1735 : ! total number of discharge values
1736 : integer(i4) :: nrunoff
1737 :
1738 : ! timepoint counter
1739 : integer(i4) :: tt
1740 :
1741 : ! mask to get lowflow values
1742 0 : logical, dimension(:), allocatable :: lowflow_mask
1743 :
1744 : ! mask to get highflow values
1745 0 : logical, dimension(:), allocatable :: highflow_mask
1746 :
1747 :
1748 : ! call mhm_eval(parameterset, runoff=runoff)
1749 0 : call eval(parameterset, runoff = runoff)
1750 0 : nGaugesTotal = size(runoff, dim = 2)
1751 :
1752 0 : multi_objective_lnnse_highflow_lnnse_lowflow_2 = 0.0_dp
1753 0 : do gg = 1, nGaugesTotal
1754 : !
1755 : ! extract runoff
1756 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1757 0 : nrunoff = size(runoff_obs, dim = 1)
1758 : !print*, 'nrunoff = ',nrunoff
1759 : !
1760 : ! mask for lowflow timepoints
1761 0 : if (allocated(lowflow_mask)) deallocate(lowflow_mask)
1762 0 : allocate(lowflow_mask(nrunoff))
1763 0 : lowflow_mask = .false.
1764 0 : q_low = minval(runoff_obs, mask = runoff_obs_mask) &
1765 0 : + 0.05_dp * (maxval(runoff_obs, mask = runoff_obs_mask) - minval(runoff_obs, mask = runoff_obs_mask))
1766 0 : forall(tt = 1 : nrunoff) lowflow_mask(tt) = (runoff_obs(tt) < q_low .and. runoff_obs_mask(tt))
1767 : !print*, 'nlow = ',count(lowflow_mask)
1768 : !
1769 : ! mask for highflow timepoints
1770 0 : if (allocated(highflow_mask)) deallocate(highflow_mask)
1771 0 : allocate(highflow_mask(nrunoff))
1772 0 : highflow_mask = (.not. lowflow_mask) .and. runoff_obs_mask
1773 : !print*, 'nhigh = ',count(highflow_mask)
1774 : !
1775 : ! lnNSE highflows
1776 : multi_objective_lnnse_highflow_lnnse_lowflow_2(1) = multi_objective_lnnse_highflow_lnnse_lowflow_2(1) + &
1777 0 : lnnse(runoff_obs, runoff_agg, mask = highflow_mask)
1778 : !
1779 : ! lnNSE lowflows
1780 : multi_objective_lnnse_highflow_lnnse_lowflow_2(2) = multi_objective_lnnse_highflow_lnnse_lowflow_2(2) + &
1781 0 : lnnse(runoff_obs, runoff_agg, mask = lowflow_mask)
1782 : end do
1783 : ! objective function value which will be minimized
1784 : multi_objective_lnnse_highflow_lnnse_lowflow_2(:) = 1.0_dp &
1785 0 : - multi_objective_lnnse_highflow_lnnse_lowflow_2(:) / real(nGaugesTotal, dp)
1786 :
1787 : ! write(*,*) 'multi_objective_lnnse_highflow_lnnse_lowflow_2 = ',multi_objective_lnnse_highflow_lnnse_lowflow_2
1788 :
1789 : ! clean up
1790 0 : deallocate(runoff_agg, runoff_obs)
1791 0 : deallocate(runoff_obs_mask)
1792 :
1793 0 : END FUNCTION multi_objective_lnnse_highflow_lnnse_lowflow_2
1794 :
1795 : ! ------------------------------------------------------------------
1796 :
1797 : ! NAME
1798 : ! multi_objective_ae_fdc_lsv_nse_djf
1799 :
1800 : ! PURPOSE
1801 : !> \brief Multi-objective function with absolute error of Flow Duration Curves
1802 : !> low-segment volume and nse of DJF's discharge.
1803 :
1804 : !> \details The objective function only depends on a parameter vector.
1805 : !> The model will be called with that parameter vector and
1806 : !> the model output is subsequently compared to observed data.
1807 :
1808 : !> The first objective is using the routine "FlowDurationCurves" from "mo_signatures" to determine the
1809 : !> low-segment volume of the FDC. The objective is the absolute difference between the observed volume
1810 : !> and the simulated volume.
1811 :
1812 : !> For the second objective the discharge of the winter months December, January and February are extracted
1813 : !> from the time series. The objective is then the Nash-Sutcliffe efficiency NSE of the observed winter
1814 : !> discharge against the simulated winter discharge.
1815 :
1816 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1817 : !> To calibrate this objective you need a multi-objective optimizer like PA-DDS.
1818 :
1819 : ! INTENT(IN)
1820 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1821 : !> \param[in] "procedure(eval_interface) :: eval"
1822 :
1823 : ! RETURN
1824 : !> \return real(dp), dimension(2) :: multi_objective_ae_fdc_lsv_nse_djf — objective function value
1825 : !> (which will be e.g. minimized by an optimization routine like PA-DDS)
1826 :
1827 : ! HISTORY
1828 : !> \authors Juliane Mai
1829 :
1830 : !> \date Feb 2016
1831 :
1832 : ! Modifications:
1833 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1834 :
1835 0 : FUNCTION multi_objective_ae_fdc_lsv_nse_djf(parameterset, eval)
1836 :
1837 0 : use mo_common_mhm_mrm_variables, only : evalPer
1838 : use mo_errormeasures, only : nse
1839 : use mo_julian, only : dec2date
1840 : use mo_mrm_global_variables, only : gauge, nMeasPerDay
1841 : use mo_mrm_signatures, only : FlowDurationCurve
1842 :
1843 : implicit none
1844 :
1845 : real(dp), dimension(:), intent(in) :: parameterset
1846 :
1847 : procedure(eval_interface), INTENT(IN), pointer :: eval
1848 :
1849 : real(dp), dimension(2) :: multi_objective_ae_fdc_lsv_nse_djf
1850 :
1851 : ! modelled runoff for a given parameter set
1852 0 : real(dp), dimension(:, :), allocatable :: runoff
1853 :
1854 : ! gauges counter
1855 : integer(i4) :: gg
1856 :
1857 : ! total number of gauges
1858 : integer(i4) :: nGaugesTotal
1859 :
1860 : ! domain ID of gauge
1861 : integer(i4) :: iDomain
1862 :
1863 : ! measured runoff
1864 0 : real(dp), dimension(:), allocatable :: runoff_obs
1865 :
1866 : ! aggregated simulated runoff
1867 0 : real(dp), dimension(:), allocatable :: runoff_agg
1868 :
1869 : ! mask for aggregated measured runoff
1870 0 : logical, dimension(:), allocatable :: runoff_obs_mask
1871 :
1872 : ! total number of discharge values
1873 : integer(i4) :: nrunoff
1874 :
1875 : ! timepoint counter
1876 : integer(i4) :: tt
1877 :
1878 : ! month of current time step
1879 : integer(i4) :: month
1880 :
1881 : ! Fractional Julian day of current time step
1882 0 : real(dp) :: current_time
1883 :
1884 : ! mask to get lowflow values
1885 0 : logical, dimension(:), allocatable :: djf_mask
1886 :
1887 : ! quantiles for FDC
1888 0 : real(dp), dimension(10) :: quantiles
1889 :
1890 : ! number of quantiles
1891 : integer(i4) :: nquantiles
1892 :
1893 : ! FDC of simulated or observed discharge
1894 0 : real(dp), dimension(size(quantiles)) :: fdc
1895 :
1896 : ! low-segment volume of FDC of simulated discharge
1897 0 : real(dp) :: lsv_mod
1898 :
1899 : ! low-segment volume of FDC of observed discharge
1900 0 : real(dp) :: lsv_obs
1901 :
1902 :
1903 : ! call mhm_eval(parameterset, runoff=runoff)
1904 0 : call eval(parameterset, runoff = runoff)
1905 0 : nGaugesTotal = size(runoff, dim = 2)
1906 0 : nquantiles = size(quantiles)
1907 0 : forall(tt = 1 : nquantiles) quantiles(tt) = real(tt, dp) / real(nquantiles, dp)
1908 :
1909 0 : multi_objective_ae_fdc_lsv_nse_djf = 0.0_dp
1910 0 : do gg = 1, nGaugesTotal
1911 : !
1912 : ! extract runoff
1913 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1914 0 : nrunoff = size(runoff_obs, dim = 1)
1915 : !
1916 : ! mask DJF timepoints
1917 0 : if (allocated(djf_mask)) deallocate(djf_mask)
1918 0 : allocate(djf_mask(nrunoff))
1919 0 : djf_mask = .false.
1920 :
1921 0 : iDomain = gauge%domainId(gg)
1922 0 : do tt = 1, nrunoff
1923 0 : current_time = evalPer(iDomain)%julStart + (tt - 1) * 1.0_dp / real(nMeasPerDay, dp)
1924 0 : call dec2date(current_time, mm = month)
1925 0 : if ((month == 1 .or. month == 2 .or. month == 12) .and. runoff_obs_mask(tt)) djf_mask(tt) = .True.
1926 : end do
1927 : !
1928 : ! Absolute error of low-segment volume of FDC
1929 0 : fdc = FlowDurationCurve(runoff_obs, quantiles, mask = runoff_obs_mask, low_segment_volume = lsv_obs)
1930 0 : fdc = FlowDurationCurve(runoff_agg, quantiles, mask = runoff_obs_mask, low_segment_volume = lsv_mod)
1931 0 : fdc = fdc * 1.0_dp ! only to avoid warning of unused variable
1932 : !
1933 : ! Absolute distance between low-segment volumes
1934 : multi_objective_ae_fdc_lsv_nse_djf(1) = multi_objective_ae_fdc_lsv_nse_djf(1) + &
1935 0 : abs(lsv_obs - lsv_mod)
1936 : !
1937 : ! NSE of DJF discharge
1938 : multi_objective_ae_fdc_lsv_nse_djf(2) = multi_objective_ae_fdc_lsv_nse_djf(2) + &
1939 0 : nse(runoff_obs, runoff_agg, mask = djf_mask)
1940 : end do
1941 : ! objective function value which will be minimized
1942 : multi_objective_ae_fdc_lsv_nse_djf(1) = &
1943 0 : multi_objective_ae_fdc_lsv_nse_djf(1) / real(nGaugesTotal, dp)
1944 : multi_objective_ae_fdc_lsv_nse_djf(2) = 1.0_dp &
1945 0 : - multi_objective_ae_fdc_lsv_nse_djf(2) / real(nGaugesTotal, dp)
1946 :
1947 : call message('multi_objective_ae_fdc_lsv_nse_djf = [', &
1948 : num2str(multi_objective_ae_fdc_lsv_nse_djf(1)), ', ', &
1949 0 : num2str(multi_objective_ae_fdc_lsv_nse_djf(2)), ']')
1950 :
1951 : ! clean up
1952 0 : deallocate(runoff_agg, runoff_obs)
1953 0 : deallocate(runoff_obs_mask)
1954 :
1955 0 : END FUNCTION multi_objective_ae_fdc_lsv_nse_djf
1956 :
1957 : ! ------------------------------------------------------------------
1958 :
1959 : ! NAME
1960 : ! objective_power6_nse_lnnse
1961 :
1962 : ! PURPOSE
1963 : !> \brief Objective function of combined NSE and lnNSE with power of 5
1964 : !> i.e. the p-norm with p=5.
1965 :
1966 : !> \details The objective function only depends on a parameter vector.
1967 : !> The model will be called with that parameter vector and
1968 : !> the model output is subsequently compared to observed data.
1969 : !> Therefore, the Nash-Sutcliffe model efficiency coefficient \f$ NSE \f$
1970 : !> \f[ NSE = 1 - \frac{\sum_{i=1}^N (Q_{obs}(i) - Q_{model}(i))^2}
1971 : !> {\sum_{i=1}^N (Q_{obs}(i) - \bar{Q_{obs}})^2} \f]
1972 : !> and the logarithmic Nash-Sutcliffe model efficiency coefficient \f$ lnNSE \f$
1973 : !> \f[ lnNSE = 1 - \frac{\sum_{i=1}^N (\ln Q_{obs}(i) - \ln Q_{model}(i))^2}
1974 : !> {\sum_{i=1}^N (\ln Q_{obs}(i) - \bar{\ln Q_{obs}})^2} \f]
1975 : !> are calculated and added up equally weighted:
1976 : !> \f[ obj\_value = \sqrt[6]{(1-NSE)^6 + (1-lnNSE)^6} \f]
1977 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
1978 :
1979 : ! INTENT(IN)
1980 : !> \param[in] "real(dp), dimension(:) :: parameterset"
1981 : !> \param[in] "procedure(eval_interface) :: eval"
1982 :
1983 : ! RETURN
1984 : !> \return real(dp) :: objective_power6_nse_lnnse — objective function value
1985 : !> (which will be e.g. minimized by an optimization routine like DDS)
1986 :
1987 : ! HISTORY
1988 : !> \authors Juliane Mai and Matthias Cuntz
1989 :
1990 : !> \date March 2014
1991 :
1992 : ! Modifications:
1993 : ! Stephan Thober Jan 2015 - introduced extract_runoff
1994 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2) to not interfere with mRM
1995 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1996 :
1997 0 : FUNCTION objective_power6_nse_lnnse(parameterset, eval)
1998 :
1999 0 : use mo_errormeasures, only : lnnse, nse
2000 :
2001 : implicit none
2002 :
2003 : real(dp), dimension(:), intent(in) :: parameterset
2004 :
2005 : procedure(eval_interface), INTENT(IN), pointer :: eval
2006 :
2007 : real(dp) :: objective_power6_nse_lnnse
2008 :
2009 : ! modelled runoff for a given parameter set
2010 : ! dim1=nTimeSteps, dim2=nGauges
2011 0 : real(dp), allocatable, dimension(:, :) :: runoff
2012 :
2013 : ! gauges counter
2014 : integer(i4) :: gg
2015 :
2016 : integer(i4) :: nGaugesTotal
2017 :
2018 : ! aggregated simulated runoff
2019 0 : real(dp), dimension(:), allocatable :: runoff_agg
2020 :
2021 : ! measured runoff
2022 0 : real(dp), dimension(:), allocatable :: runoff_obs
2023 :
2024 : ! mask for measured runoff
2025 0 : logical, dimension(:), allocatable :: runoff_obs_mask
2026 :
2027 : real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2028 :
2029 :
2030 0 : call eval(parameterset, runoff = runoff)
2031 0 : nGaugesTotal = size(runoff, dim = 2)
2032 :
2033 0 : objective_power6_nse_lnnse = 0.0_dp
2034 0 : do gg = 1, nGaugesTotal
2035 : ! extract runoff
2036 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2037 : ! NSE + lnNSE
2038 : objective_power6_nse_lnnse = objective_power6_nse_lnnse + &
2039 : ((1.0_dp - nse(runoff_obs, runoff_agg, mask = runoff_obs_mask))**6 + &
2040 0 : (1.0_dp - lnnse(runoff_obs, runoff_agg, mask = runoff_obs_mask))**6)**onesixth
2041 : end do
2042 : #ifndef MPI
2043 : ! objective function value which will be minimized
2044 0 : objective_power6_nse_lnnse = objective_power6_nse_lnnse / real(nGaugesTotal, dp)
2045 :
2046 0 : call message('objective_power6_nse_lnnse = ', num2str(objective_power6_nse_lnnse))
2047 : ! pause
2048 : #endif
2049 :
2050 0 : deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
2051 :
2052 0 : END FUNCTION objective_power6_nse_lnnse
2053 :
2054 : ! ------------------------------------------------------------------
2055 :
2056 : ! NAME
2057 : ! objective_kge
2058 :
2059 : ! PURPOSE
2060 : !> \brief Objective function of KGE.
2061 :
2062 : !> \details The objective function only depends on a parameter vector.
2063 : !> The model will be called with that parameter vector and
2064 : !> the model output is subsequently compared to observed data.
2065 :
2066 : !> Therefore, the Kling-Gupta model efficiency coefficient \f$ KGE \f$
2067 : !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f]
2068 : !> where
2069 : !> \f$ r \f$ = Pearson product-moment correlation coefficient
2070 : !> \f$ \alpha \f$ = ratio of similated mean to observed mean
2071 : !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation
2072 : !> is calculated and the objective function is
2073 : !> \f[ obj\_value = 1.0 - KGE \f]
2074 : !> \f$(1-KGE)\f$ is the objective since we always apply minimization methods.
2075 : !> The minimal value of \f$(1-KGE)\f$ is 0 for the optimal KGE of 1.0.
2076 :
2077 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
2078 :
2079 : ! INTENT(IN)
2080 : !> \param[in] "real(dp), dimension(:) :: parameterset"
2081 : !> \param[in] "procedure(eval_interface) :: eval"
2082 :
2083 : ! RETURN
2084 : !> \return real(dp) :: objective_kge — objective function value
2085 : !> (which will be e.g. minimized by an optimization routine like DDS)
2086 :
2087 : ! HISTORY
2088 : !> \authors Rohini Kumar
2089 :
2090 : !> \date August 2014
2091 :
2092 : ! Modifications:
2093 : ! Stephan Thober Jan 2015 - introduced extract_runoff
2094 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2) to not interfere with mRM
2095 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
2096 :
2097 0 : FUNCTION objective_kge(parameterset, eval)
2098 :
2099 0 : use mo_errormeasures, only : kge
2100 :
2101 : implicit none
2102 :
2103 : real(dp), dimension(:), intent(in) :: parameterset
2104 :
2105 : procedure(eval_interface), INTENT(IN), pointer :: eval
2106 :
2107 : real(dp) :: objective_kge
2108 :
2109 : ! modelled runoff for a given parameter set
2110 : ! dim1=nTimeSteps, dim2=nGauges
2111 0 : real(dp), allocatable, dimension(:, :) :: runoff
2112 :
2113 : ! gauges counter
2114 : integer(i4) :: gg
2115 :
2116 : integer(i4) :: nGaugesTotal
2117 :
2118 : ! aggregated simulated runoff
2119 0 : real(dp), dimension(:), allocatable :: runoff_agg
2120 :
2121 : ! measured runoff
2122 0 : real(dp), dimension(:), allocatable :: runoff_obs
2123 :
2124 : ! mask for measured runoff
2125 0 : logical, dimension(:), allocatable :: runoff_obs_mask
2126 :
2127 :
2128 : !
2129 0 : call eval(parameterset, runoff = runoff)
2130 0 : nGaugesTotal = size(runoff, dim = 2)
2131 :
2132 0 : objective_kge = 0.0_dp
2133 0 : do gg = 1, nGaugesTotal
2134 : ! extract runoff
2135 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2136 : ! KGE
2137 : objective_kge = objective_kge + &
2138 0 : kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
2139 : end do
2140 : #ifndef MPI
2141 : ! objective_kge = objective_kge + kge(gauge%Q, runoff_model_agg, runoff_model_agg_mask)
2142 0 : objective_kge = 1.0_dp - objective_kge / real(nGaugesTotal, dp)
2143 :
2144 0 : call message('objective_kge (i.e., 1 - KGE) = ', num2str(objective_kge))
2145 : #endif
2146 : ! pause
2147 :
2148 0 : deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
2149 :
2150 0 : END FUNCTION objective_kge
2151 :
2152 : ! ------------------------------------------------------------------
2153 :
2154 : ! NAME
2155 : ! objective_multiple_gauges_kge_power6
2156 :
2157 : ! PURPOSE
2158 : !> \brief combined objective function based on KGE raised to the power 6
2159 :
2160 : !> \details The objective function only depends on a parameter vector.
2161 : !> The model will be called with that parameter vector and
2162 : !> the model output is subsequently compared to observed data.
2163 :
2164 : !> Therefore, the Kling-Gupta model efficiency coefficient \f$ KGE \f$ for a given gauging station
2165 : !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f]
2166 : !> where
2167 : !> \f$ r \f$ = Pearson product-moment correlation coefficient
2168 : !> \f$ \alpha \f$ = ratio of similated mean to observed mean
2169 : !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation
2170 : !> is calculated and the objective function for a given gauging station (\f$ i \f$) is
2171 : !> \f[ \phi_{i} = 1.0 - KGE_{i} \f]
2172 : !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods.
2173 : !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.
2174 :
2175 : !> Finally, the overall \f$ OF \f$ is estimated based on the power-6 norm to
2176 : !> combine the \f$ \phi_{i} \f$ from all gauging stations (\f$ N \f$).
2177 : !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 } \f].
2178 :
2179 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
2180 :
2181 : ! INTENT(IN)
2182 : !> \param[in] "real(dp), dimension(:) :: parameterset"
2183 : !> \param[in] "procedure(eval_interface) :: eval"
2184 :
2185 : ! RETURN
2186 : !> \return real(dp) :: objective_multiple_gauges_kge_power6 — objective function value
2187 : !> (which will be e.g. minimized by an optimization routine like DDS)
2188 :
2189 : ! HISTORY
2190 : !> \authors Rohini Kumar
2191 :
2192 : !> \date March 2015
2193 :
2194 : ! Modifications:
2195 : ! Stephan Thober Aug 2015 - substituted nGaugesTotal variable with size(runoff, dim=2)
2196 : ! to not interfere with mRM
2197 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
2198 :
2199 0 : FUNCTION objective_multiple_gauges_kge_power6(parameterset, eval)
2200 :
2201 0 : use mo_errormeasures, only : kge
2202 :
2203 : implicit none
2204 :
2205 : real(dp), dimension(:), intent(in) :: parameterset
2206 :
2207 : procedure(eval_interface), INTENT(IN), pointer :: eval
2208 :
2209 : real(dp) :: objective_multiple_gauges_kge_power6
2210 :
2211 : real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2212 :
2213 : ! modelled runoff for a given parameter set
2214 : ! dim1=nTimeSteps, dim2=nGauges
2215 0 : real(dp), allocatable, dimension(:, :) :: runoff
2216 :
2217 : ! gauges counter
2218 : integer(i4) :: gg
2219 :
2220 : integer(i4) :: nGaugesTotal
2221 :
2222 : ! aggregated simulated runoff
2223 0 : real(dp), dimension(:), allocatable :: runoff_agg
2224 :
2225 : ! measured runoff
2226 0 : real(dp), dimension(:), allocatable :: runoff_obs
2227 :
2228 : ! mask for measured runoff
2229 0 : logical, dimension(:), allocatable :: runoff_obs_mask
2230 :
2231 :
2232 : !
2233 0 : call eval(parameterset, runoff = runoff)
2234 0 : nGaugesTotal = size(runoff, dim = 2)
2235 :
2236 0 : objective_multiple_gauges_kge_power6 = 0.0_dp
2237 0 : do gg = 1, nGaugesTotal
2238 : ! extract runoff
2239 0 : call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2240 : ! KGE
2241 : objective_multiple_gauges_kge_power6 = objective_multiple_gauges_kge_power6 + &
2242 0 : ((1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)) / real(nGaugesTotal, dp))**6
2243 : end do
2244 : #ifndef MPI
2245 0 : objective_multiple_gauges_kge_power6 = objective_multiple_gauges_kge_power6**onesixth
2246 0 : call message('objective_multiple_gauges_kge_power6 = ', num2str(objective_multiple_gauges_kge_power6))
2247 : #endif
2248 :
2249 0 : deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
2250 :
2251 0 : END FUNCTION objective_multiple_gauges_kge_power6
2252 :
2253 : ! ------------------------------------------------------------------
2254 :
2255 : ! NAME
2256 : ! objective_weighted_nse
2257 :
2258 : ! PURPOSE
2259 : !> \brief Objective function of weighted NSE.
2260 :
2261 : !> \details The objective function only depends on a parameter vector.
2262 : !> The model will be called with that parameter vector and
2263 : !> the model output is subsequently compared to observed data.
2264 : !> Therefore, the weighted Nash-Sutcliffe model efficiency coefficient \f$ NSE \f$
2265 : !> \f[ wNSE = 1 - \frac{\sum_{i=1}^N Q_{obs}(i) * (Q_{obs}(i) - Q_{model}(i))^2}
2266 : !> {\sum_{i=1}^N Q_{obs}(i) * (Q_{obs}(i) - \bar{Q_{obs}})^2} \f]
2267 : !> is calculated and the objective function is
2268 : !> \f[ obj\_value = 1- wNSE \f]
2269 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
2270 :
2271 : ! INTENT(IN)
2272 : !> \param[in] "real(dp), dimension(:) :: parameterset"
2273 : !> \param[in] "procedure(eval_interface) :: eval"
2274 :
2275 : ! RETURN
2276 : !> \return real(dp) :: objective_weighted_nse — objective function value
2277 : !> (which will be e.g. minimized by an optimization routine like DDS)
2278 :
2279 : ! HISTORY
2280 : !> \authors Stephan Thober, Bjoern Guse
2281 :
2282 : !> \date May 2018
2283 :
2284 : ! Modifications:
2285 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
2286 :
2287 0 : FUNCTION objective_weighted_nse(parameterset, eval)
2288 :
2289 0 : use mo_errormeasures, only : wnse
2290 :
2291 : implicit none
2292 :
2293 : real(dp), dimension(:), intent(in) :: parameterset
2294 :
2295 : procedure(eval_interface), INTENT(IN), pointer :: eval
2296 :
2297 : real(dp) :: objective_weighted_nse
2298 :
2299 : ! modelled runoff for a given parameter set
2300 : ! dim1=nTimeSteps, dim2=nGauges
2301 0 : real(dp), allocatable, dimension(:,:) :: runoff
2302 :
2303 : ! gauges counter
2304 : integer(i4) :: gg
2305 :
2306 : integer(i4) :: nGaugesTotal
2307 :
2308 : ! aggregated simulated runoff
2309 0 : real(dp), dimension(:), allocatable :: runoff_agg
2310 :
2311 : ! measured runoff
2312 0 : real(dp), dimension(:), allocatable :: runoff_obs
2313 :
2314 : ! mask for aggregated measured runoff
2315 0 : logical, dimension(:), allocatable :: runoff_obs_mask
2316 :
2317 :
2318 0 : call eval(parameterset, runoff=runoff)
2319 0 : nGaugesTotal = size(runoff, dim=2)
2320 :
2321 0 : objective_weighted_nse = 0.0_dp
2322 0 : do gg=1, nGaugesTotal
2323 : !
2324 0 : call extract_runoff( gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask )
2325 : !
2326 : objective_weighted_nse = objective_weighted_nse + &
2327 0 : wnse( runoff_obs, runoff_agg, mask=runoff_obs_mask)
2328 : end do
2329 : #ifndef MPI
2330 : ! objective_nse = objective_nse + nse(gauge%Q, runoff_model_agg) !, runoff_model_agg_mask)
2331 0 : objective_weighted_nse = 1.0_dp - objective_weighted_nse / real(nGaugesTotal,dp)
2332 :
2333 0 : call message('objective_weighted_nse (i.e., 1 - wNSE) = ', num2str(objective_weighted_nse))
2334 : ! pause
2335 : #endif
2336 :
2337 0 : deallocate( runoff_agg, runoff_obs, runoff_obs_mask )
2338 :
2339 0 : END FUNCTION objective_weighted_nse
2340 :
2341 : ! ------------------------------------------------------------------
2342 :
2343 : ! NAME
2344 : ! objective_sse_boxcox
2345 :
2346 : ! PURPOSE
2347 : !> \brief Objective function of sum of squared errors of transformed streamflow.
2348 :
2349 : !> \details The objective function only depends on a parameter vector.
2350 : !> The model will be called with that parameter vector and
2351 : !> the model output is subsequently compared to observed data.
2352 : !> Therefore, the sum of squared error \f$ tSSE \f$
2353 : !> \f[ tSSE = \sum_{i=1}^N (z(Q_{obs}(i), \lambda) - z(Q_{model}(i), \lambda))^2 \f]
2354 : !> is calculated where \f$ z \f$ is the transform and given by
2355 : !> \f[ z(x, \lambda) = \frac{x^\lambda -1}{\lambda} \f] for \f$ \lambda \f$ unequal to zero and
2356 : !> \f[ z(x, \lambda) = log x \f] for \f$ \lambda \f$ equal to zero.
2357 : !> The objective function is
2358 : !> \f[ obj\_value = tSSE \f]
2359 : !> The observed data \f$ Q_{obs} \f$ are global in this module.
2360 : !>
2361 : !> The boxcox transformation uses a parameter of 0.2, suggested by
2362 : !> Woldemeskel et al. Hydrol Earth Syst Sci, 2018 vol. 22 (12) pp. 6257-6278.
2363 : !> "Evaluating post-processing approaches for monthly and seasonal streamflow forecasts."
2364 : !> https://www.hydrol-earth-syst-sci.net/22/6257/2018/
2365 :
2366 :
2367 : ! INTENT(IN)
2368 : !> \param[in] "real(dp), dimension(:) :: parameterset"
2369 : !> \param[in] "procedure(eval_interface) :: eval"
2370 :
2371 : ! RETURN
2372 : !> \return real(dp) :: objective_sse_boxcox — objective function value
2373 : !> (which will be e.g. minimized by an optimization routine like DDS)
2374 :
2375 : ! HISTORY
2376 : !> \authors Stephan Thober, Dmitri Kavetski
2377 :
2378 : !> \date Aug 2019
2379 :
2380 0 : FUNCTION objective_sse_boxcox(parameterset, eval)
2381 :
2382 0 : use mo_errormeasures, only: SSE
2383 : use mo_boxcox, only: boxcox
2384 :
2385 : implicit none
2386 :
2387 : real(dp), dimension(:), intent(in) :: parameterset
2388 :
2389 : procedure(eval_interface), INTENT(IN), pointer :: eval
2390 :
2391 : real(dp) :: objective_sse_boxcox
2392 :
2393 : ! modelled runoff for a given parameter set
2394 : ! dim1=nTimeSteps, dim2=nGauges
2395 0 : real(dp), allocatable, dimension(:,:) :: runoff
2396 :
2397 : ! gauges counter
2398 : integer(i4) :: gg
2399 :
2400 : integer(i4) :: nGaugesTotal
2401 :
2402 : ! aggregated simulated runoff
2403 0 : real(dp), dimension(:), allocatable :: runoff_agg
2404 :
2405 : ! measured runoff
2406 0 : real(dp), dimension(:), allocatable :: runoff_obs
2407 :
2408 : ! mask for aggregated measured runoff
2409 0 : logical, dimension(:), allocatable :: runoff_obs_mask
2410 :
2411 : ! boxcox parameter
2412 0 : real(dp) :: lambda
2413 :
2414 : ! set box cox transformation to 0.2
2415 : ! suggested by:
2416 : ! Woldemeskel et al. Hydrol Earth Syst Sci, 2018 vol. 22 (12) pp. 6257-6278.
2417 : ! "Evaluating post-processing approaches for monthly and seasonal streamflow forecasts."
2418 : ! https://www.hydrol-earth-syst-sci.net/22/6257/2018/
2419 0 : lambda = 0.2_dp
2420 :
2421 0 : call eval(parameterset, runoff=runoff)
2422 0 : nGaugesTotal = size(runoff, dim=2)
2423 :
2424 0 : objective_sse_boxcox = 0.0_dp
2425 0 : do gg=1, nGaugesTotal
2426 : !
2427 0 : call extract_runoff( gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask )
2428 : !
2429 : objective_sse_boxcox = objective_sse_boxcox + &
2430 0 : SSE( boxcox(runoff_obs, lambda, mask=runoff_obs_mask), boxcox(runoff_agg, lambda), mask=runoff_obs_mask)
2431 : end do
2432 : #ifndef MPI
2433 : ! objective_nse = objective_nse + nse(gauge%Q, runoff_model_agg) !, runoff_model_agg_mask)
2434 0 : objective_sse_boxcox = objective_sse_boxcox / real(nGaugesTotal,dp)
2435 :
2436 0 : call message('objective_sse_boxcox = ', num2str(objective_sse_boxcox))
2437 : ! pause
2438 : #endif
2439 :
2440 0 : deallocate( runoff_agg, runoff_obs, runoff_obs_mask )
2441 :
2442 0 : END FUNCTION objective_sse_boxcox
2443 :
2444 :
2445 : ! ------------------------------------------------------------------
2446 :
2447 : ! NAME
2448 : ! extract_runoff
2449 :
2450 : ! PURPOSE
2451 : !> \brief extracts runoff data from global variables
2452 :
2453 : !> \details extracts simulated and measured runoff from global variables,
2454 : !> such that they overlay exactly. For measured runoff, only the runoff
2455 : !> during the evaluation period are cut, not succeeding nodata values.
2456 : !> For simulated runoff, warming days as well as succeeding nodata values
2457 : !> are neglected and the simulated runoff is aggregated to the resolution
2458 : !> of the observed runoff.
2459 : !> see use in this module above
2460 :
2461 : ! INTENT(IN)
2462 : !> \param[in] "integer(i4) :: gaugeId" current gauge Id
2463 : !> \param[in] "real(dp), dimension(:, :) :: runoff" simulated runoff
2464 :
2465 : ! INTENT(OUT)
2466 : !> \param[out] "real(dp), dimension(:) :: runoff_agg" aggregated simulated
2467 : !> \param[out] "real(dp), dimension(:) :: runoff_obs" extracted measured
2468 : !> \param[out] "logical, dimension(:) :: runoff_obs_mask" mask of no data values
2469 :
2470 : ! HISTORY
2471 : !> \authors Stephan Thober
2472 :
2473 : !> \date Jan 2015
2474 :
2475 : ! Modifications:
2476 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
2477 :
2478 25 : subroutine extract_runoff(gaugeId, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2479 :
2480 0 : use mo_common_mhm_mrm_variables, only : evalPer, nTstepDay, warmingDays
2481 : use mo_mrm_global_variables, only : gauge, nMeasPerDay
2482 : use mo_utils, only : ge
2483 :
2484 : implicit none
2485 :
2486 : ! current gauge Id
2487 : integer(i4), intent(in) :: gaugeId
2488 :
2489 : ! simulated runoff
2490 : real(dp), dimension(:, :), intent(in) :: runoff
2491 :
2492 : ! aggregated simulated
2493 : real(dp), dimension(:), allocatable, intent(out) :: runoff_agg
2494 :
2495 : ! extracted measured
2496 : real(dp), dimension(:), allocatable, intent(out) :: runoff_obs
2497 :
2498 : ! mask of no data values
2499 : logical, dimension(:), allocatable, intent(out) :: runoff_obs_mask
2500 :
2501 : ! domain id
2502 : integer(i4) :: iDomain
2503 :
2504 : ! timestep counter
2505 : integer(i4) :: tt
2506 :
2507 : ! length of extracted time series
2508 : integer(i4) :: length
2509 :
2510 : ! between simulated and measured time scale
2511 : integer(i4) :: factor
2512 :
2513 : ! simulated Timesteps per Day
2514 : integer(i4) :: TPD_sim
2515 :
2516 : ! observed Timesteps per Day
2517 : integer(i4) :: TPD_obs
2518 :
2519 25 : real(dp), dimension(:), allocatable :: dummy
2520 :
2521 :
2522 : ! copy time resolution to local variables
2523 25 : TPD_sim = nTstepDay
2524 25 : TPD_obs = nMeasPerDay
2525 :
2526 : ! check if modelled timestep is an integer multiple of measured timesteps
2527 25 : if (modulo(TPD_sim, TPD_obs) .eq. 0) then
2528 25 : factor = TPD_sim / TPD_obs
2529 : else
2530 0 : call error_message(' Error: Number of modelled datapoints is no multiple of measured datapoints per day')
2531 : end if
2532 :
2533 : ! extract domain Id from gauge Id
2534 25 : iDomain = gauge%domainId(gaugeId)
2535 :
2536 : ! get length of evaluation period times TPD_obs
2537 25 : length = (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1) * TPD_obs
2538 :
2539 : ! extract measurements
2540 25 : if (allocated(runoff_obs)) deallocate(runoff_obs)
2541 75 : allocate(runoff_obs(length))
2542 19037 : runoff_obs = gauge%Q(1 : length, gaugeId)
2543 :
2544 : ! create mask of observed runoff
2545 25 : if (allocated(runoff_obs_mask)) deallocate(runoff_obs_mask)
2546 75 : allocate(runoff_obs_mask(length))
2547 19012 : runoff_obs_mask = .false.
2548 19012 : forall(tt = 1 : length) runoff_obs_mask(tt) = ge(runoff_obs(tt), 0.0_dp)
2549 :
2550 : ! extract and aggregate simulated runoff
2551 25 : if (allocated(runoff_agg)) deallocate(runoff_agg)
2552 50 : allocate(runoff_agg(length))
2553 : ! remove warming days
2554 25 : length = (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1) * TPD_sim
2555 75 : allocate(dummy(length))
2556 455738 : dummy = runoff(warmingDays(iDomain) * TPD_sim + 1 : warmingDays(iDomain) * TPD_sim + length, gaugeId)
2557 : ! aggregate runoff
2558 25 : length = (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1) * TPD_obs
2559 0 : forall(tt = 1 : length) runoff_agg(tt) = sum(dummy((tt - 1) * factor + 1 : tt * factor)) / &
2560 474700 : real(factor, dp)
2561 : ! clean up
2562 25 : deallocate(dummy)
2563 :
2564 25 : end subroutine extract_runoff
2565 :
2566 :
2567 : END MODULE mo_mrm_objective_function_runoff
|
