5.13.2-dev0
mHM
The mesoscale Hydrological Model
Loading...
Searching...
No Matches
mo_objective_function.F90
Go to the documentation of this file.
1!> \file mo_objective_function.f90
2!> \brief \copybrief mo_objective_function
3!> \details \copydetails mo_objective_function
4
5! ToDo: change comment for OF 15
6
7!> \brief Objective Functions for Optimization of mHM.
8!> \details This module provides a wrapper for several objective functions used to optimize mHM against various
9!! variables.
10!! If the objective is only regarding runoff move it to mRM/mo_mrm_objective_function_runoff.f90.
11!! If it contains besides runoff another variable like TWS implement it here.
12!!
13!! All the objective functions are supposed to be minimized!
14!! - (10) SO: SM: 1.0 - KGE of catchment average soilmoisture
15!! - (11) SO: SM: 1.0 - Pattern dissimilarity (PD) of spatially distributed soil moisture
16!! - (12) SO: SM: Sum of squared errors (SSE) of spatially distributed standard score (normalization) of soil moisture
17!! - (13) SO: SM: 1.0 - average temporal correlation of spatially distributed soil moisture
18!! - (15) SO: Q + TWS: [1.0-KGE(Q)]*RMSE(domain_avg_TWS) - objective function using Q and domain average (standard score) TWS
19!! - (17) SO: N: 1.0 - KGE of spatio-temporal neutron data, catchment-average
20!! - (27) SO: ET: 1.0 - KGE of catchment average evapotranspiration
21!> \changelog
22!! - Oldrich Rakovec Oct 2015
23!! - added obj. func. 15 (objective_kge_q_rmse_tws) and extract_domain_avg_tws routine, former basin_avg
24!! - Robert Schweppe Jun 2018
25!! - refactoring and reformatting
26!> \authors Juliane Mai
27!> \date Dec 2012
28!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
29!! mHM is released under the LGPLv3+ license \license_note
30!> \ingroup f_mhm
31MODULE mo_objective_function
32
33 ! This module provides objective functions for optimization of the UFZ CHS mesoscale hydrologic model mHM.
34
35 ! Written Juliane Mai, Dec 2012
36 ! Modified Stephan Thober, Oct 2015 moved all runoff only related objectives to mRM
37
38 USE mo_kind, ONLY : i4, dp
39 use mo_optimization_utils, only : eval_interface
40 use mo_message, only : message, error_message
41
42 IMPLICIT NONE
43
44 PRIVATE
45
46 PUBLIC :: objective
47#ifdef MPI
48 PUBLIC :: objective_master, objective_subprocess ! objective function wrapper for soil moisture only
49#endif
50
51 ! ------------------------------------------------------------------
52
53CONTAINS
54
55 ! ------------------------------------------------------------------
56
57 ! NAME
58 ! objective
59
60 ! PURPOSE
61 !> \brief Wrapper for objective functions.
62
63 !> \details The functions selects the objective function case defined in a namelist,
64 !> i.e. the global variable \e opti\_function.
65 !> It return the objective function value for a specific parameter set.
66
67 ! INTENT(IN)
68 !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
69 !> \param[in] "procedure(eval_interface) :: eval"
70
71 ! INTENT(IN), OPTIONAL
72 !> \param[in] "real(dp), optional :: arg1"
73
74 ! INTENT(OUT), OPTIONAL
75 !> \param[out] "real(dp), optional :: arg2"
76 !> \param[out] "real(dp), optional :: arg3"
77
78 ! RETURN
79 !> \return real(dp) :: objective — objective function value
80 !> (which will be e.g. minimized by an optimization routine like DDS)
81
82 ! HISTORY
83 !> \authors Juliane Mai
84
85 !> \date Dec 2012
86
87 ! Modifications:
88 ! Stephan Thober Oct 2015 - moved all runoff related objective functions to mRM
89 ! Robert Schweppe Jun 2018 - refactoring and reformatting
90
91 FUNCTION objective(parameterset, eval, arg1, arg2, arg3)
92
93 use mo_common_constants, only : nodata_dp
94 use mo_common_mhm_mrm_variables, only : opti_function
95
96 implicit none
97
98 REAL(dp), DIMENSION(:), INTENT(IN) :: parameterset
99
100 procedure(eval_interface), INTENT(IN), POINTER :: eval
101
102 real(dp), optional, intent(in) :: arg1
103
104 real(dp), optional, intent(out) :: arg2
105
106 real(dp), optional, intent(out) :: arg3
107
108 REAL(dp) :: objective
109
110 real(dp), dimension(6) :: multiple_objective
111
112
113 if (present(arg1) .or. present(arg2) .or. present(arg3)) then
114 call error_message("Error mo_objective_function: Received unexpected argument, check optimization settings")
115 end if
116
117 ! set these to nan so compiler does not complain
118 if (present(arg2)) then
119 arg2 = nodata_dp
120 end if
121 if (present(arg3)) then
122 arg3 = nodata_dp
123 end if
124
125 select case (opti_function)
126 case (10)
127 ! KGE of catchment average SM
128 objective = objective_sm_kge_catchment_avg(parameterset, eval)
129 case (11)
130 ! pattern dissimilarity (PD) of SM fields
131 objective = objective_sm_pd(parameterset, eval)
132 case (12)
133 ! sum of squared errors of standard_score SM
134 objective = objective_sm_sse_standard_score(parameterset, eval)
135 case (13)
136 ! soil moisture correlation - temporal
137 objective = objective_sm_corr(parameterset, eval)
138 case (15)
139 ! KGE for Q * RMSE for domain_avg TWS (standarized scored)
140 objective = objective_kge_q_rmse_tws(parameterset, eval)
141 case (17)
142 ! KGE of catchment average SM
143 objective = objective_neutrons_kge_catchment_avg(parameterset, eval)
144 case (27)
145 ! KGE of catchment average ET
146 objective = objective_et_kge_catchment_avg(parameterset, eval)
147 case (28)
148 ! KGE for Q + SSE for SM (standarized scored)
149 objective = objective_kge_q_sm_corr(parameterset, eval)
150 case (29)
151 ! KGE for Q + KGE of catchment average ET
152 objective = objective_kge_q_et(parameterset, eval)
153 case (30)
154 ! KGE for Q * RMSE for domain_avg ET (standarized scored)
155 objective = objective_kge_q_rmse_et(parameterset, eval)
156 case (33)
157 multiple_objective = objective_q_et_tws_kge_catchment_avg(parameterset, eval)
158 objective = multiple_objective(1)
159 case (34)
160 ! KGE for Q * Absolute-Error for BFI
161 objective = objective_kge_q_bfi(parameterset, eval)
162
163 case default
164 call error_message("Error objective: opti_function not implemented yet.")
165 end select
166
167 END FUNCTION objective
168
169 ! ------------------------------------------------------------------
170
171 ! NAME
172 ! objective_master
173
174 ! PURPOSE
175 !> \brief Wrapper for objective functions.
176
177 !> \details The functions sends the parameterset to the subprocess, receives
178 !> the partial objective and calculates the final objective
179 ! INTENT(IN)
180 !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
181 !> \param[in] "procedure(eval_interface) :: eval"
182
183 ! INTENT(IN), OPTIONAL
184 !> \param[in] "real(dp), optional :: arg1"
185
186 ! INTENT(OUT), OPTIONAL
187 !> \param[out] "real(dp), optional :: arg2"
188 !> \param[out] "real(dp), optional :: arg3"
189
190 ! RETURN
191 !> \return real(dp) :: objective — objective function value
192 !> (which will be e.g. minimized by an optimization routine like DDS)
193
194 ! HISTORY
195 !> \authors Juliane Mai
196
197 !> \date Dec 2012
198
199 ! Modifications:
200 ! Stephan Thober Oct 2015 - moved all runoff related objective functions to mRM
201 ! Robert Schweppe Jun 2018 - refactoring and reformatting
202 ! Maren Kaluza Jun 2019 - parallel version
203
204#ifdef MPI
205 FUNCTION objective_master(parameterset, eval, arg1, arg2, arg3)
206
207 use mo_common_constants, only : nodata_dp
208 use mo_common_mhm_mrm_variables, only : opti_function
209 use mo_common_mpi_tools, only : distribute_parameterset
210 use mo_common_variables, only : domainmeta
211 use mo_string_utils, only : num2str
212 use mpi_f08
213
214 implicit none
215
216 REAL(dp), DIMENSION(:), INTENT(IN) :: parameterset
217
218 procedure(eval_interface), INTENT(IN), POINTER :: eval
219
220 real(dp), optional, intent(in) :: arg1
221
222 real(dp), optional, intent(out) :: arg2
223
224 real(dp), optional, intent(out) :: arg3
225
226 REAL(dp) :: objective_master
227
228 REAL(dp) :: partial_objective
229
230 real(dp), dimension(6) :: multiple_partial_objective
231
232 real(dp), dimension(6) :: multiple_master_objective
233
234 ! for sixth root
235 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
236
237 integer(i4) :: iproc, nproc
238
239 integer(i4) :: ierror
240
241 type(mpi_status) :: status
242
243
244 if (present(arg1) .or. present(arg2) .or. present(arg3)) then
245 call error_message("Error mo_objective_function: Received unexpected argument, check optimization settings")
246 end if
247
248 ! set these to nan so compiler does not complain
249 if (present(arg2)) then
250 arg2 = nodata_dp
251 end if
252 if (present(arg3)) then
253 arg3 = nodata_dp
254 end if
255
256 call distribute_parameterset(parameterset)
257 select case (opti_function)
258 case (10 : 13, 17, 27 : 29)
259 call mpi_comm_size(domainmeta%comMaster, nproc, ierror)
260 objective_master = 0.0_dp
261 do iproc = 1, nproc - 1
262 call mpi_recv(partial_objective, 1, mpi_double_precision, iproc, 0, domainmeta%comMaster, status, ierror)
263 objective_master = objective_master + partial_objective
264 end do
265 objective_master = objective_master**onesixth
266 case (15)
267 ! KGE for Q * RMSE for domain_avg TWS (standarized scored)
268 call error_message("case 15, objective_kge_q_rmse_tws not implemented in parallel yet")
269 case (30)
270 ! KGE for Q * RMSE for domain_avg ET (standarized scored)
271 ! objective_master = objective_kge_q_rmse_et(parameterset, eval)
272 call message("case 30, objective_kge_q_rmse_et not implemented in parallel yet")
273 case(33)
274 call mpi_comm_size(domainmeta%comMaster, nproc, ierror)
275 objective_master = 0.0_dp
276 multiple_master_objective(:) = 0.0_dp
277 do iproc = 1, nproc - 1
278 call mpi_recv(multiple_partial_objective, 6, mpi_double_precision, iproc, 0, domainmeta%comMaster, status, ierror)
279 multiple_master_objective = multiple_master_objective + multiple_partial_objective
280 end do
281 objective_master = objective_master + &
282 (multiple_master_objective(1)+multiple_master_objective(2)+multiple_master_objective(3))
283 objective_master = (objective_master/multiple_master_objective(4))**onesixth
284
285 case default
286 call error_message("Error objective_master: opti_function not implemented yet.")
287 end select
288
289 select case (opti_function)
290 case(10)
291 call message(' objective_sm_kge_catchment_avg = ', num2str(objective_master, '(F9.5)'))
292 case(11)
293 call message(' objective_sm_pd = ', num2str(objective_master, '(F9.5)'))
294 case(12)
295 call message(' objective_sm_sse_standard_score = ', num2str(objective_master, '(E12.5)'))
296 case(13)
297 call message(' objective_sm_corr = ', num2str(objective_master, '(F9.5)'))
298 case(17)
299 call message(' objective_neutrons_kge_catchment_avg = ', num2str(objective_master, '(F9.5)'))
300 case(27)
301 call message(' objective_et_kge_catchment_avg = ', num2str(objective_master, '(F9.5)'))
302 case(28)
303 call message(' objective_kge_q_sm_corr = ', num2str(objective_master, '(F9.5)'))
304 case(29)
305 call message(' objective_kge_q_et = ', num2str(objective_master, '(F9.5)'))
306 case(33)
307 call message(' objective_q_et_tws_kge_catchment_avg = ', num2str(objective_master, '(F9.5)'))
308 case default
309 call error_message("Error objective_master: opti_function not implemented yet, this part of the code should never execute.")
310 end select
311
312 END FUNCTION objective_master
313
314 ! ------------------------------------------------------------------
315
316 ! NAME
317 ! objective_subprocess
318
319 ! PURPOSE
320 !> \brief Wrapper for objective functions.
321
322 !> \details The function receives the parameterset from the master
323 !> process, selects the objective function case defined in a namelist,
324 !> i.e. the global variable \e opti\_function.
325 !> It returns the partial objective function value for a specific parameter set.
326
327 ! INTENT(IN)
328 !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
329 !> \param[in] "procedure(eval_interface) :: eval"
330
331 ! INTENT(IN), OPTIONAL
332 !> \param[in] "real(dp), optional :: arg1"
333
334 ! INTENT(OUT), OPTIONAL
335 !> \param[out] "real(dp), optional :: arg2"
336 !> \param[out] "real(dp), optional :: arg3"
337
338 ! RETURN
339 !> \return real(dp) :: objective — objective function value
340 !> (which will be e.g. minimized by an optimization routine like DDS)
341
342 ! HISTORY
343 !> \authors Juliane Mai
344
345 !> \date Dec 2012
346
347 ! Modifications:
348 ! Stephan Thober Oct 2015 - moved all runoff related objective functions to mRM
349 ! Robert Schweppe Jun 2018 - refactoring and reformatting
350 ! Maren Kaluza Jun 2019 - parallel version
351
352 subroutine objective_subprocess(eval, arg1, arg2, arg3)
353
354 use mo_common_constants, only : nodata_dp
355 use mo_common_mhm_mrm_variables, only : opti_function
356 use mo_common_mpi_tools, only : get_parameterset
357 use mo_common_variables, only : domainmeta
358 use mpi_f08
359
360 implicit none
361
362 procedure(eval_interface), INTENT(IN), POINTER :: eval
363
364 real(dp), optional, intent(in) :: arg1
365
366 real(dp), optional, intent(out) :: arg2
367
368 real(dp), optional, intent(out) :: arg3
369
370 REAL(dp) :: partial_objective
371
372 real(dp), dimension(6) :: multiple_partial_objective
373
374 REAL(dp), DIMENSION(:), allocatable :: parameterset
375
376 integer(i4) :: ierror
377
378 type(mpi_status) :: status
379
380 logical :: do_obj_loop
381
382 do ! a do loop without condition runs until exit
383 call mpi_recv(do_obj_loop, 1, mpi_logical, 0, 0, domainmeta%comMaster, status, ierror)
384
385 if (.not. do_obj_loop) exit
386
387 if (present(arg1) .or. present(arg2) .or. present(arg3)) then
388 call error_message("Error mo_objective_function: Received unexpected argument, check optimization settings")
389 end if
390
391 ! set these to nan so compiler does not complain
392 if (present(arg2)) then
393 arg2 = nodata_dp
394 end if
395 if (present(arg3)) then
396 arg3 = nodata_dp
397 end if
398 call get_parameterset(parameterset)
399 select case (opti_function)
400 case (10)
401 ! KGE of catchment average SM
402 partial_objective = objective_sm_kge_catchment_avg(parameterset, eval)
403 case (11)
404 ! pattern dissimilarity (PD) of SM fields
405 partial_objective = objective_sm_pd(parameterset, eval)
406 case (12)
407 ! sum of squared errors of standard_score SM
408 partial_objective = objective_sm_sse_standard_score(parameterset, eval)
409 case (13)
410 ! soil moisture correlation - temporal
411 partial_objective = objective_sm_corr(parameterset, eval)
412 case (15)
413 ! KGE for Q * RMSE for domain_avg TWS (standarized scored)
414 ! partial_objective = objective_kge_q_rmse_tws(parameterset, eval)
415 call error_message("Error objective_subprocess: case 15 not supported with MPI.")
416 case (17)
417 ! KGE of catchment average SM
418 partial_objective = objective_neutrons_kge_catchment_avg(parameterset, eval)
419 case (27)
420 ! KGE of catchment average ET
421 partial_objective = objective_et_kge_catchment_avg(parameterset, eval)
422 case (28)
423 ! KGE for Q + SSE for SM (standarized scored)
424 partial_objective = objective_kge_q_sm_corr(parameterset, eval)
425 case (29)
426 ! KGE for Q + KGE of catchment average ET
427 partial_objective = objective_kge_q_et(parameterset, eval)
428 case (30)
429 ! KGE for Q * RMSE for domain_avg ET (standarized scored)
430 ! partial_objective = objective_kge_q_rmse_et(parameterset, eval)
431 call error_message("Error objective_subprocess: case 30 not supported with MPI.")
432 case(33)
433 multiple_partial_objective = objective_q_et_tws_kge_catchment_avg(parameterset, eval)
434 case default
435 call error_message("Error objective_subprocess: opti_function not implemented yet.")
436 end select
437
438 select case (opti_function)
439 case (10 : 13, 17, 27 : 29)
440 call mpi_send(partial_objective,1, mpi_double_precision,0,0,domainmeta%comMaster,ierror)
441 case(33)
442 call mpi_send(multiple_partial_objective, 6, mpi_double_precision,0,0,domainmeta%comMaster,ierror)
443 case default
444 call error_message("Error objective_subprocess: this part should not be executed -> error in the code.")
445 end select
446
447 deallocate(parameterset)
448 end do
449
450 END subroutine objective_subprocess
451
452#endif
453 ! ------------------------------------------------------------------
454
455 ! NAME
456 ! objective_sm_kge_catchment_avg
457
458 ! PURPOSE
459 !> \brief Objective function for soil moisture.
460
461 !> \details The objective function only depends on a parameter vector.
462 !> The model will be called with that parameter vector and
463 !> the model output is subsequently compared to observed data.
464
465 !> Therefore, the Kling-Gupta model efficiency \f$ KGE \f$ of the catchment average
466 !> soil mloisture (SM) is calculated
467 !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f]
468 !> where
469 !> \f$ r \f$ = Pearson product-moment correlation coefficient
470 !> \f$ \alpha \f$ = ratio of simulated mean to observed mean SM
471 !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation
472 !> is calculated and the objective function for a given domain \f$ i \f$ is
473 !> \f[ \phi_{i} = 1.0 - KGE_{i} \f]
474 !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods.
475 !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.
476
477 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
478 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
479 !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f]
480 !> The observed data L1_sm, L1_sm_mask are global in this module.
481
482 ! INTENT(IN)
483 !> \param[in] "real(dp), dimension(:) :: parameterset"
484 !> \param[in] "procedure(eval_interface) :: eval"
485
486 ! RETURN
487 !> \return real(dp) :: objective_sm_kge_catchment_avg — objective function value
488 !> (which will be e.g. minimized by an optimization routine like DDS)
489
490 ! HISTORY
491 !> \authors Matthias Zink
492
493 !> \date May 2015
494
495 ! Modifications:
496 ! Robert Schweppe Jun 2018 - refactoring and reformatting
497
498 FUNCTION objective_sm_kge_catchment_avg(parameterset, eval)
499
500 use mo_optimization_types, only : optidata_sim
501 use mo_common_constants, only : nodata_dp
502 use mo_common_variables, only : level1, domainmeta
503 use mo_errormeasures, only : kge
504 use mo_global_variables, only : l1_smobs
505 use mo_moment, only : average
506 use mo_string_utils, only : num2str
507
508 implicit none
509
510 real(dp), dimension(:), intent(in) :: parameterset
511
512 procedure(eval_interface), INTENT(IN), POINTER :: eval
513
514 real(dp) :: objective_sm_kge_catchment_avg
515
516 ! domain loop counter
517 integer(i4) :: idomain
518
519 ! time loop counter
520 integer(i4) :: itime
521
522 ! number of time steps in simulated SM
523 integer(i4) :: n_time_steps
524
525 ! ncells1 of level 1
526 integer(i4) :: ncells1
527
528 ! number of invalid timesteps
529 real(dp) :: invalid_times
530#ifndef MPI
531 ! for sixth root
532 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
533#endif
534
535 ! spatial average of observed soil moisture
536 real(dp), dimension(:), allocatable :: sm_catch_avg_domain
537
538 ! spatial avergae of modeled soil moisture
539 real(dp), dimension(:), allocatable :: sm_opti_catch_avg_domain
540
541 type(optidata_sim), dimension(:), allocatable :: smoptisim
542
543 ! mask for valid sm catchment avg time steps
544 logical, dimension(:), allocatable :: mask_times
545
546
547 allocate(smoptisim(domainmeta%nDomains))
548 call eval(parameterset, smoptisim = smoptisim)
549
550 ! initialize some variables
551 objective_sm_kge_catchment_avg = 0.0_dp
552
553 ! loop over domain - for applying power law later on
554 do idomain = 1, domainmeta%nDomains
555
556 ! get domain information
557 ncells1 = level1(idomain)%ncells
558
559 ! allocate
560 allocate(mask_times(size(smoptisim(idomain)%dataSim, dim = 2)))
561 allocate(sm_catch_avg_domain(size(smoptisim(idomain)%dataSim, dim = 2)))
562 allocate(sm_opti_catch_avg_domain(size(smoptisim(idomain)%dataSim, dim = 2)))
563
564 ! initalize
565 mask_times = .true.
566 sm_catch_avg_domain = nodata_dp
567 sm_opti_catch_avg_domain = nodata_dp
568
569 invalid_times = 0.0_dp
570 ! calculate catchment average soil moisture
571 n_time_steps = size(smoptisim(idomain)%dataSim, dim = 2)
572 do itime = 1, n_time_steps
573
574 ! check for enough data points in timesteps for KGE calculation
575 ! more then 10 percent avaiable in current field
576 if (count(l1_smobs(idomain)%maskObs(:, itime)) .LE. (0.10_dp * real(ncells1, dp))) then
577 invalid_times = invalid_times + 1.0_dp
578 mask_times(itime) = .false.
579 cycle
580 end if
581 sm_catch_avg_domain(itime) = average(l1_smobs(idomain)%dataObs(:, itime), mask = l1_smobs(idomain)%maskObs(:, itime))
582 sm_opti_catch_avg_domain(itime) = average(smoptisim(idomain)%dataSim(:, itime), mask = l1_smobs(idomain)%maskObs(:, itime))
583 end do
584
585 ! user information about invalid times
586 if (invalid_times .GT. 0.5_dp) then
587 call message(.LT.' WARNING: objective_sm: Detected invalid timesteps ( 10 valid data points).')
588 call message(' Fraction of invalid timesteps: ', &
589 num2str(invalid_times / real(n_time_steps, dp), '(F4.2)'))
590 end if
591
592
593 ! calculate average soil moisture KGE over all domains with power law
594 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
595 objective_sm_kge_catchment_avg = objective_sm_kge_catchment_avg + &
596 ((1.0_dp - kge(sm_catch_avg_domain, sm_opti_catch_avg_domain, mask = mask_times)) / &
597 real(domainmeta%overallnumberofdomains, dp))**6
598
599 ! deallocate
600 deallocate(mask_times)
601 deallocate(sm_catch_avg_domain)
602 deallocate(sm_opti_catch_avg_domain)
603 call smoptisim(idomain)%destroy()
604 end do
605 deallocate(smoptisim)
606
607#ifndef MPI
608 objective_sm_kge_catchment_avg = objective_sm_kge_catchment_avg**onesixth
609
610 call message(' objective_sm_kge_catchment_avg = ', num2str(objective_sm_kge_catchment_avg, '(F9.5)'))
611#endif
612
613
614 END FUNCTION objective_sm_kge_catchment_avg
615
616 ! ------------------------------------------------------------------
617
618 ! NAME
619 ! objective_q_et_tws_kge_catchment_avg
620
621 ! PURPOSE
622 !> \brief Objective function for et, tws and q.
623
624 !> \details The feature of this objective function is the
625 !> separation of the eval call into four
626 !> calls, each with another index list. The subroutine eval then only
627 !> uses the indices from that index list internally instead of having loops
628 !> over all domains. The integer array domainMeta%optidata decides which
629 !> indices to use. Therefore the array is split into disjunct subsets, and,
630 !> if chosen wisely in the namelist, also covers all domains.
631 !>
632 !> With this the eval calls sum up in a way that for each domain eval was
633 !> called at most once, but for different opti_data.
634
635 ! HISTORY
636 !> \authors Maren Kaluza
637
638 !> \date July 2019
639
640 ! Modifications:
641
642 FUNCTION objective_q_et_tws_kge_catchment_avg(parameterset, eval)
643
644 use mo_optimization_types, only : optidata_sim
645 use mo_common_constants, only : nodata_dp
646 use mo_common_variables, only : domainmeta
647 use mo_global_variables, only : l1_etobs, l1_twsaobs
648 use mo_errormeasures, only : kge
649 use mo_moment, only : average
650 use mo_string_utils, only : num2str
651 use mo_mrm_objective_function_runoff, only : extract_runoff
652
653 implicit none
654
655 !> the parameterset passed to the eval subroutine
656 real(dp), dimension(:), intent(in) :: parameterset
657 !> the eval subroutine called by this objective function
658 procedure(eval_interface), INTENT(IN), POINTER :: eval
659 !> the return value of the objective function. In this case it is
660 !> an array to provide the possibility to weight the outcome accordingly
661 real(dp), dimension(6) :: objective_q_et_tws_kge_catchment_avg
662
663 !> domain loop counter
664 integer(i4) :: idomain
665
666 !> counter for short loops
667 integer(i4) :: i
668#ifndef MPI
669 !> for sixth root
670 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
671#endif
672
673 !> modelled runoff for a given parameter set
674 ! dim1=nTimeSteps, dim2=nGauges
675 real(dp), allocatable, dimension(:, :) :: runoff
676 !> number of all gauges, aquired via runoff
677 integer(i4) :: ngaugestotal
678
679 !> aggregated simulated runoff
680 real(dp), dimension(:), allocatable :: runoff_agg
681
682 !> measured runoff
683 real(dp), dimension(:), allocatable :: runoff_obs
684
685 !> mask for measured runoff
686 logical, dimension(:), allocatable :: runoff_obs_mask
687
688 !> kge_q(nGaugesTotal)
689 real(dp) :: kge_q
690
691 !> gauges counter
692 integer(i4) :: gg, icell
693
694 !> number of q domains
695 integer(i4) :: nqdomains
696
697 !> number of et domains
698 integer(i4) :: netdomains
699
700 !> number of tws domains
701 integer(i4) :: ntwsdomains
702
703 !> number of TWS and ET domains (providing both)
704 integer(i4) :: nettwsdomains
705
706 !> index array of ET domains
707 integer(i4), dimension(:), allocatable :: opti_domain_indices_et
708
709 !> index array of TWS domains
710 integer(i4), dimension(:), allocatable :: opti_domain_indices_tws
711
712 !> index array of TWS and ET domains (providing both)
713 integer(i4), dimension(:), allocatable :: opti_domain_indices_et_tws
714
715 !> index array of ET domains
716 integer(i4), dimension(:), allocatable :: opti_domain_indices_q
717
718 !> simulated et
719 type(optidata_sim), dimension(:), allocatable :: etoptisim
720
721 !> simulated tws
722 type(optidata_sim), dimension(:), allocatable :: twsoptisim
723
724 !> simulated twsa (anomaly)
725 type(optidata_sim), dimension(:), allocatable :: twsaoptisim
726
727 real(dp) :: kge_tws
728
729 real(dp) :: kge_et
730
731 integer(i4) :: numberofsummands
732
733
734 ! initialize some variables
735 objective_q_et_tws_kge_catchment_avg(:) = 0.0_dp
736 kge_tws = 0.0_dp
737 kge_et = 0.0_dp
738 kge_q = 0.0_dp
739 numberofsummands = 0
740 !--------------------------------------------
741 ! ET & TWS
742 !--------------------------------------------
743 ! eval runs to get simulated output for et and tws
744 ! before each eval call we generate an index list of the domains for which
745 ! eval should be called. Read details for further information
746 call init_indexarray_for_opti_data(domainmeta, 6, nettwsdomains, opti_domain_indices_et_tws)
747 if (nettwsdomains > 0) then
748 allocate( etoptisim(domainmeta%nDomains))
749 allocate(twsoptisim(domainmeta%nDomains))
750 allocate(twsaoptisim(domainmeta%nDomains))
751 call eval(parameterset, opti_domain_indices = opti_domain_indices_et_tws, &
752 twsoptisim = twsoptisim, etoptisim = etoptisim)
753 ! for all domains that have ET and TWS
754 do i = 1, size(opti_domain_indices_et_tws)
755 idomain = opti_domain_indices_et_tws(i)
756 call convert_tws_to_twsa(twsoptisim(idomain), l1_twsaobs(idomain), twsaoptisim(idomain))
757 do icell = 1, size(l1_etobs(idomain)%maskObs(:, :), dim = 1)
758 kge_et = kge_et + &
759 (1.0_dp - kge(l1_etobs(idomain)%dataObs(icell, :), etoptisim(idomain)%dataSim(icell, :),&
760 mask = l1_etobs(idomain)%maskObs(icell, :)))**6
761 numberofsummands = numberofsummands + 1
762 end do
763 do icell = 1, size(l1_twsaobs(idomain)%maskObs(:, :), dim = 1)
764 kge_tws = kge_tws + &
765 (1.0_dp - kge(l1_twsaobs(idomain)%dataObs(icell, :), twsaoptisim(idomain)%dataSim(icell, :),&
766 mask = l1_twsaobs(idomain)%maskObs(icell, :)))**6
767 numberofsummands = numberofsummands + 1
768 end do
769 ! deallocate
770 call etoptisim(idomain)%destroy()
771 call twsoptisim(idomain)%destroy()
772 call twsaoptisim(idomain)%destroy()
773 end do
774 deallocate(etoptisim)
775 deallocate(twsoptisim)
776 deallocate(twsaoptisim)
777 ! write(0,*) 'nEtTwsDomains, kge_tws', nEtTwsDomains, kge_tws
778 ! write(0,*) 'nEtTwsDomains, kge_et', nEtTwsDomains, kge_et
779 end if
780 !--------------------------------------------
781 ! TWS
782 !--------------------------------------------
783 ! eval runs to get simulated output for tws
784 ! before each eval call we generate an index list of the domains for which
785 ! eval should be called. Read details for further information
786 call init_indexarray_for_opti_data(domainmeta, 3, ntwsdomains, opti_domain_indices_tws)
787 if (ntwsdomains > 0) then
788 allocate(twsoptisim(domainmeta%nDomains))
789 allocate(twsaoptisim(domainmeta%nDomains))
790 call eval(parameterset, opti_domain_indices = opti_domain_indices_tws, twsoptisim = twsoptisim)
791 ! for all domains that have ET and TWS
792 do i = 1, size(opti_domain_indices_tws)
793 idomain = opti_domain_indices_tws(i)
794 call convert_tws_to_twsa(twsoptisim(idomain), l1_twsaobs(idomain), twsaoptisim(idomain))
795 do icell = 1, size(l1_twsaobs(idomain)%maskObs(:, :), dim = 1)
796 kge_tws = kge_tws + &
797 (1.0_dp - kge(l1_twsaobs(idomain)%dataObs(icell, :), twsaoptisim(idomain)%dataSim(icell, :),&
798 mask = l1_twsaobs(idomain)%maskObs(icell, :)))**6
799 numberofsummands = numberofsummands + 1
800 end do
801 call twsoptisim(idomain)%destroy()
802 end do
803 deallocate(twsoptisim)
804 ! write(0,*) 'nTwsDomains, kge_tws', nTwsDomains, kge_tws
805 end if
806 objective_q_et_tws_kge_catchment_avg(2) = kge_tws
807
808 !--------------------------------------------
809 ! ET
810 !--------------------------------------------
811 ! eval runs to get simulated output for et
812 ! before each eval call we generate an index list of the domains for which
813 ! eval should be called. Read details for further information
814 call init_indexarray_for_opti_data(domainmeta, 5, netdomains, opti_domain_indices_et)
815 if (netdomains > 0) then
816 allocate(etoptisim(domainmeta%nDomains))
817 call eval(parameterset, opti_domain_indices = opti_domain_indices_et, etoptisim = etoptisim)
818 ! for all domains that have ET and TWS
819 do i = 1, size(opti_domain_indices_et)
820 idomain = opti_domain_indices_et(i)
821 do icell = 1, size(l1_etobs(idomain)%maskObs(:, :), dim = 1)
822 kge_et = kge_et + &
823 (1.0_dp - kge(l1_etobs(idomain)%dataObs(icell, :), etoptisim(idomain)%dataSim(icell, :),&
824 mask = l1_etobs(idomain)%maskObs(icell, :)))**6
825 numberofsummands = numberofsummands + 1
826 end do
827 call etoptisim(idomain)%destroy()
828 end do
829 deallocate(etoptisim)
830 ! write(0,*) 'nEtDomains, kge_et', nEtDomains, kge_et
831 end if
832 objective_q_et_tws_kge_catchment_avg(3) = kge_et
833
834 !--------------------------------------------
835 ! RUNOFF
836 !--------------------------------------------
837 ! eval runs to get simulated output for runoff
838 ! before the eval call we generate an index list of the domains for which
839 ! eval should be called. Read details for further information
840 ! ToDo: The arrays for qTin, qTout, will be rewritten in the other calls when
841 ! Q is not called last. Change that for more flexibility
842 call init_indexarray_for_opti_data(domainmeta, 1, nqdomains, opti_domain_indices_q)
843
844 if (nqdomains > 0) then
845 call eval(parameterset, opti_domain_indices = opti_domain_indices_q, runoff = runoff)
846 ngaugestotal = size(runoff, dim = 2)
847 do gg = 1, ngaugestotal
848
849 ! extract runoff
850 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
851
852 kge_q = kge_q + &
853 (1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask))**6
854 numberofsummands = numberofsummands + 1
855 deallocate (runoff_agg, runoff_obs, runoff_obs_mask)
856 end do
857 ! write(0,*) 'nQDomains, kge_q', nQDomains, kge_q
858 end if
859 objective_q_et_tws_kge_catchment_avg(1) = kge_q
860
861 objective_q_et_tws_kge_catchment_avg(4) = real(numberofsummands, dp)
862
863
864#ifndef MPI
865 objective_q_et_tws_kge_catchment_avg(1) = ((kge_q+kge_et+kge_tws)/real(numberofsummands, dp))**onesixth
866
867 call message(' objective_q_et_tws_kge_catchment_avg = ', &
868 num2str(objective_q_et_tws_kge_catchment_avg(1), '(F9.5)'))
869#endif
870
871 END FUNCTION objective_q_et_tws_kge_catchment_avg
872
873 ! ------------------------------------------------------------------
874
875 ! NAME
876 ! init_indexarray_for_opti_data
877
878 ! PURPOSE
879 !> \brief creates an index array of the inidices of the domains eval
880 !> should MPI process.
881 !
882 !> \details The data type domainMeta contains an array optidata of size
883 !> domainMeta%nDomains, telling us, which domains should be
884 !> optimized with which opti_data. This subroutine splits all
885 !> domains assigned to a process and returns an index list
886 !> corresponding to the value of domainMeta%optidata.
887 !>
888 !> The index array opti_domain_indices can then be passed
889 !> as an optional argument to the eval subroutine. The
890 !> eval then instead of using loops over all domains only
891 !> uses the passed indices.
892 !>
893 !> This subroutine also returns the size of that array since it
894 !> helps with the calculations of the optimization in the end.
895
896 ! HISTORY
897 !> \authors Maren Kaluza
898
899 !> \date July 2019
900 subroutine init_indexarray_for_opti_data(domainMeta, optidataOption, nOptiDomains, opti_domain_indices)
901 use mo_common_types, only: domain_meta
902 !> meta data for all domains assigned to that process
903 type(domain_meta), intent(in) :: domainMeta
904 !> which opti data should be used in the eval called after calling this subroutine
905 integer(i4), intent(in) :: optidataOption
906 !> number of domains that will be optimized in the following eval call
907 integer(i4), intent(out) :: nOptiDomains
908 !> the indices of the domains that are to be processed in the following eval call
909 integer(i4), dimension(:), allocatable, intent(out) :: opti_domain_indices
910
911 !> domain loop counter
912 integer(i4) :: iDomain, i
913
914 if (allocated(opti_domain_indices)) deallocate(opti_domain_indices)
915 ! count domains on MPI process that use optidata
916 noptidomains = 0
917 do idomain = 1, domainmeta%nDomains
918 if (domainmeta%optidata(idomain) == optidataoption) noptidomains = noptidomains + 1
919 end do
920 ! write indices of these domains into an array
921 if (noptidomains > 0) then
922 allocate(opti_domain_indices(noptidomains))
923 i = 0
924 do idomain = 1, domainmeta%nDomains
925 if (domainmeta%optidata(idomain) == optidataoption) then
926 i = i + 1
927 opti_domain_indices(i) = idomain
928 end if
929 end do
930 end if
931 end subroutine init_indexarray_for_opti_data
932 ! ------------------------------------------------------------------
933
934 ! NAME
935 ! objective_sm_corr
936
937 ! PURPOSE
938 !> \brief Objective function for soil moisture.
939
940 !> \details The objective function only depends on a parameter vector.
941 !> The model will be called with that parameter vector and
942 !> the model output is subsequently compared to observed data.
943
944 !> Therefore the Pearson correlation between observed and modeled soil
945 !> moisture on each grid cell \f$ j \f$ is compared
946 !> \f[ r_j = r^2(SM_{obs}^j, SM_{sim}^j) \f]
947 !> where
948 !> \f$ r^2\f$ = Pearson correlation coefficient,
949 !> \f$ SM_{obs} \f$ = observed soil moisture,
950 !> \f$ SM_{sim} \f$ = simulated soil moisture.
951 !> The observed data \f$ SM_{obs} \f$ are global in this module.
952
953 !> The the correlation is spatially averaged as
954 !> \f[ \phi_{i} = \frac{1}{K} \cdot \sum_{j=1}^K r_j \f]
955 !> where \f$ K \f$ denotes the number of valid cells in the study domain.
956 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
957 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
958 !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 }. \f]
959 !> The observed data L1_sm, L1_sm_mask are global in this module.
960
961 ! INTENT(IN)
962 !> \param[in] "real(dp), dimension(:) :: parameterset"
963 !> \param[in] "procedure(eval_interface) :: eval"
964
965 ! RETURN
966 !> \return real(dp) :: objective_sm_corr — objective function value
967 !> (which will be e.g. minimized by an optimization routine like DDS)
968
969 ! HISTORY
970 !> \authors Matthias Zink
971
972 !> \date March 2015
973
974 ! Modifications:
975 ! Robert Schweppe Jun 2018 - refactoring and reformatting
976
977 FUNCTION objective_sm_corr(parameterset, eval)
978
979 use mo_optimization_types, only : optidata_sim
980 use mo_common_variables, only : level1, domainmeta
981 use mo_global_variables, only : l1_smobs
982 use mo_moment, only : correlation
983 use mo_string_utils, only : num2str
984
985 implicit none
986
987 real(dp), dimension(:), intent(in) :: parameterset
988
989 procedure(eval_interface), INTENT(IN), POINTER :: eval
990
991 real(dp) :: objective_sm_corr
992
993 ! domain loop counter
994 integer(i4) :: idomain
995
996 ! cell loop counter
997 integer(i4) :: icell
998
999 ! ncells1 of level 1
1000 integer(i4) :: ncells1
1001
1002 ! number of invalid cells in catchment
1003 real(dp) :: invalid_cells
1004
1005 ! domains wise objectives
1006 real(dp) :: objective_sm_corr_domain
1007
1008#ifndef MPI
1009 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1010#endif
1011
1012 type(optidata_sim), dimension(:), allocatable :: smoptisim
1013
1014
1015 allocate(smoptisim(domainmeta%nDomains))
1016 call eval(parameterset, smoptisim = smoptisim)
1017
1018 ! initialize some variables
1019 objective_sm_corr = 0.0_dp
1020
1021 ! loop over domain - for applying power law later on
1022 do idomain = 1, domainmeta%nDomains
1023
1024 ! init
1025 objective_sm_corr_domain = 0.0_dp
1026 ! get domain information
1027 ncells1 = level1(idomain)%ncells
1028
1029 invalid_cells = 0.0_dp
1030 ! temporal correlation is calculated on individual gridd cells
1031
1032 do icell = 1, size(l1_smobs(idomain)%maskObs(:, :), dim = 1)
1033
1034 ! check for enough data points in time for correlation
1035 if (count(l1_smobs(idomain)%maskObs(icell, :)) .LE. 0.10_dp * real(size(l1_smobs(idomain)%dataObs(:, :), dim = 2), dp)) then
1036 invalid_cells = invalid_cells + 1.0_dp
1037 cycle
1038 end if
1039 objective_sm_corr_domain = objective_sm_corr_domain + &
1040 correlation(l1_smobs(idomain)%dataObs(icell, :), smoptisim(idomain)%dataSim(icell, :), &
1041 mask = l1_smobs(idomain)%maskObs(icell, :))
1042 end do
1043
1044 ! user information about invalid cells
1045 if (invalid_cells .GT. 0.5_dp) then
1046 call message(.LT.' WARNING: objective_sm: Detected invalid cells in study area ( 10 valid data points).')
1047 call message(' Fraction of invalid cells: ', &
1048 num2str(invalid_cells / real(ncells1, dp), '(F4.2)'))
1049 end if
1050
1051
1052 ! calculate average soil moisture correlation over all domains with power law
1053 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1054 objective_sm_corr = objective_sm_corr + &
1055 ((1.0_dp - objective_sm_corr_domain / real(ncells1, dp)) / real(domainmeta%overallNumberOfDomains, dp))**6
1056 end do
1057#ifndef MPI
1058 objective_sm_corr = objective_sm_corr**onesixth
1059
1060 call message(' objective_sm_corr = ', num2str(objective_sm_corr, '(F9.5)'))
1061#endif
1062
1063 END FUNCTION objective_sm_corr
1064
1065 ! ------------------------------------------------------------------
1066
1067 ! NAME
1068 ! objective_sm_pd
1069
1070 ! PURPOSE
1071 !> \brief Objective function for soil moisture.
1072
1073 !> \details The objective function only depends on a parameter vector.
1074 !> The model will be called with that parameter vector and
1075 !> the model output is subsequently compared to observed data.
1076
1077 !> Therefore the Pattern Dissimilarity (PD) of observed and modeled soil
1078 !> moisture fields is calculated - aim: matching spatial patters
1079 !> \f[ E(t) = PD\left( SM_{obs}(t), SM_{sim}(t) \right) \f]
1080 !> where
1081 !> \f$ PD \f$ = pattern dissimilarity function,
1082 !> \f$ SM_{obs} \f$ = observed soil moisture,
1083 !> \f$ SM_{sim} \f$ = simulated soil moisture.
1084 !> \f$ E(t) \f$ = pattern dissimilarity at timestep \f$ t \f$.
1085 !> The the pattern dissimilaity (E) is spatially averaged as
1086 !> \f[ \phi_{i} = \frac{1}{T} \cdot \sum_{t=1}^T E_t \f]
1087 !> where \f$ T \f$ denotes the number of time steps.
1088 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1089 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1090 !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 } . \f]
1091 !> The observed data L1_sm, L1_sm_mask are global in this module.
1092 !> The observed data L1_sm, L1_sm_mask are global in this module.
1093
1094 ! INTENT(IN)
1095 !> \param[in] "real(dp), dimension(:) :: parameterset"
1096 !> \param[in] "procedure(eval_interface) :: eval"
1097
1098 ! RETURN
1099 !> \return real(dp) :: objecive_sm_pd — objective function value
1100 !> (which will be e.g. minimized by an optimization routine like DDS)
1101
1102 ! HISTORY
1103 !> \authors Matthias Zink
1104
1105 !> \date May 2015
1106
1107 ! Modifications:
1108 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1109
1110 FUNCTION objective_sm_pd(parameterset, eval)
1111
1112 use mo_optimization_types, only : optidata_sim
1113 use mo_common_constants, only : nodata_dp
1114 use mo_common_variables, only : level1, domainmeta
1115 use mo_global_variables, only : l1_smobs
1116 use mo_spatialsimilarity, only : pd
1117 use mo_string_utils, only : num2str
1118
1119 implicit none
1120
1121 real(dp), dimension(:), intent(in) :: parameterset
1122
1123 procedure(eval_interface), INTENT(IN), POINTER :: eval
1124
1125 ! objective function value
1126 real(dp) :: objective_sm_pd
1127
1128 ! domain loop counter
1129 integer(i4) :: idomain
1130
1131 ! time loop counter
1132 integer(i4) :: itime
1133
1134 ! level 1 number of culomns and rows
1135 integer(i4) :: nrows1, ncols1
1136
1137 ! for sixth root
1138#ifndef MPI
1139 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1140#endif
1141
1142 ! matrices of SM from vectorized arrays
1143 real(dp), dimension(:, :), allocatable :: mat1, mat2
1144
1145 ! pattern dissimilarity (pd) at every time step
1146 real(dp), dimension(:), allocatable :: pd_time_series
1147
1148 ! simulated soil moisture
1149 type(optidata_sim), dimension(:), allocatable :: smoptisim
1150
1151 ! mask of valid cells at level1
1152 logical, dimension(:, :), allocatable :: mask1
1153
1154 ! mask of valid sm cells
1155 logical, dimension(:, :), allocatable :: mask_sm
1156
1157 ! mask for valid sm catchment avg time steps
1158 logical, dimension(:), allocatable :: mask_times
1159
1160
1161 allocate(smoptisim(domainmeta%nDomains))
1162 call eval(parameterset, smoptisim = smoptisim)
1163
1164 ! initialize some variables
1165 objective_sm_pd = 0.0_dp
1166
1167 ! loop over domain - for applying power law later on
1168 do idomain = 1, domainmeta%nDomains
1169
1170 ! get domain information
1171 mask1 = level1(idomain)%mask
1172 ncols1 = level1(idomain)%ncols
1173 nrows1 = level1(idomain)%nrows
1174
1175 ! allocate
1176 allocate(mask_times(size(smoptisim(idomain)%dataSim, dim = 2)))
1177 allocate(pd_time_series(size(smoptisim(idomain)%dataSim, dim = 2)))
1178 allocate(mat1(nrows1, ncols1))
1179 allocate(mat2(nrows1, ncols1))
1180 allocate(mask_sm(nrows1, ncols1))
1181
1182 ! initalize
1183 mask_times = .false.
1184 pd_time_series = 0.0_dp
1185
1186 ! calculate pattern similarity criterion
1187 do itime = 1, size(smoptisim(idomain)%dataSim, dim = 2)
1188 mat1 = unpack(l1_smobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
1189 mat2 = unpack(smoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
1190 mask_sm = unpack(l1_smobs(idomain)%maskObs(:, itime), mask1, .false.)
1191 pd_time_series = pd(mat1, mat2, mask = mask_sm, valid = mask_times(itime))
1192 end do
1193
1194 if (count(mask_times) > 0_i4) then
1195 ! calculate avergae PD over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1196 objective_sm_pd = objective_sm_pd + &
1197 ((1.0_dp - sum(pd_time_series, mask = mask_times) / real(count(mask_times), dp)) / &
1198 real(domainmeta%overallnumberofdomains, dp))**6
1199 else
1200 call error_message('***ERROR: mo_objective_funtion: objective_sm_pd: No soil moisture observations available!')
1201 end if
1202
1203 ! deallocate
1204 deallocate(mask_times)
1205 deallocate(pd_time_series)
1206 deallocate(mat1)
1207 deallocate(mat2)
1208 deallocate(mask_sm)
1209 call smoptisim(idomain)%destroy()
1210 end do
1211 deallocate(smoptisim)
1212
1213#ifndef MPI
1214 objective_sm_pd = objective_sm_pd**onesixth
1215
1216 call message(' objective_sm_pd = ', num2str(objective_sm_pd, '(F9.5)'))
1217#endif
1218
1219 END FUNCTION objective_sm_pd
1220
1221 ! ------------------------------------------------------------------
1222
1223 ! NAME
1224 ! objective_sm_sse_standard_score
1225
1226 ! PURPOSE
1227 !> \brief Objective function for soil moisture.
1228
1229 !> \details The objective function only depends on a parameter vector.
1230 !> The model will be called with that parameter vector and
1231 !> the model output is subsequently compared to observed data.
1232
1233 !> Therefore the sum of squared errors (SSE) of the standard score of observed and
1234 !> modeled soil moisture is calculated. The standard score or normalization (anomaly)
1235 !> make the objctive function bias insensitive and basically the dynamics of the soil moisture
1236 !> is tried to capture by this obejective function.
1237 !> \f[ phi_i = \sum_{j=1}^K \{ standard\_score( SM_{obs}(j) )- standard\_score(SM_{sim}(j)) \}^2 \f]
1238 !> where
1239 !> \f$ standard\_score \f$ = standard score function,
1240 !> \f$ SM_{obs} \f$ = observed soil moisture,
1241 !> \f$ SM_{sim} \f$ = simulated soil moisture.
1242 !> \f$ K \f$ = valid elements in study domain.
1243 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1244 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1245 !> \f[ OF = \sqrt[6]{\sum(\phi_{i}/N)^6 }. \f]
1246 !> The observed data L1_sm, L1_sm_mask are global in this module.
1247
1248 ! INTENT(IN)
1249 !> \param[in] "real(dp), dimension(:) :: parameterset"
1250 !> \param[in] "procedure(eval_interface) :: eval"
1251
1252 ! RETURN
1253 !> \return real(dp) :: objective_sm_sse_standard_score — objective function value
1254 !> (which will be e.g. minimized by an optimization routine like DDS)
1255
1256 ! HISTORY
1257 !> \authors Matthias Zink
1258
1259 !> \date March 2015
1260
1261 ! Modifications:
1262 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1263
1264 FUNCTION objective_sm_sse_standard_score(parameterset, eval)
1265
1266 use mo_optimization_types, only : optidata_sim
1267 use mo_common_variables, only : level1, domainmeta
1268 use mo_errormeasures, only : sse
1269 use mo_global_variables, only : l1_smobs
1270 use mo_standard_score, only : standard_score
1271 use mo_string_utils, only : num2str
1272
1273 implicit none
1274
1275 real(dp), dimension(:), intent(in) :: parameterset
1276
1277 procedure(eval_interface), INTENT(IN), POINTER :: eval
1278
1279 real(dp) :: objective_sm_sse_standard_score
1280
1281 ! domain loop counter
1282 integer(i4) :: idomain
1283
1284 ! cell loop counter
1285 integer(i4) :: icell
1286
1287 ! ncells1 of level 1
1288 integer(i4) :: ncells1
1289
1290 ! number of invalid cells in catchment
1291 real(dp) :: invalid_cells
1292
1293 ! domains wise objectives
1294 real(dp) :: objective_sm_sse_standard_score_domain
1295
1296#ifndef MPI
1297 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1298#endif
1299
1300 ! simulated soil moisture
1301 type(optidata_sim), dimension(:), allocatable :: smoptisim
1302
1303
1304 call eval(parameterset, smoptisim = smoptisim)
1305
1306 ! initialize some variables
1307 objective_sm_sse_standard_score = 0.0_dp
1308
1309 ! loop over domain - for applying power law later on
1310 do idomain = 1, domainmeta%nDomains
1311
1312 ! init
1313 objective_sm_sse_standard_score_domain = 0.0_dp
1314 ! get domain information
1315 ncells1 = level1(idomain)%nCells
1316
1317 invalid_cells = 0.0_dp
1318 ! standard_score signal is calculated on individual grid cells
1319 do icell = 1, size(l1_smobs(idomain)%maskObs(:, :), dim = 1)
1320
1321 ! check for enough data points in time for statistical calculations (e.g. mean, stddev)
1322 if (count(l1_smobs(idomain)%maskObs(icell, :)) .LE. (0.10_dp * real(size(l1_smobs(idomain)%dataObs, dim = 2), dp))) then
1323 invalid_cells = invalid_cells + 1.0_dp
1324 cycle
1325 end if
1326 objective_sm_sse_standard_score_domain = objective_sm_sse_standard_score_domain + &
1327 sse(standard_score(l1_smobs(idomain)%dataObs(icell, :), mask = l1_smobs(idomain)%maskObs(icell, :)), &
1328 standard_score(smoptisim(idomain)%dataSim(icell, :), mask = l1_smobs(idomain)%maskObs(icell, :)), &
1329 mask = l1_smobs(idomain)%maskObs(icell, :))
1330
1331 end do
1332
1333 ! user information about invalid cells
1334 if (invalid_cells .GT. 0.5_dp) then
1335 call message(.LT.' WARNING: objective_sm: Detected invalid cells in study area ( 10 valid data points).')
1336 call message(' Fraction of invalid cells: ', &
1337 num2str(invalid_cells / real(ncells1, dp), '(F4.2)'))
1338 end if
1339
1340 ! calculate average soil moisture correlation over all domains with power law
1341 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1342 objective_sm_sse_standard_score = objective_sm_sse_standard_score + &
1343 (objective_sm_sse_standard_score_domain / real(domainmeta%overallNumberOfDomains, dp))**6
1344 end do
1345
1346#ifndef MPI
1347 objective_sm_sse_standard_score = objective_sm_sse_standard_score**onesixth
1348
1349 call message(' objective_sm_sse_standard_score = ', num2str(objective_sm_sse_standard_score, '(E12.5)'))
1350#endif
1351
1352 END FUNCTION objective_sm_sse_standard_score
1353
1354
1355 ! -----------------------------------------------------------------
1356
1357 ! NAME
1358 ! objective_kge_q_rmse_tws
1359
1360 ! PURPOSE
1361 !> \brief Objective function of KGE for runoff and RMSE for domain_avg TWS (standarized scores)
1362
1363 !> \details Objective function of KGE for runoff and RMSE for domain_avg TWS (standarized scores)
1364
1365 ! INTENT(IN)
1366 !> \param[in] "real(dp), dimension(:) :: parameterset"
1367 !> \param[in] "procedure(eval_interface) :: eval"
1368
1369 ! RETURN
1370 !> \return real(dp) :: objective_kge_q_rmse_tws — objective function value
1371 !> (which will be e.g. minimized by an optimization routine like DDS)
1372
1373 ! HISTORY
1374 !> \authors Oldrich Rakovec, Rohini Kumar
1375
1376 !> \date Oct. 2015
1377
1378 ! Modifications:
1379 ! Stephan Thober Oct 2015 - moved tws optimization from mo_mrm_objective_function_runoff here
1380 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1381 ! Maren Kaluza Oct 2019 - changed averaging function for tws, this will not produce the same output as before
1382
1383 FUNCTION objective_kge_q_rmse_tws(parameterset, eval)
1384
1385 use mo_optimization_types, only : optidata_sim
1386 use mo_common_constants, only : eps_dp, nodata_dp
1387 use mo_common_mhm_mrm_variables, only : evalper
1388 use mo_common_variables, only : domainmeta
1389 use mo_global_variables, only : l1_twsaobs
1390 use mo_errormeasures, only : rmse
1391 use mo_julian, only : caldat
1392 use mo_moment, only : mean
1393 use mo_standard_score, only : classified_standard_score
1394 use mo_string_utils, only : num2str
1395 use mo_temporal_aggregation, only : day2mon_average
1396 use mo_errormeasures, only : kge
1397 use mo_mrm_objective_function_runoff, only : extract_runoff
1398
1399 implicit none
1400
1401 real(dp), dimension(:), intent(in) :: parameterset
1402
1403 procedure(eval_interface), INTENT(IN), POINTER :: eval
1404
1405 real(dp) :: objective_kge_q_rmse_tws
1406
1407 ! modelled runoff for a given parameter set
1408 ! dim1=nTimeSteps, dim2=nGauges
1409 real(dp), allocatable, dimension(:, :) :: runoff
1410
1411 !> simulated tws
1412 type(optidata_sim), dimension(:), allocatable :: twsoptisim
1413
1414 ! domain counter, month counters
1415 integer(i4) :: domainid, idomain, pp, mmm
1416
1417 integer(i4) :: year, month, day
1418
1419 real(dp), dimension(domainMeta%nDomains) :: inittime
1420
1421 ! simulated tws
1422 real(dp), dimension(:), allocatable :: tws_catch_avg_domain
1423
1424 ! measured tws
1425 real(dp), dimension(:), allocatable :: tws_opti_catch_avg_domain
1426
1427 ! mask for measured tws
1428 logical, dimension(:), allocatable :: tws_obs_mask
1429
1430 ! total number of months
1431 integer(i4) :: nmonths
1432
1433 ! vector with months' classes
1434 integer(i4), dimension(:), allocatable :: month_classes
1435
1436 ! monthly values anomaly time series
1437 real(dp), DIMENSION(:), allocatable :: tws_sim_m_anom, tws_obs_m_anom
1438
1439 ! rmse_tws(domainMeta%nDomains)
1440 real(dp), dimension(:), allocatable :: rmse_tws
1441
1442 ! obj. functions
1443 real(dp) :: rmse_tws_avg, kge_q_avg
1444
1445 integer(i4) :: ngaugestotal
1446
1447 ! aggregated simulated runoff
1448 real(dp), dimension(:), allocatable :: runoff_agg
1449
1450 ! measured runoff
1451 real(dp), dimension(:), allocatable :: runoff_obs
1452
1453 ! mask for measured runoff
1454 logical, dimension(:), allocatable :: runoff_obs_mask
1455
1456 ! kge_q(nGaugesTotal)
1457 real(dp), dimension(:), allocatable :: kge_q
1458
1459 ! gauges counter
1460 integer(i4) :: gg
1461
1462 ! obtain hourly values of runoff and tws:
1463 allocate(twsoptisim(domainmeta%nDomains))
1464 call eval(parameterset, runoff = runoff, twsoptisim = twsoptisim)
1465
1466 !--------------------------------------------
1467 !! TWS
1468 !--------------------------------------------
1469
1470 ! allocate
1471 allocate(rmse_tws(domainmeta%nDomains))
1472 rmse_tws(:) = nodata_dp
1473
1474 do idomain = 1, domainmeta%nDomains
1475 if (.not. (l1_twsaobs(idomain)%timeStepInput == -2)) then
1476 call message('objective_kge_q_rmse_tws: current implementation of this subroutine only allows monthly timesteps')
1477 end if
1478 domainid = domainmeta%indices(idomain)
1479
1480 ! extract tws the same way as runoff using mrm
1481 ! Note that with the change from tws(iDomain, tt) to tws(tt, :) this
1482 ! will not work like before and also does maybe not make sense
1483 call create_domain_avg_tws(idomain, twsoptisim, tws_catch_avg_domain, tws_opti_catch_avg_domain, tws_obs_mask)
1484
1485 ! check for potentially 2 years of data
1486 if (count(tws_obs_mask) .lt. 12 * 2) then
1487 call message('objective_kge_q_rmse_tws: Length of TWS data of domain ', trim(adjustl(num2str(domainid))), &
1488 ' less than 2 years: this is not recommended')
1489 end if
1490
1491 ! get initial time of the evaluation period
1492 inittime(idomain) = real(evalper(idomain)%julStart, dp)
1493
1494 ! get calendar days, months, year
1495 call caldat(int(inittime(idomain)), yy = year, mm = month, dd = day)
1496
1497 nmonths = size(tws_obs_mask)
1498
1499 allocate (month_classes(nmonths))
1500 allocate (tws_obs_m_anom(nmonths))
1501 allocate (tws_sim_m_anom(nmonths))
1502
1503 month_classes(:) = 0
1504 tws_obs_m_anom(:) = nodata_dp
1505 tws_sim_m_anom(:) = nodata_dp
1506
1507 ! define months' classes
1508 mmm = month
1509 do pp = 1, nmonths
1510 month_classes(pp) = mmm
1511 if (mmm .LT. 12) then
1512 mmm = mmm + 1
1513 else
1514 mmm = 1
1515 end if
1516 end do
1517
1518 ! calculate standard score
1519 tws_obs_m_anom = classified_standard_score(tws_opti_catch_avg_domain, month_classes, mask = tws_obs_mask)
1520 tws_sim_m_anom = classified_standard_score(tws_catch_avg_domain, month_classes, mask = tws_obs_mask)
1521 rmse_tws(idomain) = rmse(tws_sim_m_anom, tws_obs_m_anom, mask = tws_obs_mask)
1522
1523 deallocate (month_classes)
1524 deallocate (tws_sim_m_anom)
1525 deallocate (tws_obs_m_anom)
1526 deallocate (tws_catch_avg_domain, tws_opti_catch_avg_domain, tws_obs_mask)
1527
1528 call twsoptisim(idomain)%destroy()
1529 end do
1530
1531 rmse_tws_avg = sum(rmse_tws(:), abs(rmse_tws - nodata_dp) .gt. eps_dp) / &
1532 real(count(abs(rmse_tws - nodata_dp) .gt. eps_dp), dp)
1533 deallocate(rmse_tws)
1534 deallocate(twsoptisim)
1535
1536 !--------------------------------------------
1537 !! RUNOFF
1538 !--------------------------------------------
1539 kge_q_avg = 0_dp
1540 ngaugestotal = size(runoff, dim = 2)
1541 allocate(kge_q(ngaugestotal))
1542 kge_q(:) = nodata_dp
1543
1544 do gg = 1, ngaugestotal
1545
1546 ! extract runoff
1547 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1548
1549 ! check for potentially 2 years of data
1550 pp = count(runoff_agg .ge. 0.0_dp)
1551 if (pp .lt. 365 * 2) then
1552 call message('objective_kge_q_rmse_tws: The simulation at gauge ', trim(adjustl(num2str(gg))), &
1553 ' is not long enough. Please provide at least 730 days of data.')
1554 end if
1555 ! calculate KGE for each domain:
1556 kge_q(gg) = kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1557 deallocate (runoff_agg, runoff_obs, runoff_obs_mask)
1558
1559 end do
1560
1561 ! calculate average KGE value for runoff
1562 kge_q_avg = sum(kge_q(:), abs(kge_q - nodata_dp) .gt. eps_dp) / &
1563 real(count(abs(kge_q - nodata_dp) .gt. eps_dp), dp)
1564 deallocate(kge_q)
1565
1566 !
1567 objective_kge_q_rmse_tws = rmse_tws_avg * (1._dp - kge_q_avg)
1568
1569 call message(' objective_kge_q_rmse_tws = ', num2str(objective_kge_q_rmse_tws, '(F9.5)'))
1570
1571 END FUNCTION objective_kge_q_rmse_tws
1572
1573 ! ------------------------------------------------------------------
1574
1575 ! NAME
1576 ! objective_neutrons_kge_catchment_avg
1577
1578 ! PURPOSE
1579 !> \brief Objective function for neutrons.
1580
1581 !> \details The objective function only depends on a parameter vector.
1582 !> The model will be called with that parameter vector and
1583 !> the model output is subsequently compared to observed data.
1584
1585 !> Therefore, the Kling-Gupta model efficiency \f$ KGE \f$ of the catchment average
1586 !> neutrons (N) is calculated
1587 !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f]
1588 !> where
1589 !> \f$ r \f$ = Pearson product-moment CORRELATION coefficient
1590 !> \f$ \alpha \f$ = ratio of simulated mean to observed mean SM
1591 !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation
1592 !> is calculated and the objective function for a given domain \f$ i \f$ is
1593 !> \f[ \phi_{i} = 1.0 - KGE_{i} \f]
1594 !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods.
1595 !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.
1596
1597 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1598 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1599 !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f]
1600 !> The observed data L1_neutronsdata, L1_neutronsdata_mask are global in this module.
1601
1602 ! INTENT(IN)
1603 !> \param[in] "real(dp), dimension(:) :: parameterset"
1604 !> \param[in] "procedure(eval_interface) :: eval"
1605
1606 ! RETURN
1607 !> \return real(dp) :: objective_neutrons_kge_catchment_avg — objective function value
1608 !> (which will be e.g. minimized by an optimization routine)
1609
1610 ! HISTORY
1611 !> \authors Martin Schroen
1612
1613 !> \date Jun 2015
1614
1615 ! Modifications:
1616 ! Maren Kaluza Mar 2018 - changed format string to '(I10)'
1617 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1618
1619 FUNCTION objective_neutrons_kge_catchment_avg(parameterset, eval)
1620
1621 use mo_optimization_types, only : optidata_sim
1622 use mo_common_constants, only : nodata_dp
1623 use mo_common_variables, only : domainmeta
1624 use mo_errormeasures, only : kge
1625 use mo_global_variables, only : l1_neutronsobs
1626 use mo_moment, only : average
1627 use mo_string_utils, only : num2str
1628
1629 implicit none
1630
1631 real(dp), dimension(:), intent(in) :: parameterset
1632
1633 procedure(eval_interface), INTENT(IN), POINTER :: eval
1634
1635 real(dp) :: objective_neutrons_kge_catchment_avg
1636
1637 ! domain loop counter
1638 integer(i4) :: idomain
1639
1640 ! time loop counter
1641 integer(i4) :: itime
1642
1643 ! for sixth root
1644#ifndef MPI
1645 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1646#endif
1647
1648 ! spatial average of observed neutrons
1649 real(dp), dimension(:), allocatable :: neutrons_catch_avg_domain
1650
1651 ! spatial avergae of modeled neutrons
1652 real(dp), dimension(:), allocatable :: neutrons_opti_catch_avg_domain
1653
1654 ! simulated neutrons
1655 type(optidata_sim), dimension(:), allocatable :: neutronsoptisim
1656
1657 ! mask for valid neutrons catchment avg time steps
1658 logical, dimension(:), allocatable :: mask_times
1659
1660
1661 allocate(neutronsoptisim(domainmeta%nDomains))
1662 call eval(parameterset, neutronsoptisim = neutronsoptisim)
1663
1664 ! initialize some variables
1665 objective_neutrons_kge_catchment_avg = 0.0_dp
1666
1667 ! loop over domain - for applying power law later on
1668 do idomain = 1, domainmeta%nDomains
1669
1670 ! allocate
1671 allocate(mask_times(size(neutronsoptisim(idomain)%dataSim, dim = 2)))
1672 allocate(neutrons_catch_avg_domain(size(neutronsoptisim(idomain)%dataSim, dim = 2)))
1673 allocate(neutrons_opti_catch_avg_domain(size(neutronsoptisim(idomain)%dataSim, dim = 2)))
1674
1675 ! initalize
1676 mask_times = .true.
1677 neutrons_catch_avg_domain = nodata_dp
1678 neutrons_opti_catch_avg_domain = nodata_dp
1679
1680 ! calculate catchment average soil moisture
1681 do itime = 1, size(neutronsoptisim(idomain)%dataSim, dim = 2)
1682
1683 ! check for enough data points in time for correlation
1684 if (all(.NOT. l1_neutronsobs(idomain)%maskObs(:, itime))) then
1685 call message('WARNING: neutrons data at time ', num2str(itime, '(I10)'), ' is empty.')
1686 !call message('WARNING: objective_neutrons_kge_catchment_avg: ignored current time step since less than')
1687 !call message(' 10 valid cells available in soil moisture observation')
1688 mask_times(itime) = .false.
1689 cycle
1690 end if
1691 neutrons_catch_avg_domain(itime) = average(l1_neutronsobs(idomain)%dataObs(:, itime), &
1692 mask = l1_neutronsobs(idomain)%maskObs(:, itime))
1693 neutrons_opti_catch_avg_domain(itime) = average(neutronsoptisim(idomain)%dataSim(:, itime), &
1694 mask = l1_neutronsobs(idomain)%maskObs(:, itime))
1695 end do
1696
1697 ! calculate average neutrons KGE over all domains with power law
1698 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1699 objective_neutrons_kge_catchment_avg = objective_neutrons_kge_catchment_avg + &
1700 ((1.0_dp - kge(neutrons_catch_avg_domain, neutrons_opti_catch_avg_domain, mask = mask_times)) / &
1701 real(domainmeta%overallnumberofdomains, dp))**6
1702
1703 ! deallocate
1704 deallocate(mask_times)
1705 deallocate(neutrons_catch_avg_domain)
1706 deallocate(neutrons_opti_catch_avg_domain)
1707
1708 call neutronsoptisim(idomain)%destroy()
1709 end do
1710 deallocate(neutronsoptisim)
1711
1712#ifndef MPI
1713 objective_neutrons_kge_catchment_avg = objective_neutrons_kge_catchment_avg**onesixth
1714
1715 call message(' objective_neutrons_kge_catchment_avg = ', num2str(objective_neutrons_kge_catchment_avg, '(F9.5)'))
1716#endif
1717
1718 END FUNCTION objective_neutrons_kge_catchment_avg
1719
1720 ! ------------------------------------------------------------------
1721
1722 ! NAME
1723 ! objective_et_kge_catchment_avg
1724
1725 ! PURPOSE
1726 !> \brief Objective function for evpotranspirstion (et).
1727
1728 !> \details The objective function only depends on a parameter vector.
1729 !> The model will be called with that parameter vector and
1730 !> the model output is subsequently compared to observed data.
1731
1732 !> Therefore, the Kling-Gupta model efficiency \f$ KGE \f$ of the catchment average
1733 !> evapotranspiration (et) is calculated
1734 !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f]
1735 !> where
1736 !> \f$ r \f$ = Pearson product-moment correlation coefficient
1737 !> \f$ \alpha \f$ = ratio of simulated mean to observed mean SM
1738 !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation
1739 !> is calculated and the objective function for a given domain \f$ i \f$ is
1740 !> \f[ \phi_{i} = 1.0 - KGE_{i} \f]
1741 !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods.
1742 !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.
1743
1744 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1745 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1746 !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f]
1747 !> The observed data L1_et, L1_et_mask are global in this module.
1748
1749 ! INTENT(IN)
1750 !> \param[in] "real(dp), dimension(:) :: parameterset"
1751 !> \param[in] "procedure(eval_interface) :: eval"
1752
1753 ! RETURN
1754 !> \return real(dp) :: objective_et_kge_catchment_avg — objective function value
1755 !> (which will be e.g. minimized by an optimization routine)
1756
1757 ! HISTORY
1758 !> \authors Johannes Brenner
1759
1760 !> \date Feb 2017
1761
1762 ! Modifications:
1763 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1764
1765 FUNCTION objective_et_kge_catchment_avg(parameterset, eval)
1766
1767 use mo_optimization_types, only : optidata_sim
1768 use mo_common_constants, only : nodata_dp
1769 use mo_common_variables, only : domainmeta
1770 use mo_errormeasures, only : kge
1771 use mo_moment, only : average
1772 use mo_string_utils, only : num2str
1773
1774 implicit none
1775
1776 real(dp), dimension(:), intent(in) :: parameterset
1777
1778 procedure(eval_interface), INTENT(IN), POINTER :: eval
1779
1780 real(dp) :: objective_et_kge_catchment_avg
1781
1782 ! domain loop counter
1783 integer(i4) ::idomain
1784
1785 ! for sixth root
1786#ifndef MPI
1787 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1788#endif
1789
1790 ! spatial average of observed et
1791 real(dp), dimension(:), allocatable :: et_catch_avg_domain
1792
1793 ! spatial avergae of modeled et
1794 real(dp), dimension(:), allocatable :: et_opti_catch_avg_domain
1795
1796 !> simulated et
1797 type(optidata_sim), dimension(:), allocatable :: etoptisim
1798
1799 ! mask for valid et catchment avg time steps
1800 logical, dimension(:), allocatable :: mask_times
1801
1802
1803 call eval(parameterset, etoptisim = etoptisim)
1804
1805 ! initialize some variables
1806 objective_et_kge_catchment_avg = 0.0_dp
1807
1808 ! loop over domain - for applying power law later on
1809 allocate(etoptisim(domainmeta%nDomains))
1810 do idomain = 1, domainmeta%nDomains
1811
1812 ! create et array input
1813 call create_domain_avg_et(idomain, etoptisim, et_catch_avg_domain, &
1814 et_opti_catch_avg_domain, mask_times)
1815 ! calculate average ET KGE over all domains with power law
1816 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1817
1818 objective_et_kge_catchment_avg = objective_et_kge_catchment_avg + &
1819 ((1.0_dp - kge(et_catch_avg_domain, et_opti_catch_avg_domain, mask = mask_times)) / &
1820 real(domainmeta%overallnumberofdomains, dp))**6
1821
1822 ! deallocate
1823 deallocate(mask_times)
1824 deallocate(et_catch_avg_domain)
1825 deallocate(et_opti_catch_avg_domain)
1826 call etoptisim(idomain)%destroy()
1827 end do
1828 deallocate(etoptisim)
1829
1830#ifndef MPI
1831 objective_et_kge_catchment_avg = objective_et_kge_catchment_avg**onesixth
1832
1833 call message(' objective_et_kge_catchment_avg = ', num2str(objective_et_kge_catchment_avg, '(F9.5)'))
1834#endif
1835
1836 END FUNCTION objective_et_kge_catchment_avg
1837
1838 ! -----------------------------------------------------------------
1839
1840 ! NAME
1841 ! objective_kge_q_sm_corr
1842
1843 ! PURPOSE
1844 !> \brief Objective function of KGE for runoff and correlation for SM
1845
1846 !> \details Objective function of KGE for runoff and SSE for soil moisture (standarized scores).
1847 !> Further details can be found in the documentation of objective functions
1848 !> '14 - objective_multiple_gauges_kge_power6' and '13 - objective_sm_corr'.
1849
1850 ! INTENT(IN)
1851 !> \param[in] "real(dp), dimension(:) :: parameterset"
1852 !> \param[in] "procedure(eval_interface) :: eval"
1853
1854 ! RETURN
1855 !> \return real(dp) :: objective_kge_q_sse_sm — objective function value
1856 !> (which will be e.g. minimized by an optimization routine like DDS)
1857
1858 ! HISTORY
1859 !> \authors Matthias Zink
1860
1861 !> \date Mar. 2017
1862
1863 ! Modifications:
1864 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1865
1866 FUNCTION objective_kge_q_sm_corr(parameterset, eval)
1867
1868 use mo_optimization_types, only : optidata_sim
1869 use mo_common_variables, only : level1, domainmeta
1870 use mo_global_variables, only : l1_smobs
1871 use mo_moment, only : correlation
1872 use mo_string_utils, only : num2str
1873 use mo_errormeasures, only : kge
1874 use mo_mrm_objective_function_runoff, only : extract_runoff
1875
1876 implicit none
1877
1878 real(dp), dimension(:), intent(in) :: parameterset
1879
1880 procedure(eval_interface), INTENT(IN), POINTER :: eval
1881
1882 real(dp) :: objective_kge_q_sm_corr
1883
1884 real(dp) :: objective_sm
1885
1886 real(dp) :: objective_kge
1887
1888 ! number of invalid cells in catchment
1889 real(dp) :: invalid_cells
1890
1891 ! modelled runoff for a given parameter set
1892 ! dim1=nTimeSteps, dim2=nGauges
1893 real(dp), allocatable, dimension(:, :) :: runoff
1894
1895 ! domain loop counter
1896 integer(i4) :: idomain
1897
1898 ! cell loop counter
1899 integer(i4) :: icell
1900
1901 ! ncells1 of level 1
1902 integer(i4) :: ncells1
1903
1904 ! domains wise objectives
1905 real(dp) :: objective_sm_domain
1906
1907 ! simulated soil moisture
1908 type(optidata_sim), dimension(:), allocatable :: smoptisim
1909
1910 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1911
1912 ! gauges counter
1913 integer(i4) :: gg
1914
1915 integer(i4) :: ngaugestotal
1916
1917 ! aggregated simulated runoff
1918 real(dp), dimension(:), allocatable :: runoff_agg
1919
1920 ! measured runoff
1921 real(dp), dimension(:), allocatable :: runoff_obs
1922
1923 ! mask for measured runoff
1924 logical, dimension(:), allocatable :: runoff_obs_mask
1925
1926 ! run mHM
1927 allocate(smoptisim(domainmeta%nDomains))
1928 call eval(parameterset, runoff = runoff, smoptisim = smoptisim)
1929
1930 ! -----------------------------
1931 ! SOIL MOISTURE
1932 ! -----------------------------
1933
1934 ! initialize some variables
1935 objective_sm = 0.0_dp
1936
1937 ! loop over domain - for applying power law later on
1938 do idomain = 1, domainmeta%nDomains
1939
1940 ! init
1941 objective_sm_domain = 0.0_dp
1942 ! get domain information
1943 ncells1 = level1(idomain)%nCells
1944
1945 ! correlation signal is calculated on individual grid cells
1946 invalid_cells = 0.0_dp
1947 do icell = 1, size(l1_smobs(idomain)%maskObs(:, :), dim = 1)
1948
1949 ! check for enough data points in time for statistical calculations (e.g. mean, stddev)
1950 if (count(l1_smobs(idomain)%maskObs(icell, :)) .LE. (0.10_dp * real(size(l1_smobs(idomain)%dataObs, dim = 2), dp))) then
1951 invalid_cells = invalid_cells + 1.0_dp
1952 cycle
1953 end if
1954
1955 ! calculate ojective function
1956 objective_sm_domain = objective_sm_domain + &
1957 correlation(l1_smobs(idomain)%dataObs(icell, :), smoptisim(idomain)%dataSim(icell, :), &
1958 mask = l1_smobs(idomain)%maskObs(icell, :))
1959 end do
1960
1961 ! user information about invalid cells
1962 if (invalid_cells .GT. 0.5_dp) then
1963 call message(.LT.' WARNING: objective_sm: Detected invalid cells in study area ( 10 valid data points).')
1964 call message(' Fraction of invalid cells: ', &
1965 num2str(invalid_cells / real(ncells1, dp), '(F4.2)'))
1966 end if
1967
1968 ! calculate average soil moisture objective over all domains with power law
1969 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1970 objective_sm = objective_sm + &
1971 ((1.0_dp - objective_sm_domain / real(ncells1, dp)) / real(domainmeta%overallNumberOfDomains, dp))**6
1972 call smoptisim(idomain)%destroy()
1973 end do
1974 deallocate(smoptisim)
1975
1976 ! compromise solution - sixth root
1977 objective_sm = objective_sm**onesixth
1978
1979 ! -----------------------------
1980 ! RUNOFF
1981 ! -----------------------------
1982 objective_kge = 0.0_dp
1983 ngaugestotal = size(runoff, dim = 2)
1984
1985 do gg = 1, ngaugestotal
1986
1987 ! extract runoff
1988 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1989
1990 ! KGE
1991 objective_kge = objective_kge + &
1992 ((1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)) / real(ngaugestotal, dp))**6
1993
1994 end do
1995
1996 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
1997
1998 ! compromise solution - sixth root
1999 objective_kge = objective_kge**onesixth
2000
2001 ! equal weighted compromise objective functions for discharge and soilmoisture
2002 ! ToDo: why do we use the sixth root of of objective_sm and objective_kge
2003 ! only to take the power to 6 here again, when we never need the
2004 ! intermediate results?
2005#ifdef MPI
2006 objective_kge_q_sm_corr = (objective_sm**6 + objective_kge**6)
2007#else
2008 objective_kge_q_sm_corr = (objective_sm**6 + objective_kge**6)**onesixth
2009
2010 call message(' objective_kge_q_sm_corr = ', num2str(objective_kge_q_sm_corr, '(F9.5)'))
2011#endif
2012 ! print*, "1-SM 2-Q : ", 1.0_dp-objective_sm, 1.0_dp-objective_kge ! MZMZMZMZ
2013
2014 END FUNCTION objective_kge_q_sm_corr
2015
2016
2017 ! -----------------------------------------------------------------
2018
2019 ! NAME
2020 ! objective_kge_q_et
2021
2022 ! PURPOSE
2023 !> \brief Objective function of KGE for runoff and KGE for ET
2024
2025 !> \details Objective function of KGE for runoff and KGE for ET.
2026 !> Further details can be found in the documentation of objective functions
2027 !> '14 - objective_multiple_gauges_kge_power6'.
2028
2029 ! INTENT(IN)
2030 !> \param[in] "real(dp), dimension(:) :: parameterset"
2031 !> \param[in] "procedure(eval_interface) :: eval"
2032
2033 ! RETURN
2034 !> \return real(dp) :: objective_kge_q_et — objective function value
2035 !> (which will be e.g. minimized by an optimization routine like DDS)
2036
2037 ! HISTORY
2038 !> \authors Johannes Brenner
2039
2040 !> \date July 2017
2041
2042 ! Modifications:
2043 ! Robert Schweppe Jun 2018 - refactoring and reformatting
2044
2045 FUNCTION objective_kge_q_et(parameterset, eval)
2046
2047 use mo_optimization_types, only : optidata_sim
2048 use mo_common_variables, only : level1, domainmeta
2049 use mo_errormeasures, only : kge
2050 use mo_global_variables, only : l1_etobs
2051 use mo_string_utils, only : num2str
2052 use mo_mrm_objective_function_runoff, only : extract_runoff
2053
2054 implicit none
2055
2056 real(dp), dimension(:), intent(in) :: parameterset
2057
2058 procedure(eval_interface), INTENT(IN), POINTER :: eval
2059
2060 real(dp) :: objective_kge_q_et
2061
2062 real(dp) :: objective_et
2063
2064 real(dp) :: objective_q
2065
2066 ! number of invalid cells in catchment
2067 real(dp) :: invalid_cells
2068
2069 ! modelled runoff for a given parameter set
2070 ! dim1=nTimeSteps, dim2=nGauges
2071 real(dp), allocatable, dimension(:, :) :: runoff
2072
2073 ! domain loop counter
2074 integer(i4) :: idomain
2075
2076 ! cell loop counter
2077 integer(i4) :: icell
2078
2079 ! ncells1 of level 1
2080 integer(i4) :: ncells1
2081
2082 ! domains wise objectives
2083 real(dp) :: objective_et_domain
2084
2085 !> simulated et
2086 type(optidata_sim), dimension(:), allocatable :: etoptisim
2087
2088 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2089
2090 ! gauges counter
2091 integer(i4) :: gg
2092
2093 integer(i4) :: ngaugestotal
2094
2095 ! aggregated simulated runoff
2096 real(dp), dimension(:), allocatable :: runoff_agg
2097
2098 ! measured runoff
2099 real(dp), dimension(:), allocatable :: runoff_obs
2100
2101 ! mask for measured runoff
2102 logical, dimension(:), allocatable :: runoff_obs_mask
2103
2104 ! run mHM
2105 allocate(etoptisim(domainmeta%nDomains))
2106 call eval(parameterset, runoff = runoff, etoptisim = etoptisim)
2107
2108 ! -----------------------------
2109 ! EVAPOTRANSPIRATION
2110 ! -----------------------------
2111
2112 ! initialize some variables
2113 objective_et = 0.0_dp
2114
2115 ! loop over domain - for applying power law later on
2116 do idomain = 1, domainmeta%nDomains
2117
2118 ! init
2119 objective_et_domain = 0.0_dp
2120 ! get domain information
2121 ncells1 = level1(idomain)%nCells
2122
2123 ! correlation signal is calculated on individual grid cells
2124 invalid_cells = 0.0_dp
2125 do icell = 1, size(l1_etobs(idomain)%maskObs, dim = 1)
2126
2127 ! check for enough data points in time for statistical calculations (e.g. mean, stddev)
2128 if (count(l1_etobs(idomain)%maskObs(icell, :)) .LE. &
2129 (0.10_dp * real(size(l1_etobs(idomain)%dataObs(:, :), dim = 2), dp))) then
2130 invalid_cells = invalid_cells + 1.0_dp
2131 cycle
2132 end if
2133
2134 ! calculate ojective function
2135 objective_et_domain = objective_et_domain + &
2136 kge(l1_etobs(idomain)%dataObs(icell, :), etoptisim(idomain)%dataSim(icell, :), &
2137 mask = l1_etobs(idomain)%maskObs(icell, :))
2138 end do
2139
2140 ! user information about invalid cells
2141 if (invalid_cells .GT. 0.5_dp) then
2142 call message(.LT.' WARNING: objective_et: Detected invalid cells in study area ( 10 valid data points).')
2143 call message(' Fraction of invalid cells: ', &
2144 num2str(invalid_cells / real(ncells1, dp), '(F4.2)'))
2145 end if
2146
2147 ! calculate average soil moisture objective over all domains with power law
2148 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
2149 objective_et = objective_et + &
2150 ((1.0_dp - objective_et_domain / real(ncells1, dp)) / real(domainmeta%overallNumberOfDomains, dp))**6
2151 call etoptisim(idomain)%destroy()
2152 end do
2153 deallocate(etoptisim)
2154
2155 ! compromise solution - sixth root
2156 objective_et = objective_et**onesixth
2157
2158 ! -----------------------------
2159 ! RUNOFF
2160 ! -----------------------------
2161 objective_q = 0.0_dp
2162 ngaugestotal = size(runoff, dim = 2)
2163
2164 do gg = 1, ngaugestotal
2165
2166 ! extract runoff
2167 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2168
2169 ! KGE
2170 objective_q = objective_q + &
2171 ((1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)) / real(ngaugestotal, dp))**6
2172
2173 end do
2174
2175 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
2176
2177 ! compromise solution - sixth root
2178 objective_q = objective_q**onesixth
2179
2180 ! equal weighted compromise objective functions for discharge and soilmoisture
2181 ! ToDo: why do we use the sixth root of of objective_sm and objective_kge
2182 ! only to take the power to 6 here again, when we never need the
2183 ! intermediate results?
2184#ifdef MPI
2185 objective_kge_q_et = (objective_et**6 + objective_q**6)
2186#else
2187 objective_kge_q_et = (objective_et**6 + objective_q**6)**onesixth
2188
2189 call message(' objective_kge_q_et = ', num2str(objective_kge_q_et, '(F9.5)'))
2190#endif
2191 ! print*, "1-SM 2-Q : ", 1.0_dp-objective_sm, 1.0_dp-objective_kge ! MZMZMZMZ
2192
2193 END FUNCTION objective_kge_q_et
2194
2195
2196 ! -----------------------------------------------------------------
2197
2198 ! NAME
2199 ! objective_kge_q_BFI
2200
2201 ! PURPOSE
2202 !> \brief Objective function of KGE for runoff and BFI absulute difference
2203
2204 !> \details Objective function of KGE for runoff and KGE for ET.
2205 !> Further details can be found in the documentation of objective functions
2206 !> '14 - objective_multiple_gauges_kge_power6'.
2207
2208 ! INTENT(IN)
2209 !> \param[in] "real(dp), dimension(:) :: parameterset"
2210 !> \param[in] "procedure(eval_interface) :: eval"
2211
2212 ! RETURN
2213 !> \return real(dp) :: objective_kge_q_BFI — objective function value
2214 !> (which will be e.g. minimized by an optimization routine like DDS)
2215
2216 ! HISTORY
2217 !> \authors Sebastian Müller
2218
2219 !> \date Apr 2022
2220
2221 FUNCTION objective_kge_q_bfi(parameterset, eval)
2222
2223 use mo_optimization_types, only : optidata_sim
2224 use mo_common_variables, only : level1, domainmeta
2225 use mo_errormeasures, only : kge
2226 use mo_global_variables, only : bfi_obs
2227 use mo_string_utils, only : num2str
2228 use mo_mrm_objective_function_runoff, only : extract_runoff
2229
2230 implicit none
2231
2232 real(dp), dimension(:), intent(in) :: parameterset
2233
2234 procedure(eval_interface), INTENT(IN), POINTER :: eval
2235
2236 real(dp) :: objective_kge_q_bfi
2237 real(dp) :: objective_bfi
2238 real(dp) :: objective_q
2239
2240 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2241
2242 ! number of invalid cells in catchment
2243 real(dp) :: invalid_cells
2244
2245 ! modelled runoff for a given parameter set
2246 ! dim1=nTimeSteps, dim2=nGauges
2247 real(dp), allocatable, dimension(:, :) :: runoff
2248
2249 ! domain loop counter
2250 integer(i4) :: idomain
2251
2252 !> baseflow index for each domain
2253 real(dp), dimension(:), allocatable :: bfi
2254
2255 ! counter
2256 integer(i4) :: gg, i
2257
2258 integer(i4) :: ngaugestotal
2259
2260 ! aggregated simulated runoff
2261 integer(i4), dimension(:), allocatable :: domain_ids, domain_ids_pack
2262
2263 ! aggregated simulated runoff
2264 real(dp), dimension(:), allocatable :: runoff_agg
2265
2266 ! measured runoff
2267 real(dp), dimension(:), allocatable :: runoff_obs
2268
2269 ! mask for measured runoff
2270 logical, dimension(:), allocatable :: runoff_obs_mask
2271
2272 ! run mHM
2273 allocate(bfi(domainmeta%nDomains))
2274 call eval(parameterset, runoff = runoff, bfi = bfi)
2275
2276 ! -----------------------------
2277 ! BFI
2278 ! -----------------------------
2279
2280 ! initialize some variables
2281 objective_bfi = 0.0_dp
2282
2283 if ( any(bfi_obs < 0.0_dp) ) then
2284 allocate(domain_ids(domainmeta%nDomains))
2285 allocate(domain_ids_pack(count(bfi_obs < 0.0_dp)))
2286 domain_ids = [(i, i=1,size(domain_ids))]
2287 domain_ids_pack = pack(domain_ids, mask=(bfi_obs < 0.0_dp))
2288 call error_message( &
2289 "objective_kge_q_BFI: missing BFI values for domain ", &
2290 trim(adjustl(num2str(domain_ids_pack(1)))) &
2291 )
2292 end if
2293
2294 ! loop over domain - for applying power law later on
2295 do idomain = 1, domainmeta%nDomains
2296 objective_bfi = objective_bfi + abs(bfi(idomain) - bfi_obs(idomain)) / domainmeta%nDomains
2297 end do
2298
2299 ! -----------------------------
2300 ! RUNOFF
2301 ! -----------------------------
2302 objective_q = 0.0_dp
2303 ngaugestotal = size(runoff, dim = 2)
2304
2305 do gg = 1, ngaugestotal
2306
2307 ! extract runoff
2308 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2309
2310 ! KGE
2311 objective_q = objective_q + &
2312 ((1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)) / real(ngaugestotal, dp))**6
2313
2314 end do
2315
2316 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
2317
2318 ! compromise solution - sixth root
2319 objective_q = objective_q**onesixth
2320
2321 objective_kge_q_bfi = (objective_bfi + 1._dp)*objective_q
2322 call message(' objective_kge_q_BFI = ', num2str(objective_kge_q_bfi, '(F9.5)'))
2323
2324 END FUNCTION objective_kge_q_bfi
2325
2326
2327 ! -----------------------------------------------------------------
2328
2329 ! NAME
2330 ! objective_kge_q_rmse_et
2331
2332 ! PURPOSE
2333 !> \brief Objective function of KGE for runoff and RMSE for domain_avg ET (standarized scores)
2334
2335 !> \details Objective function of KGE for runoff and RMSE for domain_avg ET (standarized scores)
2336
2337 ! INTENT(IN)
2338 !> \param[in] "real(dp), dimension(:) :: parameterset"
2339 !> \param[in] "procedure(eval_interface) :: eval"
2340
2341 ! RETURN
2342 !> \return real(dp) :: objective_kge_q_rmse_et &mdash; objective function value
2343 !> (which will be e.g. minimized by an optimization routine like DDS)
2344
2345 ! HISTORY
2346 !> \authors Johannes Brenner
2347
2348 !> \date July 2017
2349
2350 ! Modifications:
2351 ! Robert Schweppe Jun 2018 - refactoring and reformatting
2352
2353 FUNCTION objective_kge_q_rmse_et(parameterset, eval)
2354
2355 use mo_optimization_types, only : optidata_sim
2356 use mo_common_constants, only : eps_dp, nodata_dp
2357 use mo_common_mhm_mrm_variables, only : evalper
2358 use mo_common_variables, only : domainmeta
2359 use mo_errormeasures, only : rmse
2360 use mo_global_variables, only : l1_etobs
2361 use mo_julian, only : caldat
2362 use mo_moment, only : average, mean
2363 use mo_standard_score, only : classified_standard_score
2364 use mo_string_utils, only : num2str
2365 use mo_temporal_aggregation, only : day2mon_average
2366 use mo_errormeasures, only : kge
2367 use mo_mrm_objective_function_runoff, only : extract_runoff
2368
2369 implicit none
2370
2371 real(dp), dimension(:), intent(in) :: parameterset
2372
2373 procedure(eval_interface), INTENT(IN), POINTER :: eval
2374
2375 real(dp) :: objective_kge_q_rmse_et
2376
2377 ! modelled runoff for a given parameter set
2378 ! dim1=nTimeSteps, dim2=nGauges
2379 real(dp), allocatable, dimension(:, :) :: runoff
2380
2381 !> simulated et
2382 type(optidata_sim), dimension(:), allocatable :: etoptisim
2383
2384 ! time loop counter
2385 integer(i4) :: itime
2386
2387 ! domain counter, month counters
2388 integer(i4) :: idomain, pp, mmm
2389
2390 integer(i4) :: year, month, day
2391
2392 real(dp), dimension(domainMeta%nDomains) :: inittime
2393
2394 ! total number of months
2395 integer(i4) :: nmonths
2396
2397 ! vector with months' classes
2398 integer(i4), dimension(:), allocatable :: month_classes
2399
2400 ! monthly values original time series
2401 real(dp), dimension(:), allocatable :: et_sim_m, et_obs_m
2402
2403 ! monthly values anomaly time series
2404 real(dp), dimension(:), allocatable :: et_sim_m_anom, et_obs_m_anom
2405
2406 logical, dimension(:), allocatable :: et_obs_m_mask
2407
2408 ! kge_q(nGaugesTotal)
2409 real(dp), dimension(:), allocatable :: rmse_et
2410
2411 ! obj. functions
2412 real(dp) :: rmse_et_avg, kge_q_avg
2413
2414 ! spatial average of observed et
2415 real(dp), dimension(:), allocatable :: et_catch_avg_domain
2416
2417 ! spatial avergae of modeled et
2418 real(dp), dimension(:), allocatable :: et_opti_catch_avg_domain
2419
2420 ! mask for valid et catchment avg time steps
2421 logical, dimension(:), allocatable :: mask_times
2422
2423 ! rmse_et(domainMeta%nDomains)
2424 real(dp), dimension(:), allocatable :: kge_q
2425
2426 ! gauges counter
2427 integer(i4) :: gg
2428
2429 integer(i4) :: ngaugestotal
2430
2431 ! aggregated simulated runoff
2432 real(dp), dimension(:), allocatable :: runoff_agg
2433
2434 ! measured runoff
2435 real(dp), dimension(:), allocatable :: runoff_obs
2436
2437 ! mask for measured runoff
2438 logical, dimension(:), allocatable :: runoff_obs_mask
2439
2440 ! obtain simulation values of runoff (hourly) and ET
2441 ! for ET only valid cells (domains concatenated)
2442 ! et_opti: aggregate ET to needed time step for optimization (see timeStep_et_input)
2443 allocate(etoptisim(domainmeta%nDomains))
2444 call eval(parameterset, runoff = runoff, etoptisim = etoptisim)
2445
2446 !--------------------------------------------
2447 !! EVAPOTRANSPIRATION
2448 !--------------------------------------------
2449
2450 ! allocate
2451 allocate(rmse_et(domainmeta%nDomains))
2452 rmse_et(:) = nodata_dp
2453
2454 do idomain = 1, domainmeta%nDomains
2455
2456 ! allocate
2457 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
2458 allocate(et_catch_avg_domain(size(etoptisim(idomain)%dataSim, dim = 2)))
2459 allocate(et_opti_catch_avg_domain(size(etoptisim(idomain)%dataSim, dim = 2)))
2460
2461 ! initalize
2462 mask_times = .true.
2463 et_catch_avg_domain = nodata_dp
2464 et_opti_catch_avg_domain = nodata_dp
2465
2466 ! calculate catchment average evapotranspiration
2467 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
2468 ! check for enough data points in time for correlation
2469 if (all(.NOT. l1_etobs(idomain)%maskObs(:, itime))) then
2470 !write (*,*) 'WARNING: et data at time ', iTime, ' is empty.'
2471 !call message('WARNING: objective_et_kge_catchment_avg: ignored current time step since less than')
2472 !call message(' 10 valid cells available in evapotranspiration observation')
2473 mask_times(itime) = .false.
2474 cycle
2475 end if
2476 ! spatial average of observed ET
2477 et_catch_avg_domain(itime) = average(l1_etobs(idomain)%dataObs(:, itime), &
2478 mask = l1_etobs(idomain)%maskObs(:, itime))
2479 ! spatial avergae of modeled ET
2480 et_opti_catch_avg_domain(itime) = average(etoptisim(idomain)%dataSim(:, itime), &
2481 mask = l1_etobs(idomain)%maskObs(:, itime))
2482 end do
2483
2484 ! get initial time of the evaluation period
2485 inittime(idomain) = real(evalper(idomain)%julStart, dp)
2486
2487 ! get calendar days, months, year
2488 call caldat(int(inittime(idomain)), yy = year, mm = month, dd = day)
2489
2490 ! if evapotranspiration input daily
2491 select case(l1_etobs(idomain)%timeStepInput)
2492 ! JBJBJB
2493 ! daily: aggregate to monthly mean
2494 case(-1)
2495 ! calculate monthly averages from daily values of the model
2496 call day2mon_average(et_opti_catch_avg_domain, year, month, day, et_sim_m, misval = nodata_dp)
2497 ! calculate monthly averages from daily values of the observations
2498 call day2mon_average(et_catch_avg_domain, year, month, day, et_obs_m, misval = nodata_dp)
2499 !
2500 ! monthly: proceed without action
2501 case(-2)
2502 ! simulation
2503 allocate(et_sim_m(size(etoptisim(idomain)%dataSim, dim = 2)))
2504 et_sim_m = et_opti_catch_avg_domain
2505 ! observation
2506 allocate(et_obs_m(size(etoptisim(idomain)%dataSim, dim = 2)))
2507 et_obs_m = et_catch_avg_domain
2508
2509 ! yearly: ERROR stop program
2510 case(-3)
2511 call error_message('***ERROR: objective_kge_q_rmse_et: time step of evapotranspiration yearly.')
2512 end select
2513 ! remove mean from modelled time series
2514 et_sim_m(:) = et_sim_m(:) - mean(et_sim_m(:))
2515 ! remove mean from observed time series
2516 et_obs_m(:) = et_obs_m(:) - mean(et_obs_m(:))
2517 ! get number of months for given domain
2518 nmonths = size(et_obs_m)
2519 allocate (month_classes(nmonths))
2520 allocate (et_obs_m_mask(nmonths))
2521 allocate (et_obs_m_anom(nmonths))
2522 allocate (et_sim_m_anom(nmonths))
2523
2524 month_classes(:) = 0
2525 et_obs_m_mask(:) = .true.
2526 et_obs_m_anom(:) = nodata_dp
2527 et_sim_m_anom(:) = nodata_dp
2528
2529 ! define months' classes
2530 mmm = month
2531 do pp = 1, nmonths
2532 month_classes(pp) = mmm
2533 if (mmm .LT. 12) then
2534 mmm = mmm + 1
2535 else
2536 mmm = 1
2537 end if
2538 end do
2539 ! double check missing data
2540 ! define mask for missing data in observations (there are always data for simulations)
2541 where(abs(et_obs_m - nodata_dp) .lt. eps_dp) et_obs_m_mask = .false.
2542
2543 ! calculate standard score
2544 et_obs_m_anom = classified_standard_score(et_obs_m, month_classes, mask = et_obs_m_mask)
2545 et_sim_m_anom = classified_standard_score(et_sim_m, month_classes, mask = et_obs_m_mask)
2546 rmse_et(idomain) = rmse(et_sim_m_anom, et_obs_m_anom, mask = et_obs_m_mask)
2547
2548 deallocate (month_classes)
2549 deallocate (et_obs_m)
2550 deallocate (et_sim_m)
2551 deallocate (et_obs_m_mask)
2552 deallocate (et_sim_m_anom)
2553 deallocate (et_obs_m_anom)
2554 !end if
2555
2556 call etoptisim(idomain)%destroy()
2557 end do
2558
2559 rmse_et_avg = sum(rmse_et(:), abs(rmse_et - nodata_dp) .gt. eps_dp) / &
2560 real(count(abs(rmse_et - nodata_dp) .gt. eps_dp), dp)
2561 deallocate(rmse_et)
2562 deallocate(etoptisim)
2563
2564 !--------------------------------------------
2565 !! RUNOFF
2566 !--------------------------------------------
2567 kge_q_avg = 0_dp
2568 ngaugestotal = size(runoff, dim = 2)
2569 allocate(kge_q(ngaugestotal))
2570 kge_q(:) = nodata_dp
2571
2572 do gg = 1, ngaugestotal
2573
2574 ! extract runoff
2575 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2576
2577 ! check for whether to proceed with this domain or not
2578 ! potentially 3 years of data
2579 !pp = count(runoff_agg .ge. 0.0_dp )
2580 !if (pp .lt. 365*3 ) then
2581 ! deallocate (runoff_agg, runoff_obs, runoff_obs_mask)
2582 ! cycle
2583 ! else
2584 ! calculate KGE for each domain:
2585 kge_q(gg) = kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
2586 deallocate (runoff_agg, runoff_obs, runoff_obs_mask)
2587 ! end if
2588
2589 end do
2590
2591 ! calculate average KGE value for runoff
2592 kge_q_avg = sum(kge_q(:), abs(kge_q - nodata_dp) .gt. eps_dp) / &
2593 real(count(abs(kge_q - nodata_dp) .gt. eps_dp), dp)
2594 deallocate(kge_q)
2595
2596 !
2597 objective_kge_q_rmse_et = rmse_et_avg * (1._dp - kge_q_avg)
2598
2599 call message(' objective_kge_q_rmse_et = ', num2str(objective_kge_q_rmse_et, '(F9.5)'))
2600
2601 END FUNCTION objective_kge_q_rmse_et
2602
2603 subroutine create_domain_avg_tws(iDomain, twsOptiSim, tws_catch_avg_domain, &
2604 tws_opti_catch_avg_domain, mask_times)
2605 use mo_optimization_types, only : optidata_sim
2606 use mo_common_constants, only : nodata_dp
2607 use mo_global_variables, only : l1_twsaobs
2608 use mo_moment, only : average
2609 ! current domain Id
2610 integer(i4), intent(in) :: iDomain
2611
2612 ! simulated tws
2613 type(optidata_sim), dimension(:), intent(in) :: twsOptiSim
2614
2615 ! aggregated simulated
2616 real(dp), dimension(:), allocatable, intent(out) :: tws_catch_avg_domain
2617
2618 ! extracted measured
2619 real(dp), dimension(:), allocatable, intent(out) :: tws_opti_catch_avg_domain
2620
2621 ! mask of no data values
2622 logical, dimension(:), allocatable, intent(out) :: mask_times
2623
2624 ! local
2625 ! time loop counter
2626 integer(i4) :: iTime
2627
2628 ! allocate
2629 allocate(mask_times(size(twsoptisim(idomain)%dataSim, dim = 2)))
2630 allocate(tws_catch_avg_domain(size(twsoptisim(idomain)%dataSim, dim = 2)))
2631 allocate(tws_opti_catch_avg_domain(size(twsoptisim(idomain)%dataSim, dim = 2)))
2632
2633 ! initalize
2634 mask_times = .true.
2635 tws_catch_avg_domain = nodata_dp
2636 tws_opti_catch_avg_domain = nodata_dp
2637
2638 ! calculate catchment average evapotranspiration
2639 do itime = 1, size(twsoptisim(idomain)%dataSim, dim = 2)
2640
2641 ! check for enough data points in time for correlation
2642 if (all(.NOT. l1_twsaobs(idomain)%maskObs(:, itime))) then
2643 !write (*,*) 'WARNING: et data at time ', iTime, ' is empty.'
2644 !call message('WARNING: objective_et_kge_catchment_avg: ignored current time step since less than')
2645 !call message(' 10 valid cells available in evapotranspiration observation')
2646 mask_times(itime) = .false.
2647 cycle
2648 end if
2649
2650 tws_catch_avg_domain(itime) = average(l1_twsaobs(idomain)%dataObs(:, itime), &
2651 mask = l1_twsaobs(idomain)%maskObs(:, itime))
2652 tws_opti_catch_avg_domain(itime) = average(twsoptisim(idomain)%dataSim(:, itime), &
2653 mask = l1_twsaobs(idomain)%maskObs(:, itime))
2654 end do
2655
2656 end subroutine create_domain_avg_tws
2657
2658 subroutine create_domain_avg_et(iDomain, etOptiSim, et_catch_avg_domain, &
2659 et_opti_catch_avg_domain, mask_times)
2660 use mo_optimization_types, only : optidata_sim
2661 use mo_common_constants, only : nodata_dp
2662 use mo_global_variables, only : l1_etobs
2663 use mo_moment, only : average
2664 ! current domain Id
2665 integer(i4), intent(in) :: iDomain
2666
2667 ! simulated et
2668 type(optidata_sim), dimension(:), intent(in) :: etOptiSim
2669
2670 ! aggregated simulated
2671 real(dp), dimension(:), allocatable, intent(out) :: et_catch_avg_domain
2672
2673 ! extracted measured
2674 real(dp), dimension(:), allocatable, intent(out) :: et_opti_catch_avg_domain
2675
2676 ! mask of no data values
2677 logical, dimension(:), allocatable, intent(out) :: mask_times
2678
2679 ! local
2680 ! time loop counter
2681 integer(i4) :: iTime
2682
2683 ! allocate
2684 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
2685 allocate(et_catch_avg_domain(size(etoptisim(idomain)%dataSim, dim = 2)))
2686 allocate(et_opti_catch_avg_domain(size(etoptisim(idomain)%dataSim, dim = 2)))
2687
2688 ! initalize
2689 mask_times = .true.
2690 et_catch_avg_domain = nodata_dp
2691 et_opti_catch_avg_domain = nodata_dp
2692
2693 ! calculate catchment average evapotranspiration
2694 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
2695
2696 ! check for enough data points in time for correlation
2697 if (all(.NOT. l1_etobs(idomain)%maskObs(:, itime))) then
2698 !write (*,*) 'WARNING: et data at time ', iTime, ' is empty.'
2699 !call message('WARNING: objective_et_kge_catchment_avg: ignored current time step since less than')
2700 !call message(' 10 valid cells available in evapotranspiration observation')
2701 mask_times(itime) = .false.
2702 cycle
2703 end if
2704
2705 et_catch_avg_domain(itime) = average(l1_etobs(idomain)%dataObs(:, itime), &
2706 mask = l1_etobs(idomain)%maskObs(:, itime))
2707 et_opti_catch_avg_domain(itime) = average(etoptisim(idomain)%dataSim(:, itime), &
2708 mask = l1_etobs(idomain)%maskObs(:, itime))
2709 end do
2710
2711 end subroutine create_domain_avg_et
2712
2713 subroutine convert_tws_to_twsa(twsOptiSim, L1_twsaObs, twsaOptiSim)
2714 use mo_optimization_types, only : optidata_sim, optidata
2715 use mo_moment, only : average
2716 ! simulated tws
2717 type(optidata_sim), intent(in) :: twsOptiSim
2718 ! observed twsa
2719 type(optidata), intent(in) :: L1_twsaObs
2720 ! simulated twsa
2721 type(optidata_sim), intent(inout) :: twsaOptiSim
2722
2723 ! local
2724 integer(i4) :: iCell
2725 real(dp) :: twsa_av_cell
2726
2727 allocate(twsaoptisim%dataSim(size(twsoptisim%dataSim(:, :), dim = 1), size(twsoptisim%dataSim(:, :), dim = 2)))
2728
2729 do icell = 1, size(twsoptisim%dataSim(:, :), dim = 1)
2730 twsa_av_cell = average(twsoptisim%dataSim(icell, :), mask = l1_twsaobs%maskObs(icell, :))
2731 twsaoptisim%dataSim(icell, :) = twsoptisim%dataSim(icell, :) - twsa_av_cell
2732 end do
2733
2734 end subroutine convert_tws_to_twsa
2735
2736END MODULE mo_objective_function
mo_optimization_utils::eval_interface
Interface for evaluation routine.
Definition mo_optimization_utils.f90:18
mo_common_constants
Provides constants commonly used by mHM, mRM and MPR.
Definition mo_common_constants.f90:15
mo_common_constants::eps_dp
real(dp), parameter, public eps_dp
epsilon(1.0) in double precision
Definition mo_common_constants.f90:25
mo_common_constants::nodata_dp
real(dp), parameter, public nodata_dp
Definition mo_common_constants.f90:31
mo_common_mhm_mrm_variables
Provides structures needed by mHM, mRM and/or mpr.
Definition mo_common_mHM_mRM_variables.f90:19
mo_common_mhm_mrm_variables::opti_function
integer(i4), public opti_function
Definition mo_common_mHM_mRM_variables.f90:67
mo_common_mhm_mrm_variables::evalper
type(period), dimension(:), allocatable, public evalper
Definition mo_common_mHM_mRM_variables.f90:43
mo_common_mpi_tools
tools for MPI communication that are mHM or mRM specific
Definition mo_common_MPI_tools.F90:14
mo_common_types
Provides common types needed by mHM, mRM and/or mpr.
Definition mo_common_types.F90:22
mo_common_variables
Provides structures needed by mHM, mRM and/or mpr.
Definition mo_common_variables.F90:27
mo_common_variables::domainmeta
type(domain_meta), public domainmeta
Definition mo_common_variables.F90:99
mo_common_variables::level1
type(grid), dimension(:), allocatable, target, public level1
Definition mo_common_variables.F90:74
mo_global_variables
Main global variables for mHM.
Definition mo_global_variables.f90:71
mo_global_variables::l1_twsaobs
type(optidata), dimension(:), allocatable, public l1_twsaobs
this stores L1_tws, the mask, the directory of the observerd data, and the timestepInput of the simul...
Definition mo_global_variables.f90:118
mo_global_variables::l1_neutronsobs
type(optidata), dimension(:), allocatable, public l1_neutronsobs
Definition mo_global_variables.f90:114
mo_global_variables::l1_smobs
type(optidata), dimension(:), allocatable, public l1_smobs
Definition mo_global_variables.f90:112
mo_global_variables::l1_etobs
type(optidata), dimension(:), allocatable, public l1_etobs
Definition mo_global_variables.f90:116
mo_global_variables::bfi_obs
real(dp), dimension(:), allocatable, public bfi_obs
given base-flow index per domain
Definition mo_global_variables.f90:123
mo_mrm_objective_function_runoff
Objective Functions for Optimization of mHM/mRM against runoff.
Definition mo_mrm_objective_function_runoff.F90:55
mo_mrm_objective_function_runoff::extract_runoff
subroutine, public extract_runoff(gaugeid, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
extracts runoff data from global variables
Definition mo_mrm_objective_function_runoff.F90:2479
mo_mrm_objective_function_runoff::objective_kge
real(dp) function objective_kge(parameterset, eval)
Objective function of KGE.
Definition mo_mrm_objective_function_runoff.F90:2098
mo_objective_function
Objective Functions for Optimization of mHM.
Definition mo_objective_function.F90:31
mo_objective_function::objective_sm_corr
real(dp) function objective_sm_corr(parameterset, eval)
Objective function for soil moisture.
Definition mo_objective_function.F90:978
mo_objective_function::objective
real(dp) function, public objective(parameterset, eval, arg1, arg2, arg3)
Wrapper for objective functions.
Definition mo_objective_function.F90:92
mo_objective_function::create_domain_avg_et
subroutine create_domain_avg_et(idomain, etoptisim, et_catch_avg_domain, et_opti_catch_avg_domain, mask_times)
Definition mo_objective_function.F90:2660
mo_objective_function::objective_neutrons_kge_catchment_avg
real(dp) function objective_neutrons_kge_catchment_avg(parameterset, eval)
Objective function for neutrons.
Definition mo_objective_function.F90:1620
mo_objective_function::objective_sm_pd
real(dp) function objective_sm_pd(parameterset, eval)
Objective function for soil moisture.
Definition mo_objective_function.F90:1111
mo_objective_function::convert_tws_to_twsa
subroutine convert_tws_to_twsa(twsoptisim, l1_twsaobs, twsaoptisim)
Definition mo_objective_function.F90:2714
mo_objective_function::objective_kge_q_rmse_et
real(dp) function objective_kge_q_rmse_et(parameterset, eval)
Objective function of KGE for runoff and RMSE for domain_avg ET (standarized scores)
Definition mo_objective_function.F90:2354
mo_objective_function::objective_kge_q_sm_corr
real(dp) function objective_kge_q_sm_corr(parameterset, eval)
Objective function of KGE for runoff and correlation for SM.
Definition mo_objective_function.F90:1867
mo_objective_function::create_domain_avg_tws
subroutine create_domain_avg_tws(idomain, twsoptisim, tws_catch_avg_domain, tws_opti_catch_avg_domain, mask_times)
Definition mo_objective_function.F90:2605
mo_objective_function::objective_sm_sse_standard_score
real(dp) function objective_sm_sse_standard_score(parameterset, eval)
Objective function for soil moisture.
Definition mo_objective_function.F90:1265
mo_objective_function::objective_kge_q_rmse_tws
real(dp) function objective_kge_q_rmse_tws(parameterset, eval)
Objective function of KGE for runoff and RMSE for domain_avg TWS (standarized scores)
Definition mo_objective_function.F90:1384
mo_objective_function::objective_kge_q_bfi
real(dp) function objective_kge_q_bfi(parameterset, eval)
Objective function of KGE for runoff and BFI absulute difference.
Definition mo_objective_function.F90:2222
mo_objective_function::init_indexarray_for_opti_data
subroutine init_indexarray_for_opti_data(domainmeta, optidataoption, noptidomains, opti_domain_indices)
creates an index array of the inidices of the domains eval should MPI process.
Definition mo_objective_function.F90:901
mo_objective_function::objective_q_et_tws_kge_catchment_avg
real(dp) function, dimension(6) objective_q_et_tws_kge_catchment_avg(parameterset, eval)
Objective function for et, tws and q.
Definition mo_objective_function.F90:643
mo_objective_function::objective_sm_kge_catchment_avg
real(dp) function objective_sm_kge_catchment_avg(parameterset, eval)
Wrapper for objective functions.
Definition mo_objective_function.F90:499
mo_objective_function::objective_et_kge_catchment_avg
real(dp) function objective_et_kge_catchment_avg(parameterset, eval)
Objective function for evpotranspirstion (et).
Definition mo_objective_function.F90:1766
mo_objective_function::objective_kge_q_et
real(dp) function objective_kge_q_et(parameterset, eval)
Objective function of KGE for runoff and KGE for ET.
Definition mo_objective_function.F90:2046
mo_optimization_types
Type definitions for optimization routines.
Definition mo_optimization_types.f90:10
mo_optimization_utils
Utility functions, such as interface definitions, for optimization routines.
Definition mo_optimization_utils.f90:8
mo_common_types::domain_meta
DOMAIN general description.
Definition mo_common_types.F90:95
mo_optimization_types::optidata_sim
type for simulated optional data
Definition mo_optimization_types.f90:35
mo_optimization_types::optidata
optional data, such as sm, neutrons, et, tws
Definition mo_optimization_types.f90:26